124 MPI_Init(&argc,&
argv);
128 Teuchos::RCP<const Epetra_Comm> globalComm;
134 int MyPID = globalComm->MyPID();
139 Teuchos::CommandLineProcessor
CLP;
141 "This example runs a variety of stochastic Galerkin solvers.\n");
144 CLP.setOption(
"num_mesh", &n,
"Number of mesh points in each direction");
146 bool symmetric =
false;
147 CLP.setOption(
"symmetric",
"unsymmetric", &symmetric,
148 "Symmetric discretization");
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)");
153 bool rebalance_stochastic_graph =
false;
154 CLP.setOption(
"rebalance",
"no-rebalance", &rebalance_stochastic_graph,
155 "Rebalance parallel stochastic graph (requires Isorropia)");
158 CLP.setOption(
"rand_field", &randField,
160 "Random field type");
163 CLP.setOption(
"mean", &mean,
"Mean");
166 CLP.setOption(
"std_dev", &sigma,
"Standard deviation");
168 double weightCut = 1.0;
169 CLP.setOption(
"weight_cut", &weightCut,
"Weight cut");
172 CLP.setOption(
"num_kl", &num_KL,
"Number of KL terms");
175 CLP.setOption(
"order", &p,
"Polynomial order");
177 bool normalize_basis =
true;
178 CLP.setOption(
"normalize",
"unnormalize", &normalize_basis,
179 "Normalize PC basis");
182 CLP.setOption(
"sg_solver", &solve_method,
187 CLP.setOption(
"outer_krylov_method", &outer_krylov_method,
189 "Outer Krylov method (for Krylov-based SG solver)");
192 CLP.setOption(
"outer_krylov_solver", &outer_krylov_solver,
194 "Outer linear solver");
196 double outer_tol = 1e-12;
197 CLP.setOption(
"outer_tol", &outer_tol,
"Outer solver tolerance");
199 int outer_its = 1000;
200 CLP.setOption(
"outer_its", &outer_its,
"Maximum outer iterations");
203 CLP.setOption(
"inner_krylov_method", &inner_krylov_method,
205 "Inner Krylov method (for G-S, Jacobi, etc...)");
208 CLP.setOption(
"inner_krylov_solver", &inner_krylov_solver,
210 "Inner linear solver");
212 double inner_tol = 3e-13;
213 CLP.setOption(
"inner_tol", &inner_tol,
"Inner solver tolerance");
215 int inner_its = 1000;
216 CLP.setOption(
"inner_its", &inner_its,
"Maximum inner iterations");
219 CLP.setOption(
"sg_operator_method", &opMethod,
224 CLP.setOption(
"sg_prec_method", &precMethod,
226 "Preconditioner method");
228 double gs_prec_tol = 1e-1;
229 CLP.setOption(
"gs_prec_tol", &gs_prec_tol,
"Gauss-Seidel preconditioner tolerance");
232 CLP.setOption(
"gs_prec_its", &gs_prec_its,
"Maximum Gauss-Seidel preconditioner iterations");
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
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
253 <<
"\touter_krylov_method = "
255 <<
"\touter_krylov_solver = "
257 <<
"\touter_tol = " << outer_tol << std::endl
258 <<
"\touter_its = " << outer_its << std::endl
259 <<
"\tinner_krylov_method = "
261 <<
"\tinner_krylov_solver = "
263 <<
"\tinner_tol = " << inner_tol << std::endl
264 <<
"\tinner_its = " << inner_its << std::endl
265 <<
"\tsg_operator_method = " <<
sg_op_names[opMethod]
269 <<
"\tgs_prec_tol = " << gs_prec_tol << std::endl
270 <<
"\tgs_prec_its = " << gs_prec_its << std::endl;
273 bool nonlinear_expansion =
false;
275 nonlinear_expansion =
false;
277 nonlinear_expansion =
true;
281 TEUCHOS_FUNC_TIME_MONITOR(
"Total PCE Calculation Time");
284 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(num_KL);
285 for (
int i=0; i<num_KL; i++)
288 else if (randField ==
RYS)
294 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
296 int sz = basis->size();
297 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk;
298 if (nonlinear_expansion)
299 Cijk = basis->computeTripleProductTensor();
301 Cijk = basis->computeLinearTripleProductTensor();
302 Teuchos::RCP<Stokhos::OrthogPolyExpansion<int,double> > expansion =
306 std::cout <<
"Stochastic Galerkin expansion size = " << sz << std::endl;
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 =
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();
322 Teuchos::RCP<twoD_diffusion_ME> model =
324 mean, basis, nonlinear_expansion,
328 Teuchos::RCP<Teuchos::ParameterList> noxParams =
329 Teuchos::rcp(
new Teuchos::ParameterList);
332 noxParams->set(
"Nonlinear Solver",
"Line Search Based");
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 +
344 NOX::Utils::Details +
345 NOX::Utils::LinearSolverDetails +
346 NOX::Utils::Warning +
350 NOX::Utils utils(printParams);
353 Teuchos::ParameterList& searchParams = noxParams->sublist(
"Line Search");
354 searchParams.set(
"Method",
"Full Step");
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");
363 Teuchos::ParameterList& lsParams = newtonParams.sublist(
"Linear Solver");
366 Teuchos::ParameterList& stratLinSolParams =
367 newtonParams.sublist(
"Stratimikos Linear Solver");
370 Teuchos::ParameterList& stratParams =
371 stratLinSolParams.sublist(
"Stratimikos");
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);
387 Teuchos::RCP<NOX::Epetra::ModelEvaluatorInterface> det_nox_interface =
388 Teuchos::rcp(
new NOX::Epetra::ModelEvaluatorInterface(model));
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);
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");
413 else if (inner_krylov_method ==
CG) {
414 aztecOOSettings.set(
"Aztec Solver",
"CG");
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");
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");
432 &(belosParams.sublist(
"Solver Types").sublist(
"Block GMRES"));
433 if (inner_krylov_method ==
FGMRES)
434 belosSolverParams->set(
"Flexible Gmres",
true);
436 else if (inner_krylov_method ==
CG) {
437 belosParams.set(
"Solver Type",
"Block CG");
439 &(belosParams.sublist(
"Solver Types").sublist(
"Block CG"));
441 else if (inner_krylov_method ==
RGMRES) {
442 belosParams.set(
"Solver Type",
"GCRODR");
444 &(belosParams.sublist(
"Solver Types").sublist(
"GCRODR"));
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");
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);
467 det_ML.set(
"coarse: type",
"Amesos-KLU");
469 det_ML.set(
"coarse: type",
"Jacobi");
471 Teuchos::RCP<NOX::Epetra::LinearSystem> det_linsys =
472 Teuchos::rcp(
new NOX::Epetra::LinearSystemStratimikos(
473 det_printParams, det_lsParams, det_iJac,
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");
484 if (!nonlinear_expansion) {
485 sgParams->set(
"Parameter Expansion Type",
"Linear");
486 sgParams->set(
"Jacobian Expansion Type",
"Linear");
489 sgOpParams.set(
"Operator Method",
"Matrix Free");
491 sgOpParams.set(
"Operator Method",
"KL Matrix Free");
493 sgOpParams.set(
"Operator Method",
"KL Reduced Matrix Free");
495 sgOpParams.set(
"Number of KL Terms", num_KL);
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);
503 sgOpParams.set(
"Operator Method",
"Fully Assembled");
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");
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);
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");
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");
537 Teuchos::ParameterList& jacobiOpParams =
538 sgPrecParams.sublist(
"Jacobi SG Operator");
539 jacobiOpParams.set(
"Only Use Linear Terms",
true);
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");
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");
567 else if (precMethod ==
NONE) {
568 sgPrecParams.set(
"Preconditioner Method",
"None");
571 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
572 "Error! Unknown preconditioner method " << precMethod
573 <<
"." << std::endl);
576 Teuchos::RCP<Stokhos::SGModelEvaluator> sg_model =
578 expansion, sg_parallel_data,
580 EpetraExt::ModelEvaluator::InArgs sg_inArgs = sg_model->createInArgs();
581 EpetraExt::ModelEvaluator::OutArgs sg_outArgs =
582 sg_model->createOutArgs();
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];
596 sg_model->set_p_sg_init(0, *sg_p_init);
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);
605 Teuchos::RCP<NOX::Epetra::ModelEvaluatorInterface> nox_interface =
606 Teuchos::rcp(
new NOX::Epetra::ModelEvaluatorInterface(sg_model));
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;
617 Teuchos::RCP<NOX::Epetra::LinearSystem> linsys;
620 sg_outArgs.supports(EpetraExt::ModelEvaluator::OUT_ARG_WPrec);
621 Teuchos::RCP<Epetra_Operator> M;
622 Teuchos::RCP<NOX::Epetra::Interface::Preconditioner> iPrec;
624 M = sg_model->create_WPrec()->PrecOp;
625 iPrec = nox_interface;
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");
637 else if (outer_krylov_method ==
CG) {
638 aztecOOSettings.set(
"Aztec Solver",
"CG");
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");
647 Teuchos::rcp(
new NOX::Epetra::LinearSystemStratimikos(
648 printParams, stratLinSolParams, iJac, A, iPrec, M,
652 Teuchos::rcp(
new NOX::Epetra::LinearSystemStratimikos(
653 printParams, stratLinSolParams, iJac, A, *u));
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");
663 &(belosParams.sublist(
"Solver Types").sublist(
"Block GMRES"));
664 if (outer_krylov_method ==
FGMRES)
665 belosSolverParams->set(
"Flexible Gmres",
true);
667 else if (outer_krylov_method ==
CG) {
668 belosParams.set(
"Solver Type",
"Block CG");
670 &(belosParams.sublist(
"Solver Types").sublist(
"Block CG"));
672 else if (inner_krylov_method ==
RGMRES) {
673 belosParams.set(
"Solver Type",
"GCRODR");
675 &(belosParams.sublist(
"Solver Types").sublist(
"GCRODR"));
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");
685 Teuchos::rcp(
new NOX::Epetra::LinearSystemStratimikos(
686 printParams, stratLinSolParams, iJac, A, iPrec, M,
690 Teuchos::rcp(
new NOX::Epetra::LinearSystemStratimikos(
691 printParams, stratLinSolParams, iJac, A, *u));
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);
700 Teuchos::rcp(
new NOX::Epetra::LinearSystemSGGS(
701 printParams, lsParams, det_linsys, iReq, iJac,
702 basis, sg_parallel_data, A, base_map, sg_map));
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);
712 Teuchos::rcp(
new NOX::Epetra::LinearSystemSGJacobi(
713 printParams, lsParams, det_linsys, iReq, iJac,
714 basis, sg_parallel_data, A, base_map, sg_map));
718 Teuchos::RCP<NOX::Epetra::Group> grp =
719 Teuchos::rcp(
new NOX::Epetra::Group(printParams, iReq, *u, linsys));
722 Teuchos::RCP<NOX::StatusTest::Generic> statusTests =
723 NOX::StatusTest::buildStatusTests(statusParams, utils);
726 Teuchos::RCP<NOX::Solver::Generic> solver =
727 NOX::Solver::buildSolver(grp, statusTests, noxParams);
730 NOX::StatusTest::StatusType status;
732 TEUCHOS_FUNC_TIME_MONITOR(
"Total Solve Time");
733 status = solver->solve();
737 const NOX::Epetra::Group& finalGroup =
738 dynamic_cast<const NOX::Epetra::Group&
>(solver->getSolutionGroup());
740 (
dynamic_cast<const NOX::Epetra::Vector&
>(finalGroup.getX())).getEpetraVector();
743 EpetraExt::VectorToMatrixMarketFile(
"nox_solver_stochastic_solution.mm",
747 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_x_poly =
748 sg_model->create_x_sg(
View, &finalSolution);
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);
757 Teuchos::RCP<const Epetra_Vector> sg_p = sg_model->get_p_init(1);
758 Teuchos::RCP<Epetra_Vector> sg_g =
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);
766 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_g_poly =
767 sg_model->create_g_sg(0,
View, sg_g.get());
770 sg_g_poly->computeMean(g_mean);
771 sg_g_poly->computeStandardDeviation(g_std_dev);
772 std::cout.precision(16);
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;
780 if (status == NOX::StatusTest::Converged && MyPID == 0)
781 utils.out() <<
"Example Passed!" << std::endl;
785 Teuchos::TimeMonitor::summarize(std::cout);
786 Teuchos::TimeMonitor::zeroOutTimers();
790 catch (std::exception& e) {
791 std::cout << e.what() << std::endl;
793 catch (std::string& s) {
794 std::cout << s << std::endl;
797 std::cout << s << std::endl;
800 std::cout <<
"Caught unknown exception!" << std::endl;