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 ...
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)