Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
linear2d_diffusion_pce_nox_sg_solvers.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// NOX
53#include "NOX.H"
54#include "NOX_Epetra.H"
55#include "NOX_Epetra_LinearSystem_Stratimikos.H"
58
59// Stokhos Stochastic Galerkin
60#include "Stokhos_Epetra.hpp"
61
62// Utilities
63#include "Teuchos_CommandLineProcessor.hpp"
64#include "Teuchos_TimeMonitor.hpp"
65#include "Teuchos_Assert.hpp"
66
67// I/O utilities
68#include "EpetraExt_VectorOut.h"
69
70//The probability distribution of the random variables.
71double uniform_weight(const double& x){
72 return 1;
73}
74
75// Linear solvers
77const int num_krylov_solver = 2;
79const char *krylov_solver_names[] = { "AztecOO", "Belos" };
80
81// SG solver approaches
83const int num_sg_solver = 3;
85const char *sg_solver_names[] = { "Krylov", "Gauss-Seidel", "Jacobi" };
86
87// Krylov methods
89const int num_krylov_method = 4;
91const char *krylov_method_names[] = { "GMRES", "CG", "FGMRES", "RGMRES" };
92
93// Krylov operator approaches
95const int num_sg_op = 4;
97const char *sg_op_names[] = { "Matrix-Free",
98 "KL-Matrix-Free",
99 "KL-Reduced-Matrix-Free",
100 "Fully-Assembled" };
101
102// Krylov preconditioning approaches
103enum SG_Prec { MEAN, GS, AGS, AJ, ASC, KP, NONE };
104const int num_sg_prec = 7;
106const char *sg_prec_names[] = { "Mean-Based",
107 "Gauss-Seidel",
108 "Approx-Gauss-Seidel",
109 "Approx-Jacobi",
110 "Approx-Schur-Complement",
111 "Kronecker-Product",
112 "None" };
113
114// Random field types
116const int num_sg_rf = 3;
118const char *sg_rf_names[] = { "Uniform", "Rys", "Log-Normal" };
119
120int main(int argc, char *argv[]) {
121
122// Initialize MPI
123#ifdef HAVE_MPI
124 MPI_Init(&argc,&argv);
125#endif
126
127 // Create a communicator for Epetra objects
128 Teuchos::RCP<const Epetra_Comm> globalComm;
129#ifdef HAVE_MPI
130 globalComm = Teuchos::rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
131#else
132 globalComm = Teuchos::rcp(new Epetra_SerialComm);
133#endif
134 int MyPID = globalComm->MyPID();
135
136 try {
137
138 // Setup command line options
139 Teuchos::CommandLineProcessor CLP;
140 CLP.setDocString(
141 "This example runs a variety of stochastic Galerkin solvers.\n");
142
143 int n = 32;
144 CLP.setOption("num_mesh", &n, "Number of mesh points in each direction");
145
146 bool symmetric = false;
147 CLP.setOption("symmetric", "unsymmetric", &symmetric,
148 "Symmetric discretization");
149
150 int num_spatial_procs = -1;
151 CLP.setOption("num_spatial_procs", &num_spatial_procs, "Number of spatial processors (set -1 for all available procs)");
152
153 bool rebalance_stochastic_graph = false;
154 CLP.setOption("rebalance", "no-rebalance", &rebalance_stochastic_graph,
155 "Rebalance parallel stochastic graph (requires Isorropia)");
156
157 SG_RF randField = UNIFORM;
158 CLP.setOption("rand_field", &randField,
160 "Random field type");
161
162 double mean = 0.2;
163 CLP.setOption("mean", &mean, "Mean");
164
165 double sigma = 0.1;
166 CLP.setOption("std_dev", &sigma, "Standard deviation");
167
168 double weightCut = 1.0;
169 CLP.setOption("weight_cut", &weightCut, "Weight cut");
170
171 int num_KL = 2;
172 CLP.setOption("num_kl", &num_KL, "Number of KL terms");
173
174 int p = 3;
175 CLP.setOption("order", &p, "Polynomial order");
176
177 bool normalize_basis = true;
178 CLP.setOption("normalize", "unnormalize", &normalize_basis,
179 "Normalize PC basis");
180
181 SG_Solver solve_method = SG_KRYLOV;
182 CLP.setOption("sg_solver", &solve_method,
184 "SG solver method");
185
186 Krylov_Method outer_krylov_method = GMRES;
187 CLP.setOption("outer_krylov_method", &outer_krylov_method,
189 "Outer Krylov method (for Krylov-based SG solver)");
190
191 Krylov_Solver outer_krylov_solver = AZTECOO;
192 CLP.setOption("outer_krylov_solver", &outer_krylov_solver,
194 "Outer linear solver");
195
196 double outer_tol = 1e-12;
197 CLP.setOption("outer_tol", &outer_tol, "Outer solver tolerance");
198
199 int outer_its = 1000;
200 CLP.setOption("outer_its", &outer_its, "Maximum outer iterations");
201
202 Krylov_Method inner_krylov_method = GMRES;
203 CLP.setOption("inner_krylov_method", &inner_krylov_method,
205 "Inner Krylov method (for G-S, Jacobi, etc...)");
206
207 Krylov_Solver inner_krylov_solver = AZTECOO;
208 CLP.setOption("inner_krylov_solver", &inner_krylov_solver,
210 "Inner linear solver");
211
212 double inner_tol = 3e-13;
213 CLP.setOption("inner_tol", &inner_tol, "Inner solver tolerance");
214
215 int inner_its = 1000;
216 CLP.setOption("inner_its", &inner_its, "Maximum inner iterations");
217
218 SG_Op opMethod = MATRIX_FREE;
219 CLP.setOption("sg_operator_method", &opMethod,
221 "Operator method");
222
223 SG_Prec precMethod = AGS;
224 CLP.setOption("sg_prec_method", &precMethod,
226 "Preconditioner method");
227
228 double gs_prec_tol = 1e-1;
229 CLP.setOption("gs_prec_tol", &gs_prec_tol, "Gauss-Seidel preconditioner tolerance");
230
231 int gs_prec_its = 1;
232 CLP.setOption("gs_prec_its", &gs_prec_its, "Maximum Gauss-Seidel preconditioner iterations");
233
234 CLP.parse( argc, argv );
235
236 if (MyPID == 0) {
237 std::cout << "Summary of command line options:" << std::endl
238 << "\tnum_mesh = " << n << std::endl
239 << "\tsymmetric = " << symmetric << std::endl
240 << "\tnum_spatial_procs = " << num_spatial_procs << std::endl
241 << "\trebalance = " << rebalance_stochastic_graph
242 << std::endl
243 << "\trand_field = " << sg_rf_names[randField]
244 << std::endl
245 << "\tmean = " << mean << std::endl
246 << "\tstd_dev = " << sigma << std::endl
247 << "\tweight_cut = " << weightCut << std::endl
248 << "\tnum_kl = " << num_KL << std::endl
249 << "\torder = " << p << std::endl
250 << "\tnormalize_basis = " << normalize_basis << std::endl
251 << "\tsg_solver = " << sg_solver_names[solve_method]
252 << std::endl
253 << "\touter_krylov_method = "
254 << krylov_method_names[outer_krylov_method] << std::endl
255 << "\touter_krylov_solver = "
256 << krylov_solver_names[outer_krylov_solver] << std::endl
257 << "\touter_tol = " << outer_tol << std::endl
258 << "\touter_its = " << outer_its << std::endl
259 << "\tinner_krylov_method = "
260 << krylov_method_names[inner_krylov_method] << std::endl
261 << "\tinner_krylov_solver = "
262 << krylov_solver_names[inner_krylov_solver] << std::endl
263 << "\tinner_tol = " << inner_tol << std::endl
264 << "\tinner_its = " << inner_its << std::endl
265 << "\tsg_operator_method = " << sg_op_names[opMethod]
266 << std::endl
267 << "\tsg_prec_method = " << sg_prec_names[precMethod]
268 << std::endl
269 << "\tgs_prec_tol = " << gs_prec_tol << std::endl
270 << "\tgs_prec_its = " << gs_prec_its << std::endl;
271 }
272
273 bool nonlinear_expansion = false;
274 if (randField == UNIFORM || randField == RYS)
275 nonlinear_expansion = false;
276 else if (randField == LOGNORMAL)
277 nonlinear_expansion = true;
278 bool scaleOP = true;
279
280 {
281 TEUCHOS_FUNC_TIME_MONITOR("Total PCE Calculation Time");
282
283 // Create Stochastic Galerkin basis and expansion
284 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(num_KL);
285 for (int i=0; i<num_KL; i++)
286 if (randField == UNIFORM)
287 bases[i] = Teuchos::rcp(new Stokhos::LegendreBasis<int,double>(p,normalize_basis));
288 else if (randField == RYS)
289 bases[i] = Teuchos::rcp(new Stokhos::RysBasis<int,double>(p,weightCut,normalize_basis));
290 else if (randField == LOGNORMAL)
291 bases[i] = Teuchos::rcp(new Stokhos::HermiteBasis<int,double>(p,normalize_basis));
292
293 // bases[i] = Teuchos::rcp(new Stokhos::DiscretizedStieltjesBasis<int,double>("beta",p,&uniform_weight,-weightCut,weightCut,true));
294 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
295 Teuchos::rcp(new Stokhos::CompletePolynomialBasis<int,double>(bases));
296 int sz = basis->size();
297 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk;
298 if (nonlinear_expansion)
299 Cijk = basis->computeTripleProductTensor();
300 else
301 Cijk = basis->computeLinearTripleProductTensor();
302 Teuchos::RCP<Stokhos::OrthogPolyExpansion<int,double> > expansion =
304 Cijk));
305 if (MyPID == 0)
306 std::cout << "Stochastic Galerkin expansion size = " << sz << std::endl;
307
308 // Create stochastic parallel distribution
309 Teuchos::ParameterList parallelParams;
310 parallelParams.set("Number of Spatial Processors", num_spatial_procs);
311 parallelParams.set("Rebalance Stochastic Graph",
312 rebalance_stochastic_graph);
313 Teuchos::RCP<Stokhos::ParallelData> sg_parallel_data =
314 Teuchos::rcp(new Stokhos::ParallelData(basis, Cijk, globalComm,
315 parallelParams));
316 Teuchos::RCP<const EpetraExt::MultiComm> sg_comm =
317 sg_parallel_data->getMultiComm();
318 Teuchos::RCP<const Epetra_Comm> app_comm =
319 sg_parallel_data->getSpatialComm();
320
321 // Create application
322 Teuchos::RCP<twoD_diffusion_ME> model =
323 Teuchos::rcp(new twoD_diffusion_ME(app_comm, n, num_KL, sigma,
324 mean, basis, nonlinear_expansion,
325 symmetric));
326
327 // Set up NOX parameters
328 Teuchos::RCP<Teuchos::ParameterList> noxParams =
329 Teuchos::rcp(new Teuchos::ParameterList);
330
331 // Set the nonlinear solver method
332 noxParams->set("Nonlinear Solver", "Line Search Based");
333
334 // Set the printing parameters in the "Printing" sublist
335 Teuchos::ParameterList& printParams = noxParams->sublist("Printing");
336 printParams.set("MyPID", MyPID);
337 printParams.set("Output Precision", 3);
338 printParams.set("Output Processor", 0);
339 printParams.set("Output Information",
340 NOX::Utils::OuterIteration +
341 NOX::Utils::OuterIterationStatusTest +
342 NOX::Utils::InnerIteration +
343 //NOX::Utils::Parameters +
344 NOX::Utils::Details +
345 NOX::Utils::LinearSolverDetails +
346 NOX::Utils::Warning +
347 NOX::Utils::Error);
348
349 // Create printing utilities
350 NOX::Utils utils(printParams);
351
352 // Sublist for line search
353 Teuchos::ParameterList& searchParams = noxParams->sublist("Line Search");
354 searchParams.set("Method", "Full Step");
355
356 // Sublist for direction
357 Teuchos::ParameterList& dirParams = noxParams->sublist("Direction");
358 dirParams.set("Method", "Newton");
359 Teuchos::ParameterList& newtonParams = dirParams.sublist("Newton");
360 newtonParams.set("Forcing Term Method", "Constant");
361
362 // Sublist for linear solver for the Newton method
363 Teuchos::ParameterList& lsParams = newtonParams.sublist("Linear Solver");
364
365 // Alternative linear solver list for Stratimikos
366 Teuchos::ParameterList& stratLinSolParams =
367 newtonParams.sublist("Stratimikos Linear Solver");
368 // Teuchos::ParameterList& noxStratParams =
369 // stratLinSolParams.sublist("NOX Stratimikos Options");
370 Teuchos::ParameterList& stratParams =
371 stratLinSolParams.sublist("Stratimikos");
372
373 // Sublist for convergence tests
374 Teuchos::ParameterList& statusParams = noxParams->sublist("Status Tests");
375 statusParams.set("Test Type", "Combo");
376 statusParams.set("Number of Tests", 2);
377 statusParams.set("Combo Type", "OR");
378 Teuchos::ParameterList& normF = statusParams.sublist("Test 0");
379 normF.set("Test Type", "NormF");
380 normF.set("Tolerance", outer_tol);
381 normF.set("Scale Type", "Scaled");
382 Teuchos::ParameterList& maxIters = statusParams.sublist("Test 1");
383 maxIters.set("Test Type", "MaxIters");
384 maxIters.set("Maximum Iterations", 1);
385
386 // Create NOX interface
387 Teuchos::RCP<NOX::Epetra::ModelEvaluatorInterface> det_nox_interface =
388 Teuchos::rcp(new NOX::Epetra::ModelEvaluatorInterface(model));
389
390 // Create NOX linear system object
391 Teuchos::RCP<const Epetra_Vector> det_u = model->get_x_init();
392 Teuchos::RCP<Epetra_Operator> det_A = model->create_W();
393 Teuchos::RCP<NOX::Epetra::Interface::Required> det_iReq = det_nox_interface;
394 Teuchos::RCP<NOX::Epetra::Interface::Jacobian> det_iJac = det_nox_interface;
395 Teuchos::ParameterList det_printParams;
396 det_printParams.set("MyPID", MyPID);
397 det_printParams.set("Output Precision", 3);
398 det_printParams.set("Output Processor", 0);
399 det_printParams.set("Output Information", NOX::Utils::Error);
400
401 Teuchos::ParameterList det_lsParams;
402 Teuchos::ParameterList& det_stratParams =
403 det_lsParams.sublist("Stratimikos");
404 if (inner_krylov_solver == AZTECOO) {
405 det_stratParams.set("Linear Solver Type", "AztecOO");
406 Teuchos::ParameterList& aztecOOParams =
407 det_stratParams.sublist("Linear Solver Types").sublist("AztecOO").sublist("Forward Solve");
408 Teuchos::ParameterList& aztecOOSettings =
409 aztecOOParams.sublist("AztecOO Settings");
410 if (inner_krylov_method == GMRES) {
411 aztecOOSettings.set("Aztec Solver","GMRES");
412 }
413 else if (inner_krylov_method == CG) {
414 aztecOOSettings.set("Aztec Solver","CG");
415 }
416 aztecOOSettings.set("Output Frequency", 0);
417 aztecOOSettings.set("Size of Krylov Subspace", 100);
418 aztecOOParams.set("Max Iterations", inner_its);
419 aztecOOParams.set("Tolerance", inner_tol);
420 Teuchos::ParameterList& verbParams =
421 det_stratParams.sublist("Linear Solver Types").sublist("AztecOO").sublist("VerboseObject");
422 verbParams.set("Verbosity Level", "none");
423 }
424 else if (inner_krylov_solver == BELOS) {
425 det_stratParams.set("Linear Solver Type", "Belos");
426 Teuchos::ParameterList& belosParams =
427 det_stratParams.sublist("Linear Solver Types").sublist("Belos");
428 Teuchos::ParameterList* belosSolverParams = NULL;
429 if (inner_krylov_method == GMRES || inner_krylov_method == FGMRES) {
430 belosParams.set("Solver Type","Block GMRES");
431 belosSolverParams =
432 &(belosParams.sublist("Solver Types").sublist("Block GMRES"));
433 if (inner_krylov_method == FGMRES)
434 belosSolverParams->set("Flexible Gmres", true);
435 }
436 else if (inner_krylov_method == CG) {
437 belosParams.set("Solver Type","Block CG");
438 belosSolverParams =
439 &(belosParams.sublist("Solver Types").sublist("Block CG"));
440 }
441 else if (inner_krylov_method == RGMRES) {
442 belosParams.set("Solver Type","GCRODR");
443 belosSolverParams =
444 &(belosParams.sublist("Solver Types").sublist("GCRODR"));
445 }
446 belosSolverParams->set("Convergence Tolerance", inner_tol);
447 belosSolverParams->set("Maximum Iterations", inner_its);
448 belosSolverParams->set("Output Frequency",0);
449 belosSolverParams->set("Output Style",1);
450 belosSolverParams->set("Verbosity",0);
451 Teuchos::ParameterList& verbParams = belosParams.sublist("VerboseObject");
452 verbParams.set("Verbosity Level", "none");
453 }
454 det_stratParams.set("Preconditioner Type", "ML");
455 Teuchos::ParameterList& det_ML =
456 det_stratParams.sublist("Preconditioner Types").sublist("ML").sublist("ML Settings");
457 ML_Epetra::SetDefaults("SA", det_ML);
458 det_ML.set("ML output", 0);
459 det_ML.set("max levels",5);
460 det_ML.set("increasing or decreasing","increasing");
461 det_ML.set("aggregation: type", "Uncoupled");
462 det_ML.set("smoother: type","ML symmetric Gauss-Seidel");
463 det_ML.set("smoother: sweeps",2);
464 det_ML.set("smoother: pre or post", "both");
465 det_ML.set("coarse: max size", 200);
466#ifdef HAVE_ML_AMESOS
467 det_ML.set("coarse: type","Amesos-KLU");
468#else
469 det_ML.set("coarse: type","Jacobi");
470#endif
471 Teuchos::RCP<NOX::Epetra::LinearSystem> det_linsys =
472 Teuchos::rcp(new NOX::Epetra::LinearSystemStratimikos(
473 det_printParams, det_lsParams, det_iJac,
474 det_A, *det_u));
475
476 // Setup stochastic Galerkin algorithmic parameters
477 Teuchos::RCP<Teuchos::ParameterList> sgParams =
478 Teuchos::rcp(new Teuchos::ParameterList);
479 Teuchos::ParameterList& sgOpParams =
480 sgParams->sublist("SG Operator");
481 Teuchos::ParameterList& sgPrecParams =
482 sgParams->sublist("SG Preconditioner");
483
484 if (!nonlinear_expansion) {
485 sgParams->set("Parameter Expansion Type", "Linear");
486 sgParams->set("Jacobian Expansion Type", "Linear");
487 }
488 if (opMethod == MATRIX_FREE)
489 sgOpParams.set("Operator Method", "Matrix Free");
490 else if (opMethod == KL_MATRIX_FREE)
491 sgOpParams.set("Operator Method", "KL Matrix Free");
492 else if (opMethod == KL_REDUCED_MATRIX_FREE) {
493 sgOpParams.set("Operator Method", "KL Reduced Matrix Free");
494 if (randField == UNIFORM || randField == RYS)
495 sgOpParams.set("Number of KL Terms", num_KL);
496 else
497 sgOpParams.set("Number of KL Terms", basis->size());
498 sgOpParams.set("KL Tolerance", outer_tol);
499 sgOpParams.set("Sparse 3 Tensor Drop Tolerance", outer_tol);
500 sgOpParams.set("Do Error Tests", true);
501 }
502 else if (opMethod == FULLY_ASSEMBLED)
503 sgOpParams.set("Operator Method", "Fully Assembled");
504 else
505 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
506 "Error! Unknown operator method " << opMethod
507 << "." << std::endl);
508 if (precMethod == MEAN) {
509 sgPrecParams.set("Preconditioner Method", "Mean-based");
510 sgPrecParams.set("Mean Preconditioner Type", "ML");
511 Teuchos::ParameterList& precParams =
512 sgPrecParams.sublist("Mean Preconditioner Parameters");
513 precParams = det_ML;
514 }
515 else if(precMethod == GS) {
516 sgPrecParams.set("Preconditioner Method", "Gauss-Seidel");
517 sgPrecParams.sublist("Deterministic Solver Parameters") = det_lsParams;
518 sgPrecParams.set("Deterministic Solver", det_linsys);
519 sgPrecParams.set("Max Iterations", gs_prec_its);
520 sgPrecParams.set("Tolerance", gs_prec_tol);
521 }
522 else if (precMethod == AGS) {
523 sgPrecParams.set("Preconditioner Method", "Approximate Gauss-Seidel");
524 if (outer_krylov_method == CG)
525 sgPrecParams.set("Symmetric Gauss-Seidel", true);
526 sgPrecParams.set("Mean Preconditioner Type", "ML");
527 Teuchos::ParameterList& precParams =
528 sgPrecParams.sublist("Mean Preconditioner Parameters");
529 precParams = det_ML;
530 }
531 else if (precMethod == AJ) {
532 sgPrecParams.set("Preconditioner Method", "Approximate Jacobi");
533 sgPrecParams.set("Mean Preconditioner Type", "ML");
534 Teuchos::ParameterList& precParams =
535 sgPrecParams.sublist("Mean Preconditioner Parameters");
536 precParams = det_ML;
537 Teuchos::ParameterList& jacobiOpParams =
538 sgPrecParams.sublist("Jacobi SG Operator");
539 jacobiOpParams.set("Only Use Linear Terms", true);
540 }
541 else if (precMethod == ASC) {
542 sgPrecParams.set("Preconditioner Method", "Approximate Schur Complement");
543 sgPrecParams.set("Mean Preconditioner Type", "ML");
544 Teuchos::ParameterList& precParams =
545 sgPrecParams.sublist("Mean Preconditioner Parameters");
546 precParams = det_ML;
547 }
548 else if (precMethod == KP) {
549 sgPrecParams.set("Preconditioner Method", "Kronecker Product");
550 sgPrecParams.set("Only Use Linear Terms", true);
551 sgPrecParams.set("Mean Preconditioner Type", "ML");
552 Teuchos::ParameterList& meanPrecParams =
553 sgPrecParams.sublist("Mean Preconditioner Parameters");
554 meanPrecParams = det_ML;
555 sgPrecParams.set("G Preconditioner Type", "Ifpack");
556 Teuchos::ParameterList& GPrecParams =
557 sgPrecParams.sublist("G Preconditioner Parameters");
558 if (outer_krylov_method == GMRES || outer_krylov_method == FGMRES)
559 GPrecParams.set("Ifpack Preconditioner", "ILUT");
560 if (outer_krylov_method == CG)
561 GPrecParams.set("Ifpack Preconditioner", "ICT");
562 GPrecParams.set("Overlap", 1);
563 GPrecParams.set("fact: drop tolerance", 1e-4);
564 GPrecParams.set("fact: ilut level-of-fill", 1.0);
565 GPrecParams.set("schwarz: combine mode", "Add");
566 }
567 else if (precMethod == NONE) {
568 sgPrecParams.set("Preconditioner Method", "None");
569 }
570 else
571 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
572 "Error! Unknown preconditioner method " << precMethod
573 << "." << std::endl);
574
575 // Create stochastic Galerkin model evaluator
576 Teuchos::RCP<Stokhos::SGModelEvaluator> sg_model =
577 Teuchos::rcp(new Stokhos::SGModelEvaluator(model, basis, Teuchos::null,
578 expansion, sg_parallel_data,
579 sgParams, scaleOP));
580 EpetraExt::ModelEvaluator::InArgs sg_inArgs = sg_model->createInArgs();
581 EpetraExt::ModelEvaluator::OutArgs sg_outArgs =
582 sg_model->createOutArgs();
583
584 // Set up stochastic parameters
585 // The current implementation of the model doesn't actually use these
586 // values, but is hard-coded to certain uncertainty models
587 Teuchos::Array<double> point(num_KL, 1.0);
588 Teuchos::Array<double> basis_vals(sz);
589 basis->evaluateBases(point, basis_vals);
590 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_p_init =
591 sg_model->create_p_sg(0);
592 for (int i=0; i<num_KL; i++) {
593 sg_p_init->term(i,0)[i] = 0.0;
594 sg_p_init->term(i,1)[i] = 1.0 / basis_vals[i+1];
595 }
596 sg_model->set_p_sg_init(0, *sg_p_init);
597
598 // Setup stochastic initial guess
599 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_x_init =
600 sg_model->create_x_sg();
601 sg_x_init->init(0.0);
602 sg_model->set_x_sg_init(*sg_x_init);
603
604 // Create NOX interface
605 Teuchos::RCP<NOX::Epetra::ModelEvaluatorInterface> nox_interface =
606 Teuchos::rcp(new NOX::Epetra::ModelEvaluatorInterface(sg_model));
607
608 // Create NOX stochastic linear system object
609 Teuchos::RCP<const Epetra_Vector> u = sg_model->get_x_init();
610 Teuchos::RCP<const Epetra_Map> base_map = model->get_x_map();
611 Teuchos::RCP<const Epetra_Map> sg_map = sg_model->get_x_map();
612 Teuchos::RCP<Epetra_Operator> A = sg_model->create_W();
613 Teuchos::RCP<NOX::Epetra::Interface::Required> iReq = nox_interface;
614 Teuchos::RCP<NOX::Epetra::Interface::Jacobian> iJac = nox_interface;
615
616 // Build linear solver
617 Teuchos::RCP<NOX::Epetra::LinearSystem> linsys;
618 if (solve_method==SG_KRYLOV) {
619 bool has_M =
620 sg_outArgs.supports(EpetraExt::ModelEvaluator::OUT_ARG_WPrec);
621 Teuchos::RCP<Epetra_Operator> M;
622 Teuchos::RCP<NOX::Epetra::Interface::Preconditioner> iPrec;
623 if (has_M) {
624 M = sg_model->create_WPrec()->PrecOp;
625 iPrec = nox_interface;
626 }
627 stratParams.set("Preconditioner Type", "None");
628 if (outer_krylov_solver == AZTECOO) {
629 stratParams.set("Linear Solver Type", "AztecOO");
630 Teuchos::ParameterList& aztecOOParams =
631 stratParams.sublist("Linear Solver Types").sublist("AztecOO").sublist("Forward Solve");
632 Teuchos::ParameterList& aztecOOSettings =
633 aztecOOParams.sublist("AztecOO Settings");
634 if (outer_krylov_method == GMRES) {
635 aztecOOSettings.set("Aztec Solver","GMRES");
636 }
637 else if (outer_krylov_method == CG) {
638 aztecOOSettings.set("Aztec Solver","CG");
639 }
640 aztecOOSettings.set("Output Frequency", 1);
641 aztecOOSettings.set("Size of Krylov Subspace", 100);
642 aztecOOParams.set("Max Iterations", outer_its);
643 aztecOOParams.set("Tolerance", outer_tol);
644 stratLinSolParams.set("Preconditioner", "User Defined");
645 if (has_M)
646 linsys =
647 Teuchos::rcp(new NOX::Epetra::LinearSystemStratimikos(
648 printParams, stratLinSolParams, iJac, A, iPrec, M,
649 *u, true));
650 else
651 linsys =
652 Teuchos::rcp(new NOX::Epetra::LinearSystemStratimikos(
653 printParams, stratLinSolParams, iJac, A, *u));
654 }
655 else if (outer_krylov_solver == BELOS){
656 stratParams.set("Linear Solver Type", "Belos");
657 Teuchos::ParameterList& belosParams =
658 stratParams.sublist("Linear Solver Types").sublist("Belos");
659 Teuchos::ParameterList* belosSolverParams = NULL;
660 if (outer_krylov_method == GMRES || outer_krylov_method == FGMRES) {
661 belosParams.set("Solver Type","Block GMRES");
662 belosSolverParams =
663 &(belosParams.sublist("Solver Types").sublist("Block GMRES"));
664 if (outer_krylov_method == FGMRES)
665 belosSolverParams->set("Flexible Gmres", true);
666 }
667 else if (outer_krylov_method == CG) {
668 belosParams.set("Solver Type","Block CG");
669 belosSolverParams =
670 &(belosParams.sublist("Solver Types").sublist("Block CG"));
671 }
672 else if (inner_krylov_method == RGMRES) {
673 belosParams.set("Solver Type","GCRODR");
674 belosSolverParams =
675 &(belosParams.sublist("Solver Types").sublist("GCRODR"));
676 }
677 belosSolverParams->set("Convergence Tolerance", outer_tol);
678 belosSolverParams->set("Maximum Iterations", outer_its);
679 belosSolverParams->set("Output Frequency",1);
680 belosSolverParams->set("Output Style",1);
681 belosSolverParams->set("Verbosity",33);
682 stratLinSolParams.set("Preconditioner", "User Defined");
683 if (has_M)
684 linsys =
685 Teuchos::rcp(new NOX::Epetra::LinearSystemStratimikos(
686 printParams, stratLinSolParams, iJac, A, iPrec, M,
687 *u, true));
688 else
689 linsys =
690 Teuchos::rcp(new NOX::Epetra::LinearSystemStratimikos(
691 printParams, stratLinSolParams, iJac, A, *u));
692
693 }
694 }
695 else if (solve_method==SG_GS) {
696 lsParams.sublist("Deterministic Solver Parameters") = det_lsParams;
697 lsParams.set("Max Iterations", outer_its);
698 lsParams.set("Tolerance", outer_tol);
699 linsys =
700 Teuchos::rcp(new NOX::Epetra::LinearSystemSGGS(
701 printParams, lsParams, det_linsys, iReq, iJac,
702 basis, sg_parallel_data, A, base_map, sg_map));
703 }
704 else {
705 lsParams.sublist("Deterministic Solver Parameters") = det_lsParams;
706 lsParams.set("Max Iterations", outer_its);
707 lsParams.set("Tolerance", outer_tol);
708 Teuchos::ParameterList& jacobiOpParams =
709 lsParams.sublist("Jacobi SG Operator");
710 jacobiOpParams.set("Only Use Linear Terms", true);
711 linsys =
712 Teuchos::rcp(new NOX::Epetra::LinearSystemSGJacobi(
713 printParams, lsParams, det_linsys, iReq, iJac,
714 basis, sg_parallel_data, A, base_map, sg_map));
715 }
716
717 // Build NOX group
718 Teuchos::RCP<NOX::Epetra::Group> grp =
719 Teuchos::rcp(new NOX::Epetra::Group(printParams, iReq, *u, linsys));
720
721 // Create the Solver convergence test
722 Teuchos::RCP<NOX::StatusTest::Generic> statusTests =
723 NOX::StatusTest::buildStatusTests(statusParams, utils);
724
725 // Create the solver
726 Teuchos::RCP<NOX::Solver::Generic> solver =
727 NOX::Solver::buildSolver(grp, statusTests, noxParams);
728
729 // Solve the system
730 NOX::StatusTest::StatusType status;
731 {
732 TEUCHOS_FUNC_TIME_MONITOR("Total Solve Time");
733 status = solver->solve();
734 }
735
736 // Get final solution
737 const NOX::Epetra::Group& finalGroup =
738 dynamic_cast<const NOX::Epetra::Group&>(solver->getSolutionGroup());
739 const Epetra_Vector& finalSolution =
740 (dynamic_cast<const NOX::Epetra::Vector&>(finalGroup.getX())).getEpetraVector();
741
742 // Save final solution to file
743 EpetraExt::VectorToMatrixMarketFile("nox_solver_stochastic_solution.mm",
744 finalSolution);
745
746 // Save mean and variance to file
747 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_x_poly =
748 sg_model->create_x_sg(View, &finalSolution);
749 Epetra_Vector mean(*(model->get_x_map()));
750 Epetra_Vector std_dev(*(model->get_x_map()));
751 sg_x_poly->computeMean(mean);
752 sg_x_poly->computeStandardDeviation(std_dev);
753 EpetraExt::VectorToMatrixMarketFile("mean_gal.mm", mean);
754 EpetraExt::VectorToMatrixMarketFile("std_dev_gal.mm", std_dev);
755
756 // Evaluate SG responses at SG parameters
757 Teuchos::RCP<const Epetra_Vector> sg_p = sg_model->get_p_init(1);
758 Teuchos::RCP<Epetra_Vector> sg_g =
759 Teuchos::rcp(new Epetra_Vector(*(sg_model->get_g_map(0))));
760 sg_inArgs.set_p(1, sg_p);
761 sg_inArgs.set_x(Teuchos::rcp(&finalSolution,false));
762 sg_outArgs.set_g(0, sg_g);
763 sg_model->evalModel(sg_inArgs, sg_outArgs);
764
765 // Print mean and standard deviation of response
766 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_g_poly =
767 sg_model->create_g_sg(0, View, sg_g.get());
768 Epetra_Vector g_mean(*(model->get_g_map(0)));
769 Epetra_Vector g_std_dev(*(model->get_g_map(0)));
770 sg_g_poly->computeMean(g_mean);
771 sg_g_poly->computeStandardDeviation(g_std_dev);
772 std::cout.precision(16);
773 // std::cout << "\nResponse Expansion = " << std::endl;
774 // std::cout.precision(12);
775 // sg_g_poly->print(std::cout);
776 std::cout << std::endl;
777 std::cout << "Response Mean = " << std::endl << g_mean << std::endl;
778 std::cout << "Response Std. Dev. = " << std::endl << g_std_dev << std::endl;
779
780 if (status == NOX::StatusTest::Converged && MyPID == 0)
781 utils.out() << "Example Passed!" << std::endl;
782
783 }
784
785 Teuchos::TimeMonitor::summarize(std::cout);
786 Teuchos::TimeMonitor::zeroOutTimers();
787
788 }
789
790 catch (std::exception& e) {
791 std::cout << e.what() << std::endl;
792 }
793 catch (std::string& s) {
794 std::cout << s << std::endl;
795 }
796 catch (char *s) {
797 std::cout << s << std::endl;
798 }
799 catch (...) {
800 std::cout << "Caught unknown exception!" << std::endl;
801 }
802
803#ifdef HAVE_MPI
804 MPI_Finalize() ;
805#endif
806
807}
Orthogonal polynomial expansions limited to algebraic operations.
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Hermite polynomial basis.
Legendre polynomial basis.
Rys polynomial basis.
Nonlinear, stochastic Galerkin ModelEvaluator.
ModelEvaluator for a linear 2-D diffusion problem.
const SG_Solver sg_solver_values[]
int main(int argc, char *argv[])
double uniform_weight(const double &x)
const char * sg_solver_names[]
const char * krylov_method_names[]
const Krylov_Method krylov_method_values[]
const char * krylov_solver_names[]
const char * sg_prec_names[]
const SG_Prec sg_prec_values[]
const char * sg_op_names[]
const Krylov_Solver krylov_solver_values[]