59 #include "Teuchos_CommandLineProcessor.hpp" 60 #include "Teuchos_TimeMonitor.hpp" 63 #include "EpetraExt_BlockUtility.h" 64 #include "EpetraExt_RowMatrixOut.h" 76 const char *
sg_rf_names[] = {
"Uniform",
"CC-Uniform",
"Rys",
"Log-Normal" };
90 "Slow Restricted",
"Moderate Restricted",
"Unrestricted" };
95 using Teuchos::rcp_dynamic_cast;
97 using Teuchos::CommandLineProcessor;
98 using Teuchos::TimeMonitor;
99 using Teuchos::ParameterList;
100 using EpetraExt::BlockUtility;
101 using EpetraExt::ModelEvaluator;
105 MPI_Init(&argc,&
argv);
109 RCP<Epetra_Comm> Comm;
116 int MyPID = Comm->MyPID();
121 CommandLineProcessor
CLP;
123 "This example runs a stochastic collocation method.\n");
126 CLP.setOption(
"num_mesh", &n,
"Number of mesh points in each direction");
129 CLP.setOption(
"rand_field", &randField,
131 "Random field type");
134 CLP.setOption(
"mean", &mean,
"Mean");
137 CLP.setOption(
"std_dev", &sigma,
"Standard deviation");
139 double weightCut = 1.0;
140 CLP.setOption(
"weight_cut", &weightCut,
"Weight cut");
143 CLP.setOption(
"num_kl", &num_KL,
"Number of KL terms");
146 CLP.setOption(
"order", &p,
"Polynomial order");
148 bool normalize_basis =
true;
149 CLP.setOption(
"normalize",
"unnormalize", &normalize_basis,
150 "Normalize PC basis");
153 CLP.setOption(
"krylov_method", &krylov_method,
158 CLP.setOption(
"tol", &tol,
"Solver tolerance");
160 #ifdef HAVE_STOKHOS_DAKOTA 165 CLP.setOption(
"quadrature", &quad_method,
170 CLP.setOption(
"sg_growth", &sg_growth,
172 "Sparse grid growth rule");
177 std::cout <<
"Summary of command line options:" << std::endl
178 <<
"\tnum_mesh = " << n << std::endl
179 <<
"\trand_field = " <<
sg_rf_names[randField] << std::endl
180 <<
"\tmean = " << mean << std::endl
181 <<
"\tstd_dev = " << sigma << std::endl
182 <<
"\tnum_kl = " << num_KL << std::endl
183 <<
"\torder = " << p << std::endl
184 <<
"\tnormalize_basis = " << normalize_basis << std::endl
187 <<
"\ttol = " << tol << std::endl
194 bool nonlinear_expansion =
false;
196 nonlinear_expansion =
true;
199 TEUCHOS_FUNC_TIME_MONITOR(
"Total Collocation Calculation Time");
202 Array< RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(num_KL);
203 for (
int i=0; i<num_KL; i++) {
204 RCP<Stokhos::OneDOrthogPolyBasis<int,double> > b;
211 p, normalize_basis,
true));
213 else if (randField ==
RYS) {
215 p, weightCut, normalize_basis));
219 p, normalize_basis));
223 RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
227 RCP<Stokhos::Quadrature<int,double> > quad;
229 #ifdef HAVE_STOKHOS_DAKOTA 230 int sparse_grid_growth = Pecos::MODERATE_RESTRICTED_GROWTH;
232 sparse_grid_growth = Pecos::SLOW_RESTRICTED_GROWTH;
234 sparse_grid_growth = Pecos::MODERATE_RESTRICTED_GROWTH;
236 sparse_grid_growth = Pecos::UNRESTRICTED_GROWTH;
237 quad = rcp(
new Stokhos::SparseGridQuadrature<int,double>(
238 basis,p,1e-12,sparse_grid_growth));
240 TEUCHOS_TEST_FOR_EXCEPTION(
241 true, std::logic_error,
242 "Sparse grids require building with Dakota support!");
245 else if (quad_method ==
TENSOR) {
248 const Array< Array<double> >& quad_points = quad->getQuadPoints();
249 const Array<double>& quad_weights = quad->getQuadWeights();
250 int nqp = quad_weights.size();
254 basis, nonlinear_expansion);
258 RCP<Epetra_CrsMatrix> A =
262 RCP<const Epetra_Map> spatial_map = model.
get_x_map();
264 RCP<const Epetra_Map> product_map =
265 rcp(BlockUtility::GenerateBlockMap(stochastic_map, *spatial_map, *Comm),
269 for (
int i=0; i<nqp; ++i)
272 RCP<const Epetra_CrsGraph> product_graph =
273 rcp(BlockUtility::GenerateBlockGraph(stochastic_graph, spatial_graph, *Comm));
279 std::cout <<
"spatial size = " << spatial_map->NumGlobalElements()
280 <<
" quadrature size = " << nqp
281 <<
" multipoint size = " << product_map->NumGlobalElements()
286 int max_num_entries = A->MaxNumEntries();
287 Array<int> block_indices(max_num_entries);
291 for (
int qp=0; qp<nqp; qp++) {
292 TEUCHOS_FUNC_TIME_MONITOR(
"Compute MP Matrix, RHS");
295 for (
int i=0; i<num_KL; i++)
296 (*p)[i] = quad_points[qp][i];
311 const int num_rows = spatial_map->NumMyElements();
312 for (
int local_row=0; local_row<num_rows; ++local_row) {
313 const int global_row = spatial_map->GID(local_row);
314 f_mp[nqp*local_row+qp] = (*f)[local_row];
315 A->ExtractMyRowView(local_row, num_entries, values, indices);
316 for (
int j=0;
j<num_entries; ++
j) {
317 int local_col = indices[
j];
318 int global_col = A->GCID(local_col);
319 block_indices[
j] = nqp*global_col+qp;
322 nqp*global_row+qp, num_entries, values, &block_indices[0]);
332 ParameterList precParams;
333 precParams.set(
"default values",
"SA");
334 precParams.set(
"ML output", 10);
335 precParams.set(
"max levels",5);
336 precParams.set(
"increasing or decreasing",
"increasing");
337 precParams.set(
"aggregation: type",
"Uncoupled");
338 precParams.set(
"smoother: type",
"ML symmetric Gauss-Seidel");
339 precParams.set(
"smoother: sweeps",2);
340 precParams.set(
"smoother: pre or post",
"both");
341 precParams.set(
"coarse: max size", 200);
342 precParams.set(
"PDE equations",nqp);
343 #ifdef HAVE_ML_AMESOS 344 precParams.set(
"coarse: type",
"Amesos-KLU");
346 precParams.set(
"coarse: type",
"Jacobi");
348 RCP<ML_Epetra::MultiLevelPreconditioner> M_mp;
350 TEUCHOS_FUNC_TIME_MONITOR(
"Preconditioner Setup");
351 M_mp = rcp(
new ML_Epetra::MultiLevelPreconditioner(A_mp, precParams));
356 if (krylov_method ==
GMRES)
357 aztec.SetAztecOption(AZ_solver, AZ_gmres);
358 else if (krylov_method ==
CG)
359 aztec.SetAztecOption(AZ_solver, AZ_cg);
360 aztec.SetAztecOption(AZ_precond, AZ_none);
361 aztec.SetAztecOption(AZ_kspace, 100);
362 aztec.SetAztecOption(AZ_conv, AZ_r0);
363 aztec.SetAztecOption(AZ_output, 10);
364 aztec.SetUserOperator(&A_mp);
365 aztec.SetPrecOperator(M_mp.get());
371 TEUCHOS_FUNC_TIME_MONITOR(
"System Solve");
372 aztec.Iterate(1000, tol);
380 for (
int qp=0; qp<nqp; qp++) {
381 TEUCHOS_FUNC_TIME_MONITOR(
"Compute Responses");
384 for (
int i=0; i<num_KL; i++)
385 (*p)[i] = quad_points[qp][i];
388 const int num_rows = spatial_map->NumMyElements();
389 for (
int local_row=0; local_row<num_rows; ++local_row)
390 (*
x)[local_row] = x_mp[nqp*local_row+qp];
403 g2.Multiply(1.0, *
g, *
g, 0.0);
404 g_mean.
Update(quad_weights[qp], *
g, 1.0);
405 g_var.
Update(quad_weights[qp], g2, 1.0);
407 g2.Multiply(1.0, g_mean, g_mean, 0.0);
408 g_var.
Update(-1.0, g2, 1.0);
411 for (
int i=0; i<g_var.
MyLength(); i++)
414 std::cout.precision(16);
415 std::cout <<
"\nResponse Mean = " << std::endl << g_mean << std::endl;
416 std::cout <<
"Response Std. Dev. = " << std::endl << g_var << std::endl;
420 TimeMonitor::summarize(std::cout);
421 TimeMonitor::zeroOutTimers();
425 catch (std::exception& e) {
426 std::cout << e.what() << std::endl;
428 catch (std::string& s) {
429 std::cout << s << std::endl;
432 std::cout << s << std::endl;
435 std::cout <<
"Caught unknown exception!" <<std:: endl;
KOKKOS_INLINE_FUNCTION PCE< Storage > sqrt(const PCE< Storage > &a)
Hermite polynomial basis.
const Krylov_Method krylov_method_values[]
const int num_krylov_method
Teuchos::RCP< const Epetra_Map > get_f_map() const
Return residual vector map.
const char * sg_quad_names[]
int InsertGlobalIndices(int_type GlobalRow, int NumIndices, int_type *Indices)
int FillComplete(bool OptimizeDataStorage=true)
const char * sg_rf_names[]
const char * krylov_method_names[]
int Scale(double ScalarValue)
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.
int main(int argc, char *argv[])
Teuchos::RCP< const Epetra_Map > get_x_map() const
Return solution vector map.
const char * sg_growth_names[]
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
const SG_RF sg_rf_values[]
const SG_GROWTH sg_growth_values[]
Legendre polynomial basis.
const SG_Quad sg_quad_values[]
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
void evalModel(const InArgs &inArgs, const OutArgs &outArgs) const
Evaluate model on InArgs.
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.
Teuchos::RCP< const Epetra_Map > get_p_map(int l) const
Return parameter vector map.
OutArgs createOutArgs() const
Create OutArgs.
virtual int ReplaceGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
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)