Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
linear2d_diffusion_pce_mpni.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"
56#include "BelosTypes.hpp"
57
58// Stokhos Stochastic Galerkin
59#include "Stokhos_Epetra.hpp"
61
62// Timing utilities
63#include "Teuchos_TimeMonitor.hpp"
64
65// I/O utilities
66#include "EpetraExt_VectorOut.h"
67#include "EpetraExt_BlockUtility.h"
68
69int main(int argc, char *argv[]) {
70 int n = 32; // spatial discretization (per dimension)
71 int num_KL = 2; // number of KL terms
72 int p = 3; // polynomial order
73 double mu = 0.1; // mean of exponential random field
74 double s = 0.2; // std. dev. of exponential r.f.
75 bool nonlinear_expansion = false; // nonlinear expansion of diffusion coeff
76 // (e.g., log-normal)
77 bool symmetric = true; // use symmetric formulation
78 bool use_solver = true;
79 std::string solver_type = "RGMRES";
80
81// Initialize MPI
82#ifdef HAVE_MPI
83 MPI_Init(&argc,&argv);
84#endif
85
86 int MyPID;
87
88 try {
89
91
92 {
93 TEUCHOS_FUNC_TIME_MONITOR("Total PCE Calculation Time");
94
95 // Create a communicator for Epetra objects
96 Teuchos::RCP<const Epetra_Comm> globalComm;
97#ifdef HAVE_MPI
98 globalComm = Teuchos::rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
99#else
100 globalComm = Teuchos::rcp(new Epetra_SerialComm);
101#endif
102 MyPID = globalComm->MyPID();
103
104 // Create Stochastic Galerkin basis and quadrature
105 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(num_KL);
106 for (int i=0; i<num_KL; i++)
107 bases[i] = Teuchos::rcp(new Stokhos::LegendreBasis<int,double>(p, true));
108 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
110 1e-12));
111 // Teuchos::RCP<const Stokhos::Quadrature<int,double> > quad =
112 // Teuchos::rcp(new Stokhos::TensorProductQuadrature<int,double>(basis));
113 Teuchos::RCP<const Stokhos::Quadrature<int,double> > quad =
114 Teuchos::rcp(new Stokhos::SparseGridQuadrature<int,double>(basis));
115 const Teuchos::Array< Teuchos::Array<double> >& quad_points =
116 quad->getQuadPoints();
117 const Teuchos::Array<double>& quad_weights = quad->getQuadWeights();
118 const Teuchos::Array< Teuchos::Array<double> >& quad_values =
119 quad->getBasisAtQuadPoints();
120 int sz = basis->size();
121 int num_mp = quad_weights.size();
122 if (MyPID == 0)
123 std::cout << "Stochastic Galerkin expansion size = " << sz << std::endl;
124
125 // Create multi-point parallel distribution
126 int num_spatial_procs = -1;
127 if (argc > 1)
128 num_spatial_procs = std::atoi(argv[1]);
129 Teuchos::RCP<const EpetraExt::MultiComm> multi_comm =
130 Stokhos::buildMultiComm(*globalComm, num_mp, num_spatial_procs);
131 Teuchos::RCP<const Epetra_Comm> mp_comm =
132 Stokhos::getStochasticComm(multi_comm);
133 Teuchos::RCP<const Epetra_Comm> app_comm =
134 Stokhos::getSpatialComm(multi_comm);
135 Teuchos::RCP<const Epetra_Map> mp_block_map =
136 Teuchos::rcp(new Epetra_Map(num_mp, 0, *mp_comm));
137 int num_my_mp = mp_block_map->NumMyElements();
138 //mp_block_map->Print(std::cout);
139
140 Teuchos::RCP<Teuchos::ParameterList> detPrecParams =
141 Teuchos::rcp(new Teuchos::ParameterList);
142 detPrecParams->set("Preconditioner Type", "ML");
143 Teuchos::ParameterList& precParams =
144 detPrecParams->sublist("Preconditioner Parameters");
145 precParams.set("default values", "SA");
146 precParams.set("ML output", 0);
147 precParams.set("max levels",5);
148 precParams.set("increasing or decreasing","increasing");
149 precParams.set("aggregation: type", "Uncoupled");
150 precParams.set("smoother: type","ML symmetric Gauss-Seidel");
151 precParams.set("smoother: sweeps",2);
152 precParams.set("smoother: pre or post", "both");
153 precParams.set("coarse: max size", 200);
154#ifdef HAVE_ML_AMESOS
155 precParams.set("coarse: type","Amesos-KLU");
156#else
157 precParams.set("coarse: type","Jacobi");
158#endif
159
160 // Create application
161 Teuchos::RCP<twoD_diffusion_ME> model =
162 Teuchos::rcp(new twoD_diffusion_ME(app_comm, n, num_KL, mu, s, basis,
163 nonlinear_expansion, symmetric,
164 detPrecParams));
165
166 // Setup multi-point algorithmic parameters
167 Teuchos::RCP<Teuchos::ParameterList> mpParams =
168 Teuchos::rcp(new Teuchos::ParameterList);
169 Teuchos::ParameterList& mpPrecParams =
170 mpParams->sublist("MP Preconditioner");
171 mpPrecParams.set("Preconditioner Method", "Block Diagonal");
172 mpPrecParams.set("MP Preconditioner Type", "ML");
173 Teuchos::ParameterList& pointPrecParams =
174 mpPrecParams.sublist("MP Preconditioner Parameters");
175 pointPrecParams = precParams;
176
177 // Create stochastic Galerkin model evaluator
178 Teuchos::RCP<Stokhos::MPModelEvaluator> mp_model =
179 Teuchos::rcp(new Stokhos::MPModelEvaluator(model, multi_comm,
180 mp_block_map, mpParams));
181
182 // Set up multi-point parameters
183 Teuchos::RCP<Stokhos::ProductEpetraVector> mp_p_init =
184 mp_model->create_p_mp(0);
185 int my_mp_begin = mp_block_map->MinMyGID();
186 for (int j=0; j<num_my_mp; j++) {
187 for (int i=0; i<num_KL; i++) {
188 (*mp_p_init)[j][i] = quad_points[j+my_mp_begin][i];
189 }
190 }
191 mp_model->set_p_mp_init(0, *mp_p_init);
192
193 // Setup multi-point initial guess
194 Teuchos::RCP<Stokhos::ProductEpetraVector> mp_x_init =
195 mp_model->create_x_mp();
196 mp_x_init->init(0.0);
197 mp_model->set_x_mp_init(*mp_x_init);
198
199 // Set up NOX parameters
200 Teuchos::RCP<Teuchos::ParameterList> noxParams =
201 Teuchos::rcp(new Teuchos::ParameterList);
202
203 // Set the nonlinear solver method
204 noxParams->set("Nonlinear Solver", "Line Search Based");
205
206 // Set the printing parameters in the "Printing" sublist
207 Teuchos::ParameterList& printParams = noxParams->sublist("Printing");
208 printParams.set("MyPID", MyPID);
209 printParams.set("Output Precision", 3);
210 printParams.set("Output Processor", 0);
211 printParams.set("Output Information",
212 NOX::Utils::OuterIteration +
213 NOX::Utils::OuterIterationStatusTest +
214 NOX::Utils::InnerIteration +
215 NOX::Utils::LinearSolverDetails +
216 NOX::Utils::Warning +
217 NOX::Utils::Error);
218
219 // Create printing utilities
220 NOX::Utils utils(printParams);
221
222 // Sublist for line search
223 Teuchos::ParameterList& searchParams = noxParams->sublist("Line Search");
224 searchParams.set("Method", "Full Step");
225
226 // Sublist for direction
227 Teuchos::ParameterList& dirParams = noxParams->sublist("Direction");
228 dirParams.set("Method", "Newton");
229 Teuchos::ParameterList& newtonParams = dirParams.sublist("Newton");
230 newtonParams.set("Forcing Term Method", "Constant");
231
232 // Alternative linear solver list for Stratimikos
233 Teuchos::ParameterList stratLinSolParams;
234 Teuchos::ParameterList& stratParams =
235 stratLinSolParams.sublist("Stratimikos");
236
237 // Sublist for linear solver for the Newton method
238 stratParams.set("Linear Solver Type", "Belos");
239 Teuchos::ParameterList& belosParams =
240 stratParams.sublist("Linear Solver Types").sublist("Belos");
241 Teuchos::ParameterList* belosSolverParams = NULL;
242 if (solver_type == "GMRES") {
243 belosParams.set("Solver Type","Block GMRES");
244 belosSolverParams =
245 &(belosParams.sublist("Solver Types").sublist("Block GMRES"));
246 }
247 else if (solver_type == "CG") {
248 belosParams.set("Solver Type","Block CG");
249 belosSolverParams =
250 &(belosParams.sublist("Solver Types").sublist("Block CG"));
251 }
252 else if (solver_type == "RGMRES") {
253 belosParams.set("Solver Type","GCRODR");
254 belosSolverParams =
255 &(belosParams.sublist("Solver Types").sublist("GCRODR"));
256 belosSolverParams->set("Num Recycled Blocks", 20);
257 }
258 else if (solver_type == "RCG") {
259 belosParams.set("Solver Type","RCG");
260 belosSolverParams =
261 &(belosParams.sublist("Solver Types").sublist("RCG"));
262 Teuchos::RCP<const Teuchos::ParameterList> ortho_params =
263 Teuchos::rcp(new Teuchos::ParameterList);
264 belosSolverParams->set("Num Recycled Blocks", 10);
265 }
266 else
267 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
268 "Unknown solver type " << solver_type);
269 belosSolverParams->set("Convergence Tolerance", 1e-12);
270 belosSolverParams->set("Maximum Iterations", 1000);
271 if (solver_type != "CG")
272 belosSolverParams->set("Num Blocks", 100);
273 belosSolverParams->set("Block Size", 1);
274 belosSolverParams->set("Output Frequency",1);
275 belosSolverParams->set("Output Style",1);
276 //belosSolverParams->set("Verbosity",33);
277 belosSolverParams->set("Verbosity",
278 Belos::Errors +
279 Belos::Warnings +
280 Belos::StatusTestDetails);
281 stratLinSolParams.set("Preconditioner", "User Defined");
282 Teuchos::ParameterList& verboseParams =
283 belosParams.sublist("VerboseObject");
284 verboseParams.set("Verbosity Level", "medium");
285
286 // Sublist for convergence tests
287 Teuchos::ParameterList& statusParams = noxParams->sublist("Status Tests");
288 statusParams.set("Test Type", "Combo");
289 statusParams.set("Number of Tests", 2);
290 statusParams.set("Combo Type", "OR");
291 Teuchos::ParameterList& normF = statusParams.sublist("Test 0");
292 normF.set("Test Type", "NormF");
293 normF.set("Tolerance", 1e-10);
294 normF.set("Scale Type", "Scaled");
295 Teuchos::ParameterList& maxIters = statusParams.sublist("Test 1");
296 maxIters.set("Test Type", "MaxIters");
297 maxIters.set("Maximum Iterations", 1);
298
299 // Create NOX interface
300 Teuchos::RCP<NOX::Epetra::ModelEvaluatorInterface> nox_interface =
301 Teuchos::rcp(new NOX::Epetra::ModelEvaluatorInterface(mp_model));
302
303 // Create NOX linear system object
304 Teuchos::RCP<const Epetra_Vector> u = mp_model->get_x_init();
305
306 Teuchos::RCP<NOX::Epetra::LinearSystem> linsys;
307 Teuchos::RCP<Epetra_Operator> A = mp_model->create_W();
308 Teuchos::RCP<Epetra_Operator> M = mp_model->create_WPrec()->PrecOp;
309 Teuchos::RCP<NOX::Epetra::Interface::Required> iReq = nox_interface;
310 Teuchos::RCP<NOX::Epetra::Interface::Jacobian> iJac = nox_interface;
311 Teuchos::RCP<NOX::Epetra::Interface::Preconditioner> iPrec = nox_interface;
312 if (use_solver) {
313 Teuchos::ParameterList& outerSolParams =
314 newtonParams.sublist("Linear Solver");
315 outerSolParams.sublist("Deterministic Solver Parameters") =
316 stratLinSolParams;
317 outerSolParams.set("Preconditioner Strategy", "Mean");
318 Teuchos::RCP<Epetra_Operator> inner_A = model->create_W();
319 Teuchos::RCP<Epetra_Operator> inner_M = model->create_WPrec()->PrecOp;
320 Teuchos::RCP<NOX::Epetra::ModelEvaluatorInterface> inner_nox_interface =
321 Teuchos::rcp(new NOX::Epetra::ModelEvaluatorInterface(model));
322 Teuchos::RCP<NOX::Epetra::Interface::Required> inner_iReq =
323 inner_nox_interface;
324 Teuchos::RCP<NOX::Epetra::Interface::Jacobian> inner_iJac =
325 inner_nox_interface;
326 Teuchos::RCP<NOX::Epetra::Interface::Preconditioner> inner_iPrec =
327 inner_nox_interface;
328 Teuchos::RCP<const Epetra_Vector> inner_u = model->get_x_init();
329 Teuchos::RCP<NOX::Epetra::LinearSystem> inner_linsys =
330 Teuchos::rcp(new NOX::Epetra::LinearSystemStratimikos(
331 printParams,
332 stratLinSolParams,
333 inner_iJac, inner_A, inner_iPrec, inner_M,
334 *inner_u, true));
335 linsys =
336 Teuchos::rcp(new NOX::Epetra::LinearSystemMPBD(printParams,
337 outerSolParams,
338 inner_linsys,
339 iReq, iJac, A,
340 model->get_x_map()));
341 }
342 else {
343 newtonParams.sublist("Stratimikos Linear Solver") = stratLinSolParams;
344 linsys =
345 Teuchos::rcp(new NOX::Epetra::LinearSystemStratimikos(printParams,
346 stratLinSolParams,
347 iJac, A, iPrec, M,
348 *u, true));
349 }
350
351 // Build NOX group
352 Teuchos::RCP<NOX::Epetra::Group> grp =
353 Teuchos::rcp(new NOX::Epetra::Group(printParams, iReq, *u, linsys));
354
355 // Create the Solver convergence test
356 Teuchos::RCP<NOX::StatusTest::Generic> statusTests =
357 NOX::StatusTest::buildStatusTests(statusParams, utils);
358
359 // Create the solver
360 Teuchos::RCP<NOX::Solver::Generic> solver =
361 NOX::Solver::buildSolver(grp, statusTests, noxParams);
362
363 // Solve the system
364 NOX::StatusTest::StatusType status = solver->solve();
365
366 // Get final solution
367 const NOX::Epetra::Group& finalGroup =
368 dynamic_cast<const NOX::Epetra::Group&>(solver->getSolutionGroup());
369 const Epetra_Vector& finalSolution =
370 (dynamic_cast<const NOX::Epetra::Vector&>(finalGroup.getX())).getEpetraVector();
371 Teuchos::RCP<const Epetra_Vector> mp_x =
372 Teuchos::rcp(&finalSolution, false);
373
374 // Evaluate final responses
375 Teuchos::RCP<const Epetra_Vector> mp_p = mp_model->get_p_init(1);
376 Teuchos::RCP<Epetra_Vector> mp_g =
377 Teuchos::rcp(new Epetra_Vector(*(mp_model->get_g_map(0))));
378 EpetraExt::ModelEvaluator::InArgs mp_inArgs = mp_model->createInArgs();
379 EpetraExt::ModelEvaluator::OutArgs mp_outArgs = mp_model->createOutArgs();
380 mp_inArgs.set_x(mp_x);
381 mp_inArgs.set_p(1, mp_p);
382 mp_outArgs.set_g(0, mp_g);
383 mp_model->evalModel(mp_inArgs, mp_outArgs);
384
385 // Import x and g
386 Teuchos::RCP<Epetra_LocalMap> mp_local_map =
387 Teuchos::rcp(new Epetra_LocalMap(num_mp, 0, *globalComm));
388 Teuchos::RCP<const Epetra_Map> mp_local_x_map =
389 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(
390 *(model->get_x_map()), *mp_local_map, *globalComm));
391 Teuchos::RCP<const Epetra_Map> mp_local_g_map =
392 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(
393 *(model->get_g_map(0)), *mp_local_map, *globalComm));
394 Epetra_Import mp_x_importer(*mp_local_x_map, mp_x->Map());
395 Epetra_Import mp_g_importer(*mp_local_g_map, mp_g->Map());
396 Epetra_Vector mp_local_x(*mp_local_x_map);
397 Epetra_Vector mp_local_g(*mp_local_g_map);
398 mp_local_x.Import(*mp_x, mp_x_importer, Add);
399 mp_local_g.Import(*mp_g, mp_g_importer, Add);
400
401 // Compute PC expansions
403 mp_local_map,
404 model->get_x_map(),
405 mp_local_x_map,
406 multi_comm, View, mp_local_x);
408 mp_local_map,
409 model->get_g_map(0),
410 mp_local_g_map,
411 multi_comm, View, mp_local_g);
412 Teuchos::RCP<Epetra_LocalMap> sg_block_map =
413 Teuchos::rcp(new Epetra_LocalMap(sz, 0, *globalComm));
414 Stokhos::EpetraVectorOrthogPoly sg_x(basis, sg_block_map,
415 model->get_x_map(), multi_comm);
416 Stokhos::EpetraVectorOrthogPoly sg_g(basis, sg_block_map,
417 model->get_g_map(0), multi_comm);
418 for (int i=0; i<sz; i++) {
419 sg_x[i].PutScalar(0.0);
420 sg_g[i].PutScalar(0.0);
421 for (int j=0; j<num_mp; j++) {
422 sg_x[i].Update(quad_weights[j]*quad_values[j][i], mp_x_vec[j], 1.0);
423 sg_g[i].Update(quad_weights[j]*quad_values[j][i], mp_g_vec[j], 1.0);
424 }
425 }
426
427 // Save SG solution to file
428 EpetraExt::VectorToMatrixMarketFile("ni_stochastic_solution.mm",
429 *(sg_x.getBlockVector()));
430
431 // Save mean and variance to file
432 Epetra_Vector mean(*(model->get_x_map()));
433 Epetra_Vector std_dev(*(model->get_x_map()));
434 sg_x.computeMean(mean);
435 sg_x.computeStandardDeviation(std_dev);
436 EpetraExt::VectorToMatrixMarketFile("mean_gal.mm", mean);
437
438 // Print mean and standard deviation of responses
439 Epetra_Vector g_mean(*(model->get_g_map(0)));
440 Epetra_Vector g_std_dev(*(model->get_g_map(0)));
441 sg_g.computeMean(g_mean);
442 sg_g.computeStandardDeviation(g_std_dev);
443 // std::cout << "\nResponse Expansion = " << std::endl;
444 // std::cout.precision(12);
445 // sg_g.print(std::cout);
446 std::cout << std::endl;
447 std::cout << "Response Mean = " << std::endl << g_mean << std::endl;
448 std::cout << "Response Std. Dev. = " << std::endl << g_std_dev << std::endl;
449
450 if (status == NOX::StatusTest::Converged && MyPID == 0)
451 utils.out() << "Example Passed!" << std::endl;
452
453 }
454
455 Teuchos::TimeMonitor::summarize(std::cout);
456 Teuchos::TimeMonitor::zeroOutTimers();
457
458 }
459
460 catch (std::exception& e) {
461 std::cout << e.what() << std::endl;
462 }
463 catch (std::string& s) {
464 std::cout << s << std::endl;
465 }
466 catch (char *s) {
467 std::cout << s << std::endl;
468 }
469 catch (...) {
470 std::cout << "Caught unknown exception!" <<std:: endl;
471 }
472
473#ifdef HAVE_MPI
474 MPI_Finalize() ;
475#endif
476
477}
Add
int Import(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
static void SetTracebackMode(int TracebackModeValue)
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
A container class storing an orthogonal polynomial whose coefficients are vectors,...
void computeMean(Epetra_Vector &v) const
Compute mean.
void computeStandardDeviation(Epetra_Vector &v) const
Compute standard deviation.
Legendre polynomial basis.
Multi-point model evaluator.
A container class for products of Epetra_Vector's.
Teuchos::RCP< EpetraExt::BlockVector > getBlockVector()
Get block vector.
ModelEvaluator for a linear 2-D diffusion problem.
int main(int argc, char *argv[])
Teuchos::RCP< const Epetra_Comm > getStochasticComm(const Teuchos::RCP< const EpetraExt::MultiComm > &globalMultiComm)
Teuchos::RCP< const Epetra_Comm > getSpatialComm(const Teuchos::RCP< const EpetraExt::MultiComm > &globalMultiComm)
Teuchos::RCP< const EpetraExt::MultiComm > buildMultiComm(const Epetra_Comm &globalComm, int num_global_stochastic_blocks, int num_spatial_procs=-1)