106 MPI_Init(&argc,&
argv);
110 Teuchos::RCP<Epetra_Comm> Comm;
117 int MyPID = Comm->MyPID();
122 Teuchos::CommandLineProcessor
CLP;
124 "This example runs a stochastic collocation method.\n");
127 CLP.setOption(
"num_mesh", &n,
"Number of mesh points in each direction");
130 CLP.setOption(
"rand_field", &randField,
132 "Random field type");
135 CLP.setOption(
"mean", &mean,
"Mean");
138 CLP.setOption(
"std_dev", &sigma,
"Standard deviation");
140 double weightCut = 1.0;
141 CLP.setOption(
"weight_cut", &weightCut,
"Weight cut");
144 CLP.setOption(
"num_kl", &num_KL,
"Number of KL terms");
147 CLP.setOption(
"order", &p,
"Polynomial order");
149 bool normalize_basis =
true;
150 CLP.setOption(
"normalize",
"unnormalize", &normalize_basis,
151 "Normalize PC basis");
154 CLP.setOption(
"krylov_solver", &krylov_solver,
159 CLP.setOption(
"krylov_method", &krylov_method,
164 CLP.setOption(
"prec_strategy", &precStrategy,
166 "Preconditioner strategy");
169 CLP.setOption(
"tol", &tol,
"Solver tolerance");
172 CLP.setOption(
"sg_growth", &sg_growth,
174 "Sparse grid growth rule");
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
194 <<
"\ttol = " << tol << std::endl
199 bool nonlinear_expansion =
false;
201 nonlinear_expansion =
true;
204 TEUCHOS_FUNC_TIME_MONITOR(
"Total Collocation Calculation Time");
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;
212 p, normalize_basis));
217 p, normalize_basis,
true));
219 else if (randField ==
RYS) {
221 p, weightCut, normalize_basis));
225 p, normalize_basis));
229 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
233 int sparse_grid_growth = Pecos::MODERATE_RESTRICTED_GROWTH;
235 sparse_grid_growth = Pecos::SLOW_RESTRICTED_GROWTH;
237 sparse_grid_growth = Pecos::MODERATE_RESTRICTED_GROWTH;
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();
249 basis, nonlinear_expansion);
252 Teuchos::RCP<Epetra_Vector> p =
254 Teuchos::RCP<Epetra_Vector> x =
256 Teuchos::RCP<Epetra_Vector>
f =
258 Teuchos::RCP<Epetra_Vector>
g =
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();
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");
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");
304 &(belosParams.sublist(
"Solver Types").sublist(
"Block GMRES"));
305 if (krylov_method ==
FGMRES)
306 belosSolverParams->set(
"Flexible Gmres",
true);
308 else if (krylov_method ==
RGMRES) {
309 belosParams.set(
"Solver Type",
"GCRODR");
311 &(belosParams.sublist(
"Solver Types").sublist(
"GCRODR"));
312 belosSolverParams->set(
"Num Recycled Blocks", 10);
314 else if (krylov_method ==
CG) {
315 belosParams.set(
"Solver Type",
"Block CG");
317 &(belosParams.sublist(
"Solver Types").sublist(
"Block CG"));
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");
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);
340 precParams.set(
"coarse: type",
"Amesos-KLU");
342 precParams.set(
"coarse: type",
"Jacobi");
346 Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
347 linearSolverBuilder.setParameterList(Teuchos::rcp(&stratParams,
false));
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);
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();
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));
373 Teuchos::RCP<Thyra::SolveCriteria<double> > solveCriteria;
374 if (!(krylov_solver ==
BELOS && krylov_method ==
CG)) {
376 Thyra::SolveMeasureType solveMeasure(
377 Thyra::SOLVE_MEASURE_NORM_RESIDUAL,
378 Thyra::SOLVE_MEASURE_NORM_INIT_RESIDUAL);
380 Teuchos::rcp(
new Thyra::SolveCriteria<double>(solveMeasure, tol));
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());
392 x_mean.PutScalar(0.0);
393 x_var.PutScalar(0.0);
395 for (
int qp=0; qp<nqp; qp++) {
396 TEUCHOS_FUNC_TIME_MONITOR(
"Collocation Loop");
399 for (
int i=0; i<num_KL; i++)
400 (*p)[i] = quad_points[qp][i];
410 model.evalModel(inArgs, outArgs);
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());
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());
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());
433 std::cout <<
"Collocation point " << qp+1 <<
'/' << nqp <<
": "
434 << solveStatus.message << std::endl;
438 x_thyra = Teuchos::null;
439 f_thyra = Teuchos::null;
442 outArgs2.set_g(0,
g);
443 model.evalModel(inArgs, outArgs2);
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);
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);
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]);
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;
470 EpetraExt::VectorToMatrixMarketFile(
"mean_col.mm", x_mean);
471 EpetraExt::VectorToMatrixMarketFile(
"std_dev_col.mm", x_var);
475 Teuchos::TimeMonitor::summarize(std::cout);
476 Teuchos::TimeMonitor::zeroOutTimers();
480 catch (std::exception& e) {
481 std::cout << e.what() << std::endl;
483 catch (std::string& s) {
484 std::cout << s << std::endl;
487 std::cout << s << std::endl;
490 std::cout <<
"Caught unknown exception!" <<std:: endl;