47 #include "Epetra_MpiComm.h" 49 #include "Epetra_SerialComm.h" 54 #include "NOX_Epetra.H" 55 #include "NOX_Epetra_LinearSystem_Stratimikos.H" 56 #include "BelosTypes.hpp" 63 #include "Teuchos_TimeMonitor.hpp" 66 #include "EpetraExt_VectorOut.h" 67 #include "EpetraExt_BlockUtility.h" 75 bool nonlinear_expansion =
false;
77 bool symmetric =
true;
78 bool use_solver =
true;
79 std::string solver_type =
"RGMRES";
83 MPI_Init(&argc,&
argv);
90 Epetra_Object::SetTracebackMode(1);
93 TEUCHOS_FUNC_TIME_MONITOR(
"Total PCE Calculation Time");
96 Teuchos::RCP<const Epetra_Comm> globalComm;
98 globalComm = Teuchos::rcp(
new Epetra_MpiComm(MPI_COMM_WORLD));
100 globalComm = Teuchos::rcp(
new Epetra_SerialComm);
102 MyPID = globalComm->MyPID();
105 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(num_KL);
106 for (
int i=0; i<num_KL; i++)
108 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<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();
123 std::cout <<
"Stochastic Galerkin expansion size = " << sz << std::endl;
126 int num_spatial_procs = -1;
128 num_spatial_procs = std::atoi(
argv[1]);
129 Teuchos::RCP<const EpetraExt::MultiComm> multi_comm =
131 Teuchos::RCP<const Epetra_Comm> mp_comm =
133 Teuchos::RCP<const Epetra_Comm> app_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();
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");
157 precParams.set(
"coarse: type",
"Jacobi");
161 Teuchos::RCP<twoD_diffusion_ME> model =
163 nonlinear_expansion, symmetric,
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;
178 Teuchos::RCP<Stokhos::MPModelEvaluator> mp_model =
180 mp_block_map, mpParams));
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];
191 mp_model->set_p_mp_init(0, *mp_p_init);
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);
200 Teuchos::RCP<Teuchos::ParameterList> noxParams =
201 Teuchos::rcp(
new Teuchos::ParameterList);
204 noxParams->set(
"Nonlinear Solver",
"Line Search Based");
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 +
220 NOX::Utils utils(printParams);
223 Teuchos::ParameterList& searchParams = noxParams->sublist(
"Line Search");
224 searchParams.set(
"Method",
"Full Step");
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");
233 Teuchos::ParameterList stratLinSolParams;
234 Teuchos::ParameterList& stratParams =
235 stratLinSolParams.sublist(
"Stratimikos");
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");
245 &(belosParams.sublist(
"Solver Types").sublist(
"Block GMRES"));
247 else if (solver_type ==
"CG") {
248 belosParams.set(
"Solver Type",
"Block CG");
250 &(belosParams.sublist(
"Solver Types").sublist(
"Block CG"));
252 else if (solver_type ==
"RGMRES") {
253 belosParams.set(
"Solver Type",
"GCRODR");
255 &(belosParams.sublist(
"Solver Types").sublist(
"GCRODR"));
256 belosSolverParams->set(
"Num Recycled Blocks", 20);
258 else if (solver_type ==
"RCG") {
259 belosParams.set(
"Solver Type",
"RCG");
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);
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);
277 belosSolverParams->set(
"Verbosity",
280 Belos::StatusTestDetails);
281 stratLinSolParams.set(
"Preconditioner",
"User Defined");
282 Teuchos::ParameterList& verboseParams =
283 belosParams.sublist(
"VerboseObject");
284 verboseParams.set(
"Verbosity Level",
"medium");
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);
300 Teuchos::RCP<NOX::Epetra::ModelEvaluatorInterface> nox_interface =
301 Teuchos::rcp(
new NOX::Epetra::ModelEvaluatorInterface(mp_model));
304 Teuchos::RCP<const Epetra_Vector> u = mp_model->get_x_init();
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;
313 Teuchos::ParameterList& outerSolParams =
314 newtonParams.sublist(
"Linear Solver");
315 outerSolParams.sublist(
"Deterministic Solver Parameters") =
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 =
324 Teuchos::RCP<NOX::Epetra::Interface::Jacobian> inner_iJac =
326 Teuchos::RCP<NOX::Epetra::Interface::Preconditioner> inner_iPrec =
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(
333 inner_iJac, inner_A, inner_iPrec, inner_M,
336 Teuchos::rcp(
new NOX::Epetra::LinearSystemMPBD(printParams,
340 model->get_x_map()));
343 newtonParams.sublist(
"Stratimikos Linear Solver") = stratLinSolParams;
345 Teuchos::rcp(
new NOX::Epetra::LinearSystemStratimikos(printParams,
352 Teuchos::RCP<NOX::Epetra::Group> grp =
353 Teuchos::rcp(
new NOX::Epetra::Group(printParams, iReq, *u, linsys));
356 Teuchos::RCP<NOX::StatusTest::Generic> statusTests =
357 NOX::StatusTest::buildStatusTests(statusParams, utils);
360 Teuchos::RCP<NOX::Solver::Generic> solver =
361 NOX::Solver::buildSolver(grp, statusTests, noxParams);
364 NOX::StatusTest::StatusType status = solver->solve();
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);
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);
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);
406 multi_comm, View, mp_local_x);
411 multi_comm, View, mp_local_g);
412 Teuchos::RCP<Epetra_LocalMap> sg_block_map =
413 Teuchos::rcp(
new Epetra_LocalMap(sz, 0, *globalComm));
415 model->get_x_map(), multi_comm);
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);
428 EpetraExt::VectorToMatrixMarketFile(
"ni_stochastic_solution.mm",
429 *(sg_x.getBlockVector()));
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);
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);
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;
450 if (status == NOX::StatusTest::Converged && MyPID == 0)
451 utils.out() <<
"Example Passed!" << std::endl;
455 Teuchos::TimeMonitor::summarize(std::cout);
456 Teuchos::TimeMonitor::zeroOutTimers();
460 catch (std::exception& e) {
461 std::cout << e.what() << std::endl;
463 catch (std::string& s) {
464 std::cout << s << std::endl;
467 std::cout << s << std::endl;
470 std::cout <<
"Caught unknown exception!" <<std:: endl;
int main(int argc, char *argv[])
Teuchos::RCP< const Epetra_Comm > getStochasticComm(const Teuchos::RCP< const EpetraExt::MultiComm > &globalMultiComm)
Multi-point model evaluator.
Teuchos::RCP< const EpetraExt::MultiComm > buildMultiComm(const Epetra_Comm &globalComm, int num_global_stochastic_blocks, int num_spatial_procs=-1)
ModelEvaluator for a linear 2-D diffusion problem.
A container class storing an orthogonal polynomial whose coefficients are vectors, operators, or in general any type that would have an expensive copy constructor.
A container class for products of Epetra_Vector's.
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Legendre polynomial basis.
Teuchos::RCP< const Epetra_Comm > getSpatialComm(const Teuchos::RCP< const EpetraExt::MultiComm > &globalMultiComm)