Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
linear2d_diffusion_collocation_strat.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// ModelEvaluator implementing our problem
43#include "twoD_diffusion_ME.hpp"
44
45// Epetra communicator
46#ifdef HAVE_MPI
47#include "Epetra_MpiComm.h"
48#else
49#include "Epetra_SerialComm.h"
50#endif
51
52// Stratimikos stuff
53#include "Stratimikos_DefaultLinearSolverBuilder.hpp"
54#include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
55#include "Thyra_PreconditionerBase.hpp"
56#include "Thyra_PreconditionerFactoryBase.hpp"
57#include "Thyra_EpetraThyraWrappers.hpp"
58#include "Thyra_EpetraLinearOp.hpp"
59
60// Stokhos Stochastic Galerkin
61#include "Stokhos.hpp"
62
63// Timing utilities
64#include "Teuchos_CommandLineProcessor.hpp"
65#include "Teuchos_TimeMonitor.hpp"
66
67// I/O utilities
68#include "EpetraExt_VectorOut.h"
69
70// Linear solvers
72const int num_krylov_solver = 2;
74const char *krylov_solver_names[] = { "AztecOO", "Belos" };
75
76// Krylov methods
78const int num_krylov_method = 4;
80const char *krylov_method_names[] = { "GMRES", "CG", "FGMRES", "RGMRES" };
81
82// Collocation preconditioning strategies
84const int num_prec_strategy = 2;
86const char *prec_strategy_names[] = { "Mean", "Rebuild" };
87
88// Random field types
90const int num_sg_rf = 4;
92const char *sg_rf_names[] = { "Uniform", "CC-Uniform", "Rys", "Log-Normal" };
93
94// Sparse grid growth rules
96const int num_sg_growth = 3;
99const char *sg_growth_names[] = {
100 "Slow Restricted", "Moderate Restricted", "Unrestricted" };
101
102int main(int argc, char *argv[]) {
103
104// Initialize MPI
105#ifdef HAVE_MPI
106 MPI_Init(&argc,&argv);
107#endif
108
109 // Create a communicator for Epetra objects
110 Teuchos::RCP<Epetra_Comm> Comm;
111#ifdef HAVE_MPI
112 Comm = Teuchos::rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
113#else
114 Comm = Teuchos::rcp(new Epetra_SerialComm);
115#endif
116
117 int MyPID = Comm->MyPID();
118
119 try {
120
121 // Setup command line options
122 Teuchos::CommandLineProcessor CLP;
123 CLP.setDocString(
124 "This example runs a stochastic collocation method.\n");
125
126 int n = 32;
127 CLP.setOption("num_mesh", &n, "Number of mesh points in each direction");
128
129 SG_RF randField = UNIFORM;
130 CLP.setOption("rand_field", &randField,
132 "Random field type");
133
134 double mean = 0.2;
135 CLP.setOption("mean", &mean, "Mean");
136
137 double sigma = 0.1;
138 CLP.setOption("std_dev", &sigma, "Standard deviation");
139
140 double weightCut = 1.0;
141 CLP.setOption("weight_cut", &weightCut, "Weight cut");
142
143 int num_KL = 2;
144 CLP.setOption("num_kl", &num_KL, "Number of KL terms");
145
146 int p = 3;
147 CLP.setOption("order", &p, "Polynomial order");
148
149 bool normalize_basis = true;
150 CLP.setOption("normalize", "unnormalize", &normalize_basis,
151 "Normalize PC basis");
152
153 Krylov_Solver krylov_solver = AZTECOO;
154 CLP.setOption("krylov_solver", &krylov_solver,
156 "Linear solver");
157
158 Krylov_Method krylov_method = CG;
159 CLP.setOption("krylov_method", &krylov_method,
161 "Krylov method");
162
163 PrecStrategy precStrategy = MEAN;
164 CLP.setOption("prec_strategy", &precStrategy,
166 "Preconditioner strategy");
167
168 double tol = 1e-12;
169 CLP.setOption("tol", &tol, "Solver tolerance");
170
171 SG_GROWTH sg_growth = MODERATE_RESTRICTED;
172 CLP.setOption("sg_growth", &sg_growth,
174 "Sparse grid growth rule");
175
176 CLP.parse( argc, argv );
177
178 if (MyPID == 0) {
179 std::cout << "Summary of command line options:" << std::endl
180 << "\tnum_mesh = " << n << std::endl
181 << "\trand_field = " << sg_rf_names[randField] << std::endl
182 << "\tmean = " << mean << std::endl
183 << "\tstd_dev = " << sigma << std::endl
184 << "\tweight_cut = " << weightCut << std::endl
185 << "\tnum_kl = " << num_KL << std::endl
186 << "\torder = " << p << std::endl
187 << "\tnormalize_basis = " << normalize_basis << std::endl
188 << "\tkrylov_solver = " << krylov_solver_names[krylov_solver]
189 << std::endl
190 << "\tkrylov_method = " << krylov_method_names[krylov_method]
191 << std::endl
192 << "\tprec_strategy = " << prec_strategy_names[precStrategy]
193 << std::endl
194 << "\ttol = " << tol << std::endl
195 << "\tsg_growth = " << sg_growth_names[sg_growth]
196 << std::endl;
197 }
198
199 bool nonlinear_expansion = false;
200 if (randField == LOGNORMAL)
201 nonlinear_expansion = true;
202
203 {
204 TEUCHOS_FUNC_TIME_MONITOR("Total Collocation Calculation Time");
205
206 // Create Stochastic Galerkin basis
207 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(num_KL);
208 for (int i=0; i<num_KL; i++) {
209 Teuchos::RCP<Stokhos::OneDOrthogPolyBasis<int,double> > b;
210 if (randField == UNIFORM) {
211 b = Teuchos::rcp(new Stokhos::LegendreBasis<int,double>(
212 p, normalize_basis));
213 }
214 else if (randField == CC_UNIFORM) {
215 b =
217 p, normalize_basis, true));
218 }
219 else if (randField == RYS) {
220 b = Teuchos::rcp(new Stokhos::RysBasis<int,double>(
221 p, weightCut, normalize_basis));
222 }
223 else if (randField == LOGNORMAL) {
224 b = Teuchos::rcp(new Stokhos::HermiteBasis<int,double>(
225 p, normalize_basis));
226 }
227 bases[i] = b;
228 }
229 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
230 Teuchos::rcp(new Stokhos::CompletePolynomialBasis<int,double>(bases));
231
232 // Create sparse grid
233 int sparse_grid_growth = Pecos::MODERATE_RESTRICTED_GROWTH;
234 if (sg_growth == SLOW_RESTRICTED)
235 sparse_grid_growth = Pecos::SLOW_RESTRICTED_GROWTH;
236 else if (sg_growth == MODERATE_RESTRICTED)
237 sparse_grid_growth = Pecos::MODERATE_RESTRICTED_GROWTH;
238 else if (sg_growth == UNRESTRICTED)
239 sparse_grid_growth = Pecos::UNRESTRICTED_GROWTH;
240 Stokhos::SparseGridQuadrature<int,double> quad(basis,p,1e-12,sg_growth);
241 const Teuchos::Array< Teuchos::Array<double> >& quad_points =
242 quad.getQuadPoints();
243 const Teuchos::Array<double>& quad_weights =
244 quad.getQuadWeights();
245 int nqp = quad_weights.size();
246
247 // Create application
248 twoD_diffusion_ME model(Comm, n, num_KL, sigma, mean,
249 basis, nonlinear_expansion);
250
251 // Model data
252 Teuchos::RCP<Epetra_Vector> p =
253 Teuchos::rcp(new Epetra_Vector(*(model.get_p_map(0))));
254 Teuchos::RCP<Epetra_Vector> x =
255 Teuchos::rcp(new Epetra_Vector(*(model.get_x_map())));
256 Teuchos::RCP<Epetra_Vector> f =
257 Teuchos::rcp(new Epetra_Vector(*(model.get_f_map())));
258 Teuchos::RCP<Epetra_Vector> g =
259 Teuchos::rcp(new Epetra_Vector(*(model.get_g_map(0))));
260 Teuchos::RCP<Epetra_CrsMatrix> A =
261 Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(model.create_W());
262 EpetraExt::ModelEvaluator::InArgs inArgs = model.createInArgs();
263 EpetraExt::ModelEvaluator::OutArgs outArgs = model.createOutArgs();
264 EpetraExt::ModelEvaluator::OutArgs outArgs2 = model.createOutArgs();
265
266 // Data to compute/store mean & variance
267 Epetra_Vector x2(*(model.get_x_map()));
268 Epetra_Vector x_mean(*(model.get_x_map()));
269 Epetra_Vector x_var(*(model.get_x_map()));
270 Epetra_Vector g2(*(model.get_g_map(0)));
271 Epetra_Vector g_mean(*(model.get_g_map(0)));
272 Epetra_Vector g_var(*(model.get_g_map(0)));
273
274 // Create ParameterList controlling Stratimikos
275 Teuchos::ParameterList stratParams;
276 if (krylov_solver == AZTECOO) {
277 stratParams.set("Linear Solver Type", "AztecOO");
278 Teuchos::ParameterList& aztecParams =
279 stratParams.sublist("Linear Solver Types").sublist("AztecOO").sublist("Forward Solve");
280 aztecParams.set("Max Iterations", 1000);
281 aztecParams.set("Tolerance", tol);
282 Teuchos::ParameterList& aztecSettings =
283 aztecParams.sublist("AztecOO Settings");
284 if (krylov_method == GMRES)
285 aztecSettings.set("Aztec Solver", "GMRES");
286 else if (krylov_method == CG)
287 aztecSettings.set("Aztec Solver", "CG");
288 aztecSettings.set("Aztec Preconditioner", "none");
289 aztecSettings.set("Size of Krylov Subspace", 100);
290 aztecSettings.set("Convergence Test", "r0");
291 aztecSettings.set("Output Frequency", 10);
292 Teuchos::ParameterList& verbParams =
293 stratParams.sublist("Linear Solver Types").sublist("AztecOO").sublist("VerboseObject");
294 verbParams.set("Verbosity Level", "none");
295 }
296 else if (krylov_solver == BELOS) {
297 stratParams.set("Linear Solver Type", "Belos");
298 Teuchos::ParameterList& belosParams =
299 stratParams.sublist("Linear Solver Types").sublist("Belos");
300 Teuchos::ParameterList* belosSolverParams = NULL;
301 if (krylov_method == GMRES || krylov_method == FGMRES) {
302 belosParams.set("Solver Type","Block GMRES");
303 belosSolverParams =
304 &(belosParams.sublist("Solver Types").sublist("Block GMRES"));
305 if (krylov_method == FGMRES)
306 belosSolverParams->set("Flexible Gmres", true);
307 }
308 else if (krylov_method == RGMRES) {
309 belosParams.set("Solver Type","GCRODR");
310 belosSolverParams =
311 &(belosParams.sublist("Solver Types").sublist("GCRODR"));
312 belosSolverParams->set("Num Recycled Blocks", 10);
313 }
314 else if (krylov_method == CG) {
315 belosParams.set("Solver Type","Block CG");
316 belosSolverParams =
317 &(belosParams.sublist("Solver Types").sublist("Block CG"));
318 }
319
320 belosSolverParams->set("Convergence Tolerance", tol);
321 belosSolverParams->set("Maximum Iterations", 1000);
322 belosSolverParams->set("Num Blocks", 100);
323 belosSolverParams->set("Output Frequency",10);
324 Teuchos::ParameterList& verbParams = belosParams.sublist("VerboseObject");
325 verbParams.set("Verbosity Level", "none");
326 }
327 stratParams.set("Preconditioner Type", "ML");
328 Teuchos::ParameterList& precParams =
329 stratParams.sublist("Preconditioner Types").sublist("ML").sublist("ML Settings");
330 precParams.set("default values", "SA");
331 precParams.set("ML output", 0);
332 precParams.set("max levels",5);
333 precParams.set("increasing or decreasing","increasing");
334 precParams.set("aggregation: type", "Uncoupled");
335 precParams.set("smoother: type","ML symmetric Gauss-Seidel");
336 precParams.set("smoother: sweeps",2);
337 precParams.set("smoother: pre or post", "both");
338 precParams.set("coarse: max size", 200);
339#ifdef HAVE_ML_AMESOS
340 precParams.set("coarse: type","Amesos-KLU");
341#else
342 precParams.set("coarse: type","Jacobi");
343#endif
344
345 // Create Stratimikos linear solver builder
346 Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
347 linearSolverBuilder.setParameterList(Teuchos::rcp(&stratParams, false));
348
349 // Create a linear solver factory given information read from the
350 // parameter list.
351 Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<double> > lowsFactory =
352 linearSolverBuilder.createLinearSolveStrategy("");
353 Teuchos::RCP<Teuchos::FancyOStream> out =
354 Teuchos::VerboseObjectBase::getDefaultOStream();
355 lowsFactory->setOStream(out);
356 lowsFactory->setVerbLevel(Teuchos::VERB_LOW);
357
358 // Initialize the LinearOpWithSolve
359 Teuchos::RCP<Thyra::LinearOpWithSolveBase<double> > lows =
360 lowsFactory->createOp();
361 Teuchos::RCP<Thyra::PreconditionerFactoryBase<double> > precFactory =
362 lowsFactory->getPreconditionerFactory();
363 Teuchos::RCP<Thyra::PreconditionerBase<double> > M_thyra =
364 precFactory->createPrec();
365
366 // Wrap Thyra objects around Epetra objects
367 Teuchos::RCP<const Thyra::LinearOpBase<double> > A_thyra =
368 Thyra::epetraLinearOp(A);
369 Teuchos::RCP<const Thyra::LinearOpSourceBase<double> > losb =
370 rcp(new Thyra::DefaultLinearOpSource<double>(A_thyra));
371
372 // Create solver convergence criteria
373 Teuchos::RCP<Thyra::SolveCriteria<double> > solveCriteria;
374 if (!(krylov_solver == BELOS && krylov_method == CG)) {
375 // For some reason, Belos' Block-CG doesn't support this
376 Thyra::SolveMeasureType solveMeasure(
377 Thyra::SOLVE_MEASURE_NORM_RESIDUAL,
378 Thyra::SOLVE_MEASURE_NORM_INIT_RESIDUAL);
379 solveCriteria =
380 Teuchos::rcp(new Thyra::SolveCriteria<double>(solveMeasure, tol));
381 }
382
383 // Computation of prec happens here:
384 if (precStrategy == MEAN) {
385 TEUCHOS_FUNC_TIME_MONITOR("Deterministic Preconditioner Calculation");
386 *A = *(model.get_mean());
387 precFactory->initializePrec(losb, M_thyra.get());
388 Thyra::initializePreconditionedOp<double>(
389 *lowsFactory, A_thyra, M_thyra, lows.ptr());
390 }
391
392 x_mean.PutScalar(0.0);
393 x_var.PutScalar(0.0);
394 // Loop over colloction points
395 for (int qp=0; qp<nqp; qp++) {
396 TEUCHOS_FUNC_TIME_MONITOR("Collocation Loop");
397
398 // Set parameters
399 for (int i=0; i<num_KL; i++)
400 (*p)[i] = quad_points[qp][i];
401
402 // Set in/out args
403 inArgs.set_p(0, p);
404 inArgs.set_x(x);
405 outArgs.set_f(f);
406 outArgs.set_W(A);
407
408 // Evaluate model at collocation point
409 x->PutScalar(0.0);
410 model.evalModel(inArgs, outArgs);
411 f->Scale(-1.0);
412
413 Teuchos::RCP<Thyra::VectorBase<double> > x_thyra =
414 Thyra::create_Vector(x, A_thyra->domain());
415 Teuchos::RCP<const Thyra::VectorBase<double> > f_thyra =
416 Thyra::create_Vector(f, A_thyra->range());
417
418 // Compute preconditioner
419 if (precStrategy != MEAN) {
420 TEUCHOS_FUNC_TIME_MONITOR("Deterministic Preconditioner Calculation");
421 precFactory->initializePrec(losb, M_thyra.get());
422 Thyra::initializePreconditionedOp<double>(
423 *lowsFactory, A_thyra, M_thyra, lows.ptr());
424 }
425
426 // Solve linear system
427 {
428 TEUCHOS_FUNC_TIME_MONITOR("Deterministic Solve");
429 Thyra::SolveStatus<double> solveStatus =
430 lows->solve(Thyra::NOTRANS, *f_thyra, x_thyra.ptr(),
431 solveCriteria.ptr());
432 if (MyPID == 0) {
433 std::cout << "Collocation point " << qp+1 <<'/' << nqp << ": "
434 << solveStatus.message << std::endl;
435 }
436 }
437
438 x_thyra = Teuchos::null;
439 f_thyra = Teuchos::null;
440
441 // Compute responses
442 outArgs2.set_g(0, g);
443 model.evalModel(inArgs, outArgs2);
444
445 // Sum contributions to mean and variance
446 x2.Multiply(1.0, *x, *x, 0.0);
447 g2.Multiply(1.0, *g, *g, 0.0);
448 x_mean.Update(quad_weights[qp], *x, 1.0);
449 x_var.Update(quad_weights[qp], x2, 1.0);
450 g_mean.Update(quad_weights[qp], *g, 1.0);
451 g_var.Update(quad_weights[qp], g2, 1.0);
452
453 }
454 x2.Multiply(1.0, x_mean, x_mean, 0.0);
455 g2.Multiply(1.0, g_mean, g_mean, 0.0);
456 x_var.Update(-1.0, x2, 1.0);
457 g_var.Update(-1.0, g2, 1.0);
458
459 // Compute standard deviations
460 for (int i=0; i<x_var.MyLength(); i++)
461 x_var[i] = std::sqrt(x_var[i]);
462 for (int i=0; i<g_var.MyLength(); i++)
463 g_var[i] = std::sqrt(g_var[i]);
464
465 std::cout.precision(16);
466 std::cout << "\nResponse Mean = " << std::endl << g_mean << std::endl;
467 std::cout << "Response Std. Dev. = " << std::endl << g_var << std::endl;
468
469 // Save mean and variance to file
470 EpetraExt::VectorToMatrixMarketFile("mean_col.mm", x_mean);
471 EpetraExt::VectorToMatrixMarketFile("std_dev_col.mm", x_var);
472
473 }
474
475 Teuchos::TimeMonitor::summarize(std::cout);
476 Teuchos::TimeMonitor::zeroOutTimers();
477
478 }
479
480 catch (std::exception& e) {
481 std::cout << e.what() << std::endl;
482 }
483 catch (std::string& s) {
484 std::cout << s << std::endl;
485 }
486 catch (char *s) {
487 std::cout << s << std::endl;
488 }
489 catch (...) {
490 std::cout << "Caught unknown exception!" <<std:: endl;
491 }
492
493#ifdef HAVE_MPI
494 MPI_Finalize() ;
495#endif
496
497}
Legendre polynomial basis using Clenshaw-Curtis quadrature points.
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Hermite polynomial basis.
Legendre polynomial basis.
Rys polynomial basis.
ModelEvaluator for a linear 2-D diffusion problem.
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[])
const PrecStrategy prec_strategy_values[]
const char * krylov_method_names[]
const Krylov_Method krylov_method_values[]
const char * krylov_solver_names[]
const Krylov_Solver krylov_solver_values[]
const char * sg_rf_names[]
const char * sg_growth_names[]
const SG_GROWTH sg_growth_values[]
const SG_RF sg_rf_values[]
const char * prec_strategy_names[]