47 #include "Epetra_MpiComm.h" 49 #include "Epetra_SerialComm.h" 59 #include "Teuchos_CommandLineProcessor.hpp" 60 #include "Teuchos_TimeMonitor.hpp" 63 #include "EpetraExt_VectorOut.h" 81 const char *
sg_rf_names[] = {
"Uniform",
"CC-Uniform",
"Rys",
"Log-Normal" };
95 "Slow Restricted",
"Moderate Restricted",
"Unrestricted" };
101 MPI_Init(&argc,&
argv);
105 Teuchos::RCP<Epetra_Comm> Comm;
107 Comm = Teuchos::rcp(
new Epetra_MpiComm(MPI_COMM_WORLD));
109 Comm = Teuchos::rcp(
new Epetra_SerialComm);
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 =
264 Teuchos::rcp(
new Epetra_Vector(*(model.
get_p_map(0))));
265 Teuchos::RCP<Epetra_Vector>
x =
266 Teuchos::rcp(
new Epetra_Vector(*(model.
get_x_map())));
267 Teuchos::RCP<Epetra_Vector>
f =
268 Teuchos::rcp(
new Epetra_Vector(*(model.
get_f_map())));
269 Teuchos::RCP<Epetra_Vector>
g =
270 Teuchos::rcp(
new Epetra_Vector(*(model.
get_g_map(0))));
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();
279 Epetra_Vector x_mean(*(model.
get_x_map()));
280 Epetra_Vector x_var(*(model.
get_x_map()));
282 Epetra_Vector g_mean(*(model.
get_g_map(0)));
283 Epetra_Vector g_var(*(model.
get_g_map(0)));
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);
296 #ifdef HAVE_ML_AMESOS 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");
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];
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);
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++)
387 for (
int i=0; i<g_var.MyLength(); 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;
const char * prec_strategy_names[]
KOKKOS_INLINE_FUNCTION PCE< Storage > sqrt(const PCE< Storage > &a)
Hermite polynomial basis.
const Krylov_Method krylov_method_values[]
const SG_Quad sg_quad_values[]
Teuchos::RCP< Epetra_CrsMatrix > get_mean() const
Get mean matrix.
const int num_prec_strategy
Teuchos::RCP< const Epetra_Map > get_f_map() const
Return residual vector map.
const SG_GROWTH sg_growth_values[]
int main(int argc, char *argv[])
const char * sg_growth_names[]
Teuchos::RCP< const Epetra_Map > get_g_map(int j) const
Return response function map.
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
ModelEvaluator for a linear 2-D diffusion problem.
const int num_krylov_method
const SG_RF sg_rf_values[]
const char * sg_quad_names[]
Teuchos::RCP< const Epetra_Map > get_x_map() const
Return solution vector map.
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Legendre polynomial basis.
void evalModel(const InArgs &inArgs, const OutArgs &outArgs) const
Evaluate model on InArgs.
const char * sg_rf_names[]
ScalarType f(const Teuchos::Array< ScalarType > &x, double a, double b)
Legendre polynomial basis using Clenshaw-Curtis quadrature points.
Teuchos::RCP< Epetra_Operator > create_W() const
Create W = alpha*M + beta*J matrix.
InArgs createInArgs() const
Create InArgs.
const PrecStrategy prec_strategy_values[]
Teuchos::RCP< const Epetra_Map > get_p_map(int l) const
Return parameter vector map.
OutArgs createOutArgs() const
Create OutArgs.
Defines quadrature for a tensor product basis by tensor products of 1-D quadrature rules...
ScalarType g(const Teuchos::Array< ScalarType > &x, const ScalarType &y)
const char * krylov_method_names[]