Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_InterlacedOpUnitTest.cpp
Go to the documentation of this file.
1/*
2// @HEADER
3// ***********************************************************************
4//
5// Stokhos Package
6// Copyright (2009) Sandia Corporation
7//
8// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
9// license for use of this work by or on behalf of the U.S. Government.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
39//
40// ***********************************************************************
41// @HEADER
42*/
43
44#include <Teuchos_ConfigDefs.hpp>
45#include <Teuchos_UnitTestHarness.hpp>
46#include <Teuchos_TimeMonitor.hpp>
47#include <Teuchos_RCP.hpp>
48
50
51// Stokhos Stochastic Galerkin
52#include "Stokhos_Epetra.hpp"
55
56#ifdef HAVE_MPI
57#include "Epetra_MpiComm.h"
58#else
59#include "Epetra_SerialComm.h"
60#endif
61#include "Epetra_CrsGraph.h"
62#include "Epetra_Map.h"
63#include "EpetraExt_BlockUtility.h"
64#include "EpetraExt_RowMatrixOut.h"
65
66TEUCHOS_UNIT_TEST(interlaced_op, test)
67{
68#ifdef HAVE_MPI
69 Teuchos::RCP<const Epetra_Comm> comm = Teuchos::rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
70#else
71 Teuchos::RCP<const Epetra_Comm> comm = Teuchos::rcp(new Epetra_SerialComm);
72#endif
73
74 //int rank = comm->MyPID();
75 int numProc = comm->NumProc();
76
77 int num_KL = 1;
78 int porder = 5;
79 bool full_expansion = false;
80
81 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis = buildBasis(num_KL,porder);
82 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk;
83 Teuchos::RCP<Stokhos::ParallelData> sg_parallel_data;
84 Teuchos::RCP<Stokhos::OrthogPolyExpansion<int,double> > expansion;
85 {
86 if(full_expansion)
87 Cijk = basis->computeTripleProductTensor();
88 else
89 Cijk = basis->computeLinearTripleProductTensor();
90
91 Teuchos::ParameterList parallelParams;
92 parallelParams.set("Number of Spatial Processors", numProc);
93 sg_parallel_data = Teuchos::rcp(new Stokhos::ParallelData(basis, Cijk, comm,
94 parallelParams));
95
96 expansion = Teuchos::rcp(new Stokhos::AlgebraicOrthogPolyExpansion<int,double>(basis,
97 Cijk));
98 }
99 Teuchos::RCP<const EpetraExt::MultiComm> sg_comm = sg_parallel_data->getMultiComm();
100
101 // determinstic PDE graph
102 Teuchos::RCP<Epetra_Map> determRowMap = Teuchos::rcp(new Epetra_Map(-1,10,0,*comm));
103 Teuchos::RCP<Epetra_CrsGraph> determGraph = Teuchos::rcp(new Epetra_CrsGraph(Copy,*determRowMap,1));
104 for(int row=0;row<determRowMap->NumMyElements();row++) {
105 int gid = determRowMap->GID(row);
106 determGraph->InsertGlobalIndices(gid,1,&gid);
107 }
108 for(int row=1;row<determRowMap->NumMyElements()-1;row++) {
109 int gid = determRowMap->GID(row);
110 int indices[2] = {gid-1,gid+1};
111 determGraph->InsertGlobalIndices(gid,2,indices);
112 }
113 determGraph->FillComplete();
114
115 Teuchos::RCP<Teuchos::ParameterList> params = Teuchos::rcp(new Teuchos::ParameterList);
116 params->set("Scale Operator by Inverse Basis Norms", false);
117 params->set("Include Mean", true);
118 params->set("Only Use Linear Terms", false);
119
120 Teuchos::RCP<Stokhos::EpetraSparse3Tensor> epetraCijk =
121 Teuchos::rcp(new Stokhos::EpetraSparse3Tensor(basis,Cijk,sg_comm));
122 Teuchos::RCP<Stokhos::EpetraOperatorOrthogPoly> W_sg_blocks =
123 Teuchos::rcp(new Stokhos::EpetraOperatorOrthogPoly(basis, epetraCijk->getStochasticRowMap(), determRowMap, determRowMap, sg_comm));
124 for(int i=0; i<W_sg_blocks->size(); i++) {
125 Teuchos::RCP<Epetra_CrsMatrix> crsMat = Teuchos::rcp(new Epetra_CrsMatrix(Copy,*determGraph));
126 crsMat->PutScalar(1.0 + i);
127 W_sg_blocks->setCoeffPtr(i,crsMat); // allocate a bunch of matrices
128 }
129
130 Teuchos::RCP<const Epetra_Map> sg_map =
131 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(
132 *determRowMap, *(epetraCijk->getStochasticRowMap()),
133 *(epetraCijk->getMultiComm())));
134
135 // build an interlaced operator (object under test) and a benchmark
136 // fully assembled operator
138
139 Stokhos::InterlacedOperator op(sg_comm,basis,epetraCijk,determGraph,params);
140 op.PutScalar(0.0);
141 op.setupOperator(W_sg_blocks);
142
143 Stokhos::FullyAssembledOperator full_op(sg_comm,basis,epetraCijk,determGraph,sg_map,sg_map,params);
144 full_op.PutScalar(0.0);
145 full_op.setupOperator(W_sg_blocks);
146
147 // here we test interlaced operator against the fully assembled operator
149 bool result = true;
150 for(int i=0;i<100;i++) {
151 // build vector for fully assembled operator (blockwise)
152 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> x_vec_blocks =
153 Teuchos::rcp(new Stokhos::EpetraVectorOrthogPoly(basis,epetraCijk->getStochasticRowMap(),determRowMap,epetraCijk->getMultiComm()));
154 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> f_vec_blocks =
155 Teuchos::rcp(new Stokhos::EpetraVectorOrthogPoly(basis,epetraCijk->getStochasticRowMap(),determRowMap,epetraCijk->getMultiComm()));
156 Teuchos::RCP<Epetra_Vector> x_vec_blocked = x_vec_blocks->getBlockVector();
157 Teuchos::RCP<Epetra_Vector> f_vec_blocked = f_vec_blocks->getBlockVector();
158 x_vec_blocked->Random(); // build an initial vector
159 f_vec_blocked->PutScalar(0.0);
160
161 // build interlaced vectors
162 Teuchos::RCP<Epetra_Vector> x_vec_inter = Teuchos::rcp(new Epetra_Vector(op.OperatorDomainMap()));
163 Teuchos::RCP<Epetra_Vector> f_vec_inter = Teuchos::rcp(new Epetra_Vector(op.OperatorRangeMap()));
164 Teuchos::RCP<Epetra_Vector> f_vec_blk_inter = Teuchos::rcp(new Epetra_Vector(op.OperatorRangeMap()));
165 Stokhos::SGModelEvaluator_Interlaced::copyToInterlacedVector(*x_vec_blocks,*x_vec_inter); // copy random x to
166 f_vec_inter->PutScalar(0.0);
167
168 full_op.Apply(*x_vec_blocked,*f_vec_blocked);
169 op.Apply(*x_vec_inter,*f_vec_inter);
170
171 // copy blocked action to interlaced for comparison
173
174 // compute norm
175 double error = 0.0;
176 double true_norm = 0.0;
177 f_vec_blk_inter->NormInf(&true_norm);
178 f_vec_blk_inter->Update(-1.0,*f_vec_inter,1.0);
179 f_vec_blk_inter->NormInf(&error);
180
181 out << "rel error = " << error/true_norm << " ( " << true_norm << " ), ";
182 result &= (error/true_norm < 1e-14);
183 }
184 out << std::endl;
185
186 TEST_ASSERT(result);
187}
Copy
TEUCHOS_UNIT_TEST(interlaced_op, test)
Teuchos::RCP< const Stokhos::CompletePolynomialBasis< int, double > > buildBasis(int num_KL, int porder)
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const=0
virtual const Epetra_Map & OperatorDomainMap() const=0
virtual const Epetra_Map & OperatorRangeMap() const=0
Orthogonal polynomial expansions limited to algebraic operations.
A container class storing an orthogonal polynomial whose coefficients are vectors,...
A container class storing an orthogonal polynomial whose coefficients are vectors,...
An Epetra operator representing the block stochastic Galerkin operator generated by fully assembling ...
virtual int Apply(const Epetra_MultiVector &Input, Epetra_MultiVector &Result) const
Wrap Apply() to add a timer.
virtual void setupOperator(const Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > &poly)
Setup operator.
An Epetra operator representing the block stochastic Galerkin operator generated by fully assembling ...
virtual void setupOperator(const Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > &poly)
Setup operator.
static void copyToInterlacedVector(const Stokhos::EpetraVectorOrthogPoly &x_sg, Epetra_Vector &x)