Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
linear2d_diffusion.cpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Stokhos Package
5// Copyright (2009) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
38//
39// ***********************************************************************
40// @HEADER
41
42// Class implementing our problem
44
45// Epetra communicator
46#ifdef HAVE_MPI
47#include "Epetra_MpiComm.h"
48#else
49#include "Epetra_SerialComm.h"
50#endif
51
52// solver
53#include "Ifpack2_Factory.hpp"
54#include "BelosLinearProblem.hpp"
55#include "BelosTpetraAdapter.hpp"
56#include "BelosPseudoBlockCGSolMgr.hpp"
57#include "BelosPseudoBlockGmresSolMgr.hpp"
58#include "Tpetra_Core.hpp"
59#include "MatrixMarket_Tpetra.hpp"
60
61// Stokhos Stochastic Galerkin
62#include "Stokhos_Epetra.hpp"
63
64// Timing utilities
65#include "Teuchos_TimeMonitor.hpp"
66
67int main(int argc, char *argv[]) {
68 typedef double Scalar;
69 typedef double MeshScalar;
70 typedef double BasisScalar;
71 typedef int LocalOrdinal;
72 typedef int GlobalOrdinal;
73 typedef Tpetra::Map<>::node_type Node;
74
75 using Teuchos::RCP;
76 using Teuchos::rcp;
77 using Teuchos::ParameterList;
78
79 LocalOrdinal n = 32; // spatial discretization (per dimension)
80 LocalOrdinal num_KL = 2; // number of KL terms
81 LocalOrdinal p = 3; // polynomial order
82 BasisScalar mu = 0.2; // mean of exponential random field
83 BasisScalar s = 0.1; // std. dev. of exponential r.f.
84 bool nonlinear_expansion = false; // nonlinear expansion of diffusion coeff
85 // (e.g., log-normal)
86 bool symmetric = false; // use symmetric formulation
87
88// Initialize MPI
89#ifdef HAVE_MPI
90 MPI_Init(&argc,&argv);
91#endif
92
93 LocalOrdinal MyPID;
94
95 try {
96
97 {
98 TEUCHOS_FUNC_TIME_MONITOR("Total PCE Calculation Time");
99
100 // Create a communicator for Epetra objects
101 RCP<const Epetra_Comm> globalComm;
102#ifdef HAVE_MPI
103 globalComm = rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
104#else
105 globalComm = rcp(new Epetra_SerialComm);
106#endif
107 MyPID = globalComm->MyPID();
108
109 // Create Stochastic Galerkin basis and expansion
110 Teuchos::Array< RCP<const Stokhos::OneDOrthogPolyBasis<LocalOrdinal,BasisScalar> > > bases(num_KL);
111 for (LocalOrdinal i=0; i<num_KL; i++)
112 bases[i] = rcp(new Stokhos::LegendreBasis<LocalOrdinal,BasisScalar>(p,true));
113 RCP<const Stokhos::CompletePolynomialBasis<LocalOrdinal,BasisScalar> > basis =
115 1e-12));
116 LocalOrdinal sz = basis->size();
117 RCP<Stokhos::Sparse3Tensor<LocalOrdinal,BasisScalar> > Cijk;
118 if (nonlinear_expansion)
119 Cijk = basis->computeTripleProductTensor();
120 else
121 Cijk = basis->computeLinearTripleProductTensor();
122 RCP<Stokhos::OrthogPolyExpansion<LocalOrdinal,BasisScalar> > expansion =
124 Cijk));
125 if (MyPID == 0)
126 std::cout << "Stochastic Galerkin expansion size = " << sz << std::endl;
127
128 // Create stochastic parallel distribution
129 LocalOrdinal num_spatial_procs = -1;
130 if (argc > 1)
131 num_spatial_procs = std::atoi(argv[1]);
132 ParameterList parallelParams;
133 parallelParams.set("Number of Spatial Processors", num_spatial_procs);
134 // parallelParams.set("Rebalance Stochastic Graph", true);
135 // Teuchos::ParameterList& isorropia_params =
136 // parallelParams.sublist("Isorropia");
137 // isorropia_params.set("Balance objective", "nonzeros");
138 RCP<Stokhos::ParallelData> sg_parallel_data =
139 rcp(new Stokhos::ParallelData(basis, Cijk, globalComm,
140 parallelParams));
141 RCP<const EpetraExt::MultiComm> sg_comm =
142 sg_parallel_data->getMultiComm();
143 RCP<const Epetra_Comm> app_comm =
144 sg_parallel_data->getSpatialComm();
145
146 // Create Teuchos::Comm from Epetra_Comm
147 RCP< Teuchos::Comm<int> > teuchos_app_comm;
148#ifdef HAVE_MPI
149 RCP<const Epetra_MpiComm> app_mpi_comm =
150 Teuchos::rcp_dynamic_cast<const Epetra_MpiComm>(app_comm);
151 RCP<const Teuchos::OpaqueWrapper<MPI_Comm> > raw_mpi_comm =
152 Teuchos::opaqueWrapper(app_mpi_comm->Comm());
153 teuchos_app_comm = rcp(new Teuchos::MpiComm<int>(raw_mpi_comm));
154#else
155 teuchos_app_comm = rcp(new Teuchos::SerialComm<int>());
156#endif
157
158 // Create application
160 RCP<problem_type> model =
161 rcp(new problem_type(teuchos_app_comm, n, num_KL, s, mu,
162 nonlinear_expansion, symmetric));
163
164
165 // Create vectors and operators
166 typedef problem_type::Tpetra_Vector Tpetra_Vector;
167 typedef problem_type::Tpetra_CrsMatrix Tpetra_CrsMatrix;
168 typedef Tpetra::MatrixMarket::Writer<Tpetra_CrsMatrix> Writer;
169 RCP<const Tpetra_Vector> p = model->get_p_init(0);
170 RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(model->get_x_map());
171 x->putScalar(0.0);
172 RCP<Tpetra_Vector> f = Tpetra::createVector<Scalar>(model->get_f_map());
173 RCP<Tpetra_Vector> dx = Tpetra::createVector<Scalar>(model->get_x_map());
174 RCP<Tpetra_CrsMatrix> J = model->create_W();
175
176 // Create preconditioner
177 Teuchos::ParameterList precParams;
178 std::string prec_name = "RILUK";
179 precParams.set("fact: iluk level-of-fill", 1);
180 precParams.set("fact: iluk level-of-overlap", 0);
181 typedef Ifpack2::Preconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tprec;
182 Teuchos::RCP<Tprec> M;
183 Ifpack2::Factory factory;
184 M = factory.create<Tpetra_CrsMatrix>(prec_name, J);
185 M->setParameters(precParams);
186
187 // Evaluate model
188 model->computeResidual(*x, *p, *f);
189 model->computeJacobian(*x, *p, *J);
190 M->initialize();
191 M->compute();
192
193 // Print initial residual norm
194 Scalar norm_f;
195 f->norm2(Teuchos::arrayView(&norm_f,1));
196 if (MyPID == 0)
197 std::cout << "\nInitial residual norm = " << norm_f << std::endl;
198
199 // Setup Belos solver
200 RCP<ParameterList> belosParams = rcp(new ParameterList);
201 belosParams->set("Num Blocks", 20);
202 belosParams->set("Convergence Tolerance", 1e-12);
203 belosParams->set("Maximum Iterations", 1000);
204 belosParams->set("Verbosity", 33);
205 belosParams->set("Output Style", 1);
206 belosParams->set("Output Frequency", 1);
207 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
208 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
209 typedef Belos::OperatorTraits<Scalar,MV,OP> BOPT;
210 typedef Belos::MultiVecTraits<Scalar,MV> BMVT;
211 typedef Belos::LinearProblem<Scalar,MV,OP> BLinProb;
212 RCP< BLinProb > problem = rcp(new BLinProb(J, dx, f));
213 problem->setRightPrec(M);
214 problem->setProblem();
215 RCP<Belos::SolverManager<Scalar,MV,OP> > solver;
216 if (symmetric)
217 solver = rcp(new Belos::PseudoBlockCGSolMgr<Scalar,MV,OP>(problem,
218 belosParams));
219 else
220 solver = rcp(new Belos::PseudoBlockGmresSolMgr<Scalar,MV,OP>(problem,
221 belosParams));
222
223 // Solve linear system
224 Belos::ReturnType ret = solver->solve();
225
226 // Update x
227 x->update(-1.0, *dx, 1.0);
228
229 // Compute new residual & response function
230 RCP<Tpetra_Vector> g = Tpetra::createVector<Scalar>(model->get_g_map(0));
231 f->putScalar(0.0);
232 model->computeResidual(*x, *p, *f);
233 model->computeResponse(*x, *p, *g);
234
235 // Print final residual norm
236 f->norm2(Teuchos::arrayView(&norm_f,1));
237 if (MyPID == 0)
238 std::cout << "\nFinal residual norm = " << norm_f << std::endl;
239
240 // Print response
241 std::cout << "\nResponse = " << std::endl;
242 Writer::writeDense(std::cout, g);
243
244 }
245
246 Teuchos::TimeMonitor::summarize(std::cout);
247 Teuchos::TimeMonitor::zeroOutTimers();
248
249 }
250
251 catch (std::exception& e) {
252 std::cout << e.what() << std::endl;
253 }
254 catch (string& s) {
255 std::cout << s << std::endl;
256 }
257 catch (char *s) {
258 std::cout << s << std::endl;
259 }
260 catch (...) {
261 std::cout << "Caught unknown exception!" <<std:: endl;
262 }
263
264#ifdef HAVE_MPI
265 MPI_Finalize() ;
266#endif
267
268}
expr expr expr dx(i, j)
Orthogonal polynomial expansions limited to algebraic operations.
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Legendre polynomial basis.
A linear 2-D diffusion problem.
KokkosClassic::DefaultNode::DefaultNodeType Node
ScalarType g(const Teuchos::Array< ScalarType > &x, const ScalarType &y)
ScalarType f(const Teuchos::Array< ScalarType > &x, double a, double b)
int main(int argc, char *argv[])