Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
linear2d_diffusion_pce_ifpack2.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// Stokhos Stochastic Galerkin
43#include "Stokhos_Epetra.hpp"
44#include "Stokhos_Sacado.hpp"
45#include "Stokhos_Ifpack2.hpp"
46
47// Class implementing our problem
49
50// Epetra communicator
51#ifdef HAVE_MPI
52#include "Epetra_MpiComm.h"
53#else
54#include "Epetra_SerialComm.h"
55#endif
56
57// solver
58#include "Ifpack2_Factory.hpp"
59#include "BelosLinearProblem.hpp"
60#include "kokkos_pce_specializations.hpp"
61#include "BelosPseudoBlockCGSolMgr.hpp"
62#include "BelosPseudoBlockGmresSolMgr.hpp"
63#include "MatrixMarket_Tpetra.hpp"
64#include "BelosBlockGmresSolMgr.hpp"
65
66// Utilities
67#include "Teuchos_TimeMonitor.hpp"
68#include "Teuchos_CommandLineProcessor.hpp"
69
70// Scalar types
72
73// Random field types
75const int num_sg_rf = 2;
77const char *sg_rf_names[] = { "Uniform", "Log-Normal" };
78
79// Krylov methods
81const int num_krylov_method = 2;
83const char *krylov_method_names[] = { "GMRES", "CG" };
84
85// Preconditioning approaches
87const int num_sg_prec = 3;
89const char *sg_prec_names[] = { "None",
90 "Mean-Based",
91 "Stochastic" };
92
93// Stochastic division approaches
95const int num_sg_div = 5;
97const char *sg_div_names[] = { "Direct",
98 "SPD-Direct",
99 "Mean-Based",
100 "Quadrature",
101 "CG"};
102
103// Stochastic division preconditioner approaches
105const int num_sg_divprec = 5;
107const char *sg_divprec_names[] = { "None",
108 "Diag",
109 "Jacobi",
110 "GS",
111 "Schur"};
112
113
114// Option for Schur complement precond: full or diag D
116const int num_schur_option = 2;
118const char *schur_option_names[] = { "full", "diag"};
119
120// Full matrix or linear matrix (pb = dim + 1 ) used for preconditioner
122const int num_prec_option = 2;
124const char *prec_option_names[] = { "full", "linear"};
125
126int main(int argc, char *argv[]) {
127 typedef double MeshScalar;
128 typedef double BasisScalar;
129 typedef Tpetra::Map<>::node_type Node;
130 typedef Teuchos::ScalarTraits<Scalar>::magnitudeType magnitudeType;
131
132 using Teuchos::RCP;
133 using Teuchos::rcp;
134 using Teuchos::Array;
135 using Teuchos::ArrayRCP;
136 using Teuchos::ArrayView;
137 using Teuchos::ParameterList;
138
139// Initialize MPI
140#ifdef HAVE_MPI
141 MPI_Init(&argc,&argv);
142#endif
143
144 LocalOrdinal MyPID;
145
146 try {
147
148 // Create a communicator for Epetra objects
149 RCP<const Epetra_Comm> globalComm;
150#ifdef HAVE_MPI
151 globalComm = rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
152#else
153 globalComm = rcp(new Epetra_SerialComm);
154#endif
155 MyPID = globalComm->MyPID();
156
157 // Setup command line options
158 Teuchos::CommandLineProcessor CLP;
159 CLP.setDocString(
160 "This example runs an interlaced stochastic Galerkin solvers.\n");
161
162 int n = 32;
163 CLP.setOption("num_mesh", &n, "Number of mesh points in each direction");
164
165 bool symmetric = false;
166 CLP.setOption("symmetric", "unsymmetric", &symmetric,
167 "Symmetric discretization");
168
169 int num_spatial_procs = -1;
170 CLP.setOption("num_spatial_procs", &num_spatial_procs, "Number of spatial processors (set -1 for all available procs)");
171
172 SG_RF randField = UNIFORM;
173 CLP.setOption("rand_field", &randField,
175 "Random field type");
176
177 double mu = 0.2;
178 CLP.setOption("mean", &mu, "Mean");
179
180 double s = 0.1;
181 CLP.setOption("std_dev", &s, "Standard deviation");
182
183 int num_KL = 2;
184 CLP.setOption("num_kl", &num_KL, "Number of KL terms");
185
186 int order = 3;
187 CLP.setOption("order", &order, "Polynomial order");
188
189 bool normalize_basis = true;
190 CLP.setOption("normalize", "unnormalize", &normalize_basis,
191 "Normalize PC basis");
192
193 Krylov_Method solver_method = GMRES;
194 CLP.setOption("solver_method", &solver_method,
196 "Krylov solver method");
197
198 SG_Prec prec_method = STOCHASTIC;
199 CLP.setOption("prec_method", &prec_method,
201 "Preconditioner method");
202
203 SG_Div division_method = DIRECT;
204 CLP.setOption("division_method", &division_method,
206 "Stochastic division method");
207
208 SG_DivPrec divprec_method = NO;
209 CLP.setOption("divprec_method", &divprec_method,
211 "Preconditioner for division method");
212 Schur_option schur_option = diag;
213 CLP.setOption("schur_option", &schur_option,
215 "Schur option");
216 Prec_option prec_option = whole;
217 CLP.setOption("prec_option", &prec_option,
219 "Prec option");
220
221
222 double solver_tol = 1e-12;
223 CLP.setOption("solver_tol", &solver_tol, "Outer solver tolerance");
224
225 double div_tol = 1e-6;
226 CLP.setOption("div_tol", &div_tol, "Tolerance in Iterative Solver");
227
228 int prec_level = 1;
229 CLP.setOption("prec_level", &prec_level, "Level in Schur Complement Prec 0->Solve A0u0=g0 with division; 1->Form 1x1 Schur Complement");
230
231 int max_it_div = 50;
232 CLP.setOption("max_it_div", &max_it_div, "Maximum # of Iterations in Iterative Solver for Division");
233
234 bool equilibrate = true; //JJH 8/26/12 changing to true to match ETP example
235 CLP.setOption("equilibrate", "noequilibrate", &equilibrate,
236 "Equilibrate the linear system");
237
238
239 CLP.parse( argc, argv );
240
241 if (MyPID == 0) {
242 std::cout << "Summary of command line options:" << std::endl
243 << "\tnum_mesh = " << n << std::endl
244 << "\tsymmetric = " << symmetric << std::endl
245 << "\tnum_spatial_procs = " << num_spatial_procs << std::endl
246 << "\trand_field = " << sg_rf_names[randField]
247 << std::endl
248 << "\tmean = " << mu << std::endl
249 << "\tstd_dev = " << s << std::endl
250 << "\tnum_kl = " << num_KL << std::endl
251 << "\torder = " << order << std::endl
252 << "\tnormalize_basis = " << normalize_basis << std::endl
253 << "\tsolver_method = " << krylov_method_names[solver_method] << std::endl
254 << "\tprec_method = " << sg_prec_names[prec_method] << std::endl
255 << "\tdivision_method = " << sg_div_names[division_method] << std::endl
256 << "\tdiv_tol = " << div_tol << std::endl
257 << "\tdiv_prec = " << sg_divprec_names[divprec_method] << std::endl
258 << "\tprec_level = " << prec_level << std::endl
259 << "\tmax_it_div = " << max_it_div << std::endl;
260 }
261 bool nonlinear_expansion = false;
262 if (randField == UNIFORM)
263 nonlinear_expansion = false;
264 else if (randField == LOGNORMAL)
265 nonlinear_expansion = true;
266
267 {
268 TEUCHOS_FUNC_TIME_MONITOR("Total PCE Calculation Time");
269
270 // Create Stochastic Galerkin basis and expansion
271 Teuchos::Array< RCP<const Stokhos::OneDOrthogPolyBasis<LocalOrdinal,BasisScalar> > > bases(num_KL);
272 for (LocalOrdinal i=0; i<num_KL; i++)
273 if (randField == UNIFORM)
274 bases[i] = rcp(new Stokhos::LegendreBasis<LocalOrdinal,BasisScalar>(order, normalize_basis));
275 else if (randField == LOGNORMAL)
276 bases[i] = rcp(new Stokhos::HermiteBasis<int,double>(order, normalize_basis));
277 RCP<const Stokhos::CompletePolynomialBasis<LocalOrdinal,BasisScalar> > basis =
279 LocalOrdinal sz = basis->size();
280 RCP<Stokhos::Sparse3Tensor<LocalOrdinal,BasisScalar> > Cijk =
281 basis->computeTripleProductTensor();
282 RCP<const Stokhos::Quadrature<int,double> > quad =
284 RCP<ParameterList> expn_params = Teuchos::rcp(new ParameterList);
285 if (division_method == MEAN_DIV) {
286 expn_params->set("Division Strategy", "Mean-Based");
287 expn_params->set("Use Quadrature for Division", false);
288 }
289 else if (division_method == DIRECT) {
290 expn_params->set("Division Strategy", "Dense Direct");
291 expn_params->set("Use Quadrature for Division", false);
292 }
293 else if (division_method == SPD_DIRECT) {
294 expn_params->set("Division Strategy", "SPD Dense Direct");
295 expn_params->set("Use Quadrature for Division", false);
296 }
297 else if (division_method == CGD) {
298 expn_params->set("Division Strategy", "CG");
299 expn_params->set("Use Quadrature for Division", false);
300 }
301
302 else if (division_method == QUAD) {
303 expn_params->set("Use Quadrature for Division", true);
304 }
305
306 if (divprec_method == NO)
307 expn_params->set("Prec Strategy", "None");
308 else if (divprec_method == DIAG)
309 expn_params->set("Prec Strategy", "Diag");
310 else if (divprec_method == JACOBI)
311 expn_params->set("Prec Strategy", "Jacobi");
312 else if (divprec_method == GS)
313 expn_params->set("Prec Strategy", "GS");
314 else if (divprec_method == SCHUR)
315 expn_params->set("Prec Strategy", "Schur");
316
317 if (schur_option == diag)
318 expn_params->set("Schur option", "diag");
319 else
320 expn_params->set("Schur option", "full");
321 if (prec_option == linear)
322 expn_params->set("Prec option", "linear");
323
324
325 if (equilibrate)
326 expn_params->set("Equilibrate", 1);
327 else
328 expn_params->set("Equilibrate", 0);
329 expn_params->set("Division Tolerance", div_tol);
330 expn_params->set("prec_iter", prec_level);
331 expn_params->set("max_it_div", max_it_div);
332
333 RCP<Stokhos::OrthogPolyExpansion<LocalOrdinal,BasisScalar> > expansion =
335 basis, Cijk, quad, expn_params));
336
337 if (MyPID == 0)
338 std::cout << "Stochastic Galerkin expansion size = " << sz << std::endl;
339
340 // Create stochastic parallel distribution
341 ParameterList parallelParams;
342 parallelParams.set("Number of Spatial Processors", num_spatial_procs);
343 // parallelParams.set("Rebalance Stochastic Graph", true);
344 // Teuchos::ParameterList& isorropia_params =
345 // parallelParams.sublist("Isorropia");
346 // isorropia_params.set("Balance objective", "nonzeros");
347 RCP<Stokhos::ParallelData> sg_parallel_data =
348 rcp(new Stokhos::ParallelData(basis, Cijk, globalComm, parallelParams));
349 RCP<const EpetraExt::MultiComm> sg_comm =
350 sg_parallel_data->getMultiComm();
351 RCP<const Epetra_Comm> app_comm =
352 sg_parallel_data->getSpatialComm();
353
354 // Create Teuchos::Comm from Epetra_Comm
355 RCP< Teuchos::Comm<int> > teuchos_app_comm;
356#ifdef HAVE_MPI
357 RCP<const Epetra_MpiComm> app_mpi_comm =
358 Teuchos::rcp_dynamic_cast<const Epetra_MpiComm>(app_comm);
359 RCP<const Teuchos::OpaqueWrapper<MPI_Comm> > raw_mpi_comm =
360 Teuchos::opaqueWrapper(app_mpi_comm->Comm());
361 teuchos_app_comm = rcp(new Teuchos::MpiComm<int>(raw_mpi_comm));
362#else
363 teuchos_app_comm = rcp(new Teuchos::SerialComm<int>());
364#endif
365
366 // Create application
368 RCP<problem_type> model =
369 rcp(new problem_type(teuchos_app_comm, n, num_KL, s, mu,
370 nonlinear_expansion, symmetric));
371
372 // Create vectors and operators
373 typedef problem_type::Tpetra_Vector Tpetra_Vector;
374 typedef problem_type::Tpetra_CrsMatrix Tpetra_CrsMatrix;
375 typedef Tpetra::MatrixMarket::Writer<Tpetra_CrsMatrix> Writer;
376
377 RCP<Tpetra_Vector> p = Tpetra::createVector<Scalar>(model->get_p_map(0));
378 RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(model->get_x_map());
379 x->putScalar(0.0);
380 RCP<Tpetra_Vector> f = Tpetra::createVector<Scalar>(model->get_f_map());
381 RCP<Tpetra_Vector> dx = Tpetra::createVector<Scalar>(model->get_x_map());
382 RCP<Tpetra_CrsMatrix> J = model->create_W();
383 RCP<Tpetra_CrsMatrix> J0;
384 if (prec_method == MEAN)
385 J0 = model->create_W();
386
387 // Set PCE expansion of p
388 p->putScalar(0.0);
389 ArrayRCP<Scalar> p_view = p->get1dViewNonConst();
390 for (ArrayRCP<Scalar>::size_type i=0; i<p_view.size(); i++) {
391 p_view[i].reset(expansion);
392 p_view[i].copyForWrite();
393 }
394 Array<double> point(num_KL, 1.0);
395 Array<double> basis_vals(sz);
396 basis->evaluateBases(point, basis_vals);
397 if (order > 0) {
398 for (int i=0; i<num_KL; i++) {
399 p_view[i].term(i,1) = 1.0 / basis_vals[i+1];
400 }
401 }
402
403 // Create preconditioner
404 typedef Ifpack2::Preconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tprec;
405 Teuchos::RCP<Tprec> M;
406 if (prec_method != NONE) {
407 ParameterList precParams;
408 std::string prec_name = "RILUK";
409 precParams.set("fact: iluk level-of-fill", 4);
410 precParams.set("fact: iluk level-of-overlap", 0);
411
412 // std::string prec_name = "ILUT";
413 // precParams.set("fact: ilut level-of-fill", 10.0);
414 // precParams.set("fact: drop tolerance", 1.0e-16);
415 Ifpack2::Factory factory;
416 if (prec_method == MEAN) {
417 M = factory.create<Tpetra_CrsMatrix>(prec_name, J0);
418 } else if (prec_method == STOCHASTIC) {
419 M = factory.create<Tpetra_CrsMatrix>(prec_name, J);
420 }
421 M->setParameters(precParams);
422 }
423
424 // Evaluate model
425 model->computeResidual(*x, *p, *f);
426 model->computeJacobian(*x, *p, *J);
427
428 // Compute mean for mean-based preconditioner
429 if (prec_method == MEAN) {
430 size_t nrows = J->getLocalNumRows();
431 ArrayView<const LocalOrdinal> indices;
432 ArrayView<const Scalar> values;
433 J0->resumeFill();
434 for (size_t i=0; i<nrows; i++) {
435 J->getLocalRowView(i, indices, values);
436 Array<Scalar> values0(values.size());
437 for (LocalOrdinal j=0; j<values.size(); j++)
438 values0[j] = values[j].coeff(0);
439 J0->replaceLocalValues(i, indices, values0);
440 }
441 J0->fillComplete();
442 }
443
444 // compute preconditioner
445 if (prec_method != NONE) {
446 M->initialize();
447 M->compute();
448 }
449
450 // Setup Belos solver
451 RCP<ParameterList> belosParams = rcp(new ParameterList);
452 belosParams->set("Flexible Gmres", false);
453 belosParams->set("Num Blocks", 500);//20
454 belosParams->set("Convergence Tolerance", solver_tol);
455 belosParams->set("Maximum Iterations", 1000);
456 belosParams->set("Verbosity", 33);
457 belosParams->set("Output Style", 1);
458 belosParams->set("Output Frequency", 1);
459
460 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
461 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
462 typedef Belos::OperatorTraits<Scalar,MV,OP> BOPT;
463 typedef Belos::MultiVecTraits<double,MV> BMVT;
464 typedef Belos::LinearProblem<double,MV,OP> BLinProb;
465 RCP< BLinProb > problem = rcp(new BLinProb(J, dx, f));
466 if (prec_method != NONE)
467 problem->setRightPrec(M);
468 problem->setProblem();
469 RCP<Belos::SolverManager<double,MV,OP> > solver;
470 if (solver_method == CG)
471 solver = rcp(new Belos::PseudoBlockCGSolMgr<double,MV,OP>(problem, belosParams));
472 else if (solver_method == GMRES)
473 solver = rcp(new Belos::BlockGmresSolMgr<double,MV,OP>(problem, belosParams));
474
475 // Print initial residual norm
476 std::vector<double> norm_f(1);
477 BMVT::MvNorm(*f, norm_f);
478 if (MyPID == 0)
479 std::cout << "\nInitial residual norm = " << norm_f[0] << std::endl;
480
481 // Solve linear system
482 Belos::ReturnType ret = solver->solve();
483
484 if (MyPID == 0) {
485 if (ret == Belos::Converged)
486 std::cout << "Solver converged!" << std::endl;
487 else
488 std::cout << "Solver failed to converge!" << std::endl;
489 }
490
491 // Update x
492 x->update(-1.0, *dx, 1.0);
493 Writer::writeDenseFile("stochastic_solution.mm", x);
494
495 // Compute new residual & response function
496 RCP<Tpetra_Vector> g = Tpetra::createVector<Scalar>(model->get_g_map(0));
497 f->putScalar(0.0);
498 model->computeResidual(*x, *p, *f);
499 model->computeResponse(*x, *p, *g);
500
501 // Print final residual norm
502 BMVT::MvNorm(*f, norm_f);
503 if (MyPID == 0)
504 std::cout << "\nFinal residual norm = " << norm_f[0] << std::endl;
505
506 // Expected results for num_mesh=32
507 double g_mean_exp = 0.172988; // expected response mean
508 double g_std_dev_exp = 0.0380007; // expected response std. dev.
509 double g_tol = 1e-6; // tolerance on determining success
510 if (n == 8) {
511 g_mean_exp = 1.327563e-01;
512 g_std_dev_exp = 2.949064e-02;
513 }
514
515 double g_mean = g->get1dView()[0].mean();
516 double g_std_dev = g->get1dView()[0].standard_deviation();
517 std::cout << std::endl;
518 std::cout << "Response Mean = " << g_mean << std::endl;
519 std::cout << "Response Std. Dev. = " << g_std_dev << std::endl;
520 bool passed = false;
521 if (norm_f[0] < 1.0e-10 &&
522 std::abs(g_mean-g_mean_exp) < g_tol &&
523 std::abs(g_std_dev - g_std_dev_exp) < g_tol)
524 passed = true;
525 if (MyPID == 0) {
526 if (passed)
527 std::cout << "Example Passed!" << std::endl;
528 else{
529 std::cout << "Example Failed!" << std::endl;
530 std::cout << "Expected Response Mean = "<< g_mean_exp << std::endl;
531 std::cout << "Expected Response Std. Dev. = "<< g_std_dev_exp << std::endl;
532 }
533 }
534
535 }
536
537 Teuchos::TimeMonitor::summarize(std::cout);
538 Teuchos::TimeMonitor::zeroOutTimers();
539
540 }
541
542 catch (std::exception& e) {
543 std::cout << e.what() << std::endl;
544 }
545 catch (string& s) {
546 std::cout << s << std::endl;
547 }
548 catch (char *s) {
549 std::cout << s << std::endl;
550 }
551 catch (...) {
552 std::cout << "Caught unknown exception!" <<std:: endl;
553 }
554
555#ifdef HAVE_MPI
556 MPI_Finalize() ;
557#endif
558
559}
expr expr expr dx(i, j)
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Hermite polynomial basis.
Legendre polynomial basis.
Orthogonal polynomial expansions based on numerical quadrature.
Defines quadrature for a tensor product basis by tensor products of 1-D quadrature rules.
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)
const int num_krylov_method
const int num_sg_div
const char * prec_option_names[]
int main(int argc, char *argv[])
const char * sg_divprec_names[]
const int num_prec_option
const SG_Div sg_div_values[]
const int num_sg_prec
const int num_schur_option
const Prec_option Prec_option_values[]
const SG_DivPrec sg_divprec_values[]
const char * krylov_method_names[]
const char * sg_div_names[]
const Krylov_Method krylov_method_values[]
const char * sg_prec_names[]
const Schur_option Schur_option_values[]
const SG_Prec sg_prec_values[]
const char * sg_rf_names[]
const int num_sg_divprec
const int num_sg_rf
const SG_RF sg_rf_values[]
const char * schur_option_names[]