101 MPI_Init(&argc,&
argv);
105 Teuchos::RCP<Epetra_Comm> Comm;
112 int MyPID = Comm->MyPID();
117 Teuchos::CommandLineProcessor
CLP;
119 "This example runs a stochastic collocation method.\n");
122 CLP.setOption(
"num_mesh", &n,
"Number of mesh points in each direction");
125 CLP.setOption(
"rand_field", &randField,
127 "Random field type");
130 CLP.setOption(
"mean", &mean,
"Mean");
133 CLP.setOption(
"std_dev", &sigma,
"Standard deviation");
135 double weightCut = 1.0;
136 CLP.setOption(
"weight_cut", &weightCut,
"Weight cut");
139 CLP.setOption(
"num_kl", &num_KL,
"Number of KL terms");
142 CLP.setOption(
"order", &p,
"Polynomial order");
144 bool normalize_basis =
true;
145 CLP.setOption(
"normalize",
"unnormalize", &normalize_basis,
146 "Normalize PC basis");
149 CLP.setOption(
"krylov_method", &krylov_method,
154 CLP.setOption(
"prec_strategy", &precStrategy,
156 "Preconditioner strategy");
159 CLP.setOption(
"tol", &tol,
"Solver tolerance");
161#ifdef HAVE_STOKHOS_DAKOTA
166 CLP.setOption(
"quadrature", &quad_method,
171 CLP.setOption(
"sg_growth", &sg_growth,
173 "Sparse grid growth rule");
178 std::cout <<
"Summary of command line options:" << std::endl
179 <<
"\tnum_mesh = " << n << std::endl
180 <<
"\trand_field = " <<
sg_rf_names[randField] << std::endl
181 <<
"\tmean = " << mean << std::endl
182 <<
"\tstd_dev = " << sigma << std::endl
183 <<
"\tnum_kl = " << num_KL << std::endl
184 <<
"\torder = " << p << std::endl
185 <<
"\tnormalize_basis = " << normalize_basis << std::endl
190 <<
"\ttol = " << tol << std::endl
197 bool nonlinear_expansion =
false;
199 nonlinear_expansion =
true;
202 TEUCHOS_FUNC_TIME_MONITOR(
"Total Collocation Calculation Time");
205 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(num_KL);
206 for (
int i=0; i<num_KL; i++) {
207 Teuchos::RCP<Stokhos::OneDOrthogPolyBasis<int,double> > b;
210 p, normalize_basis));
215 p, normalize_basis,
true));
217 else if (randField ==
RYS) {
219 p, weightCut, normalize_basis));
223 p, normalize_basis));
227 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
231 Teuchos::RCP<Stokhos::Quadrature<int,double> > quad;
233#ifdef HAVE_STOKHOS_DAKOTA
234 int sparse_grid_growth = Pecos::MODERATE_RESTRICTED_GROWTH;
236 sparse_grid_growth = Pecos::SLOW_RESTRICTED_GROWTH;
238 sparse_grid_growth = Pecos::MODERATE_RESTRICTED_GROWTH;
240 sparse_grid_growth = Pecos::UNRESTRICTED_GROWTH;
241 quad = Teuchos::rcp(
new Stokhos::SparseGridQuadrature<int,double>(
242 basis,p,1e-12,sparse_grid_growth));
244 TEUCHOS_TEST_FOR_EXCEPTION(
245 true, std::logic_error,
246 "Sparse grids require building with Dakota support!");
249 else if (quad_method ==
TENSOR) {
253 const Teuchos::Array< Teuchos::Array<double> >& quad_points =
254 quad->getQuadPoints();
255 const Teuchos::Array<double>& quad_weights = quad->getQuadWeights();
256 int nqp = quad_weights.size();
260 basis, nonlinear_expansion);
263 Teuchos::RCP<Epetra_Vector> p =
265 Teuchos::RCP<Epetra_Vector> x =
267 Teuchos::RCP<Epetra_Vector>
f =
269 Teuchos::RCP<Epetra_Vector>
g =
271 Teuchos::RCP<Epetra_CrsMatrix> A =
272 Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(model.create_W());
273 EpetraExt::ModelEvaluator::InArgs inArgs = model.createInArgs();
274 EpetraExt::ModelEvaluator::OutArgs outArgs = model.createOutArgs();
275 EpetraExt::ModelEvaluator::OutArgs outArgs2 = model.createOutArgs();
286 Teuchos::ParameterList precParams;
287 precParams.set(
"default values",
"SA");
288 precParams.set(
"ML output", 0);
289 precParams.set(
"max levels",5);
290 precParams.set(
"increasing or decreasing",
"increasing");
291 precParams.set(
"aggregation: type",
"Uncoupled");
292 precParams.set(
"smoother: type",
"ML symmetric Gauss-Seidel");
293 precParams.set(
"smoother: sweeps",2);
294 precParams.set(
"smoother: pre or post",
"both");
295 precParams.set(
"coarse: max size", 200);
297 precParams.set(
"coarse: type",
"Amesos-KLU");
299 precParams.set(
"coarse: type",
"Jacobi");
301 bool checkFiltering =
false;
302 if (precStrategy ==
REUSE) {
303 checkFiltering =
true;
304 precParams.set(
"reuse: enable",
true);
306 Teuchos::RCP<ML_Epetra::MultiLevelPreconditioner> M =
307 Teuchos::rcp(
new ML_Epetra::MultiLevelPreconditioner(*A, precParams,
309 if (precStrategy ==
MEAN) {
310 TEUCHOS_FUNC_TIME_MONITOR(
"Deterministic Preconditioner Calculation");
311 *A = *(model.get_mean());
312 M->ComputePreconditioner();
317 if (krylov_method ==
GMRES)
318 aztec.SetAztecOption(AZ_solver, AZ_gmres);
319 else if (krylov_method ==
CG)
320 aztec.SetAztecOption(AZ_solver, AZ_cg);
321 aztec.SetAztecOption(AZ_precond, AZ_none);
322 aztec.SetAztecOption(AZ_kspace, 100);
323 aztec.SetAztecOption(AZ_conv, AZ_r0);
324 aztec.SetAztecOption(AZ_output, 0);
325 aztec.SetUserOperator(A.get());
326 aztec.SetPrecOperator(M.get());
327 aztec.SetLHS(x.get());
328 aztec.SetRHS(
f.get());
330 x_mean.PutScalar(0.0);
331 x_var.PutScalar(0.0);
333 for (
int qp=0; qp<nqp; qp++) {
334 TEUCHOS_FUNC_TIME_MONITOR(
"Collocation Loop");
336 if (qp%100 == 0 || qp == nqp-1)
337 std::cout <<
"Collocation point " << qp+1 <<
'/' << nqp <<
"\n";
340 for (
int i=0; i<num_KL; i++)
341 (*p)[i] = quad_points[qp][i];
351 model.evalModel(inArgs, outArgs);
355 if (precStrategy !=
MEAN) {
356 TEUCHOS_FUNC_TIME_MONITOR(
"Deterministic Preconditioner Calculation");
357 M->ComputePreconditioner(checkFiltering);
362 TEUCHOS_FUNC_TIME_MONITOR(
"Deterministic Solve");
363 aztec.Iterate(1000, tol);
367 outArgs2.set_g(0,
g);
368 model.evalModel(inArgs, outArgs2);
371 x2.Multiply(1.0, *x, *x, 0.0);
372 g2.Multiply(1.0, *
g, *
g, 0.0);
373 x_mean.Update(quad_weights[qp], *x, 1.0);
374 x_var.Update(quad_weights[qp], x2, 1.0);
375 g_mean.Update(quad_weights[qp], *
g, 1.0);
376 g_var.Update(quad_weights[qp], g2, 1.0);
379 x2.Multiply(1.0, x_mean, x_mean, 0.0);
380 g2.Multiply(1.0, g_mean, g_mean, 0.0);
381 x_var.Update(-1.0, x2, 1.0);
382 g_var.Update(-1.0, g2, 1.0);
385 for (
int i=0; i<x_var.MyLength(); i++)
386 x_var[i] = std::sqrt(x_var[i]);
387 for (
int i=0; i<g_var.MyLength(); i++)
388 g_var[i] = std::sqrt(g_var[i]);
390 std::cout.precision(16);
391 std::cout <<
"\nResponse Mean = " << std::endl << g_mean << std::endl;
392 std::cout <<
"Response Std. Dev. = " << std::endl << g_var << std::endl;
395 EpetraExt::VectorToMatrixMarketFile(
"mean_col.mm", x_mean);
396 EpetraExt::VectorToMatrixMarketFile(
"std_dev_col.mm", x_var);
400 Teuchos::TimeMonitor::summarize(std::cout);
401 Teuchos::TimeMonitor::zeroOutTimers();
405 catch (std::exception& e) {
406 std::cout << e.what() << std::endl;
408 catch (std::string& s) {
409 std::cout << s << std::endl;
412 std::cout << s << std::endl;
415 std::cout <<
"Caught unknown exception!" <<std:: endl;