59#include "Teuchos_TimeMonitor.hpp"
62#include "EpetraExt_VectorOut.h"
70 bool nonlinear_expansion =
false;
72 bool symmetric =
false;
74 double g_mean_exp = 0.172988;
75 double g_std_dev_exp = 0.0380007;
80 MPI_Init(&argc,&
argv);
88 TEUCHOS_FUNC_TIME_MONITOR(
"Total PCE Calculation Time");
91 Teuchos::RCP<const Epetra_Comm> globalComm;
97 MyPID = globalComm->MyPID();
100 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(num_KL);
101 for (
int i=0; i<num_KL; i++)
103 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
106 int sz = basis->size();
107 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk;
108 if (nonlinear_expansion)
109 Cijk = basis->computeTripleProductTensor();
111 Cijk = basis->computeLinearTripleProductTensor();
112 Teuchos::RCP<Stokhos::OrthogPolyExpansion<int,double> > expansion =
116 std::cout <<
"Stochastic Galerkin expansion size = " << sz << std::endl;
119 int num_spatial_procs = -1;
121 num_spatial_procs = std::atoi(
argv[1]);
122 Teuchos::ParameterList parallelParams;
123 parallelParams.set(
"Number of Spatial Processors", num_spatial_procs);
128 Teuchos::RCP<Stokhos::ParallelData> sg_parallel_data =
131 Teuchos::RCP<const EpetraExt::MultiComm> sg_comm =
132 sg_parallel_data->getMultiComm();
133 Teuchos::RCP<const Epetra_Comm> app_comm =
134 sg_parallel_data->getSpatialComm();
137 Teuchos::RCP<twoD_diffusion_ME> model =
139 nonlinear_expansion, symmetric));
142 Teuchos::RCP<Teuchos::ParameterList> sgParams =
143 Teuchos::rcp(
new Teuchos::ParameterList);
144 Teuchos::ParameterList& sgOpParams =
145 sgParams->sublist(
"SG Operator");
146 Teuchos::ParameterList& sgPrecParams =
147 sgParams->sublist(
"SG Preconditioner");
148 if (!nonlinear_expansion) {
149 sgParams->set(
"Parameter Expansion Type",
"Linear");
150 sgParams->set(
"Jacobian Expansion Type",
"Linear");
152 sgOpParams.set(
"Operator Method",
"Matrix Free");
153 sgPrecParams.set(
"Preconditioner Method",
"Mean-based");
156 sgPrecParams.set(
"Symmetric Gauss-Seidel", symmetric);
157 sgPrecParams.set(
"Mean Preconditioner Type",
"ML");
158 Teuchos::ParameterList& precParams =
159 sgPrecParams.sublist(
"Mean Preconditioner Parameters");
160 precParams.set(
"default values",
"SA");
161 precParams.set(
"ML output", 0);
162 precParams.set(
"max levels",5);
163 precParams.set(
"increasing or decreasing",
"increasing");
164 precParams.set(
"aggregation: type",
"Uncoupled");
165 precParams.set(
"smoother: type",
"ML symmetric Gauss-Seidel");
166 precParams.set(
"smoother: sweeps",2);
167 precParams.set(
"smoother: pre or post",
"both");
168 precParams.set(
"coarse: max size", 200);
170 precParams.set(
"coarse: type",
"Amesos-KLU");
172 precParams.set(
"coarse: type",
"Jacobi");
176 Teuchos::RCP<Stokhos::SGModelEvaluator> sg_model =
178 expansion, sg_parallel_data,
184 Teuchos::Array<double> point(num_KL, 1.0);
185 Teuchos::Array<double> basis_vals(sz);
186 basis->evaluateBases(point, basis_vals);
187 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_p_poly =
188 sg_model->create_p_sg(0);
189 for (
int i=0; i<num_KL; i++) {
190 sg_p_poly->term(i,0)[i] = 0.0;
191 sg_p_poly->term(i,1)[i] = 1.0 / basis_vals[i+1];
193 std::cout << *sg_p_poly << std::endl;
196 Teuchos::RCP<const Epetra_Vector> sg_p = sg_p_poly->getBlockVector();
197 Teuchos::RCP<Epetra_Vector> sg_x =
199 sg_x->PutScalar(0.0);
200 Teuchos::RCP<Epetra_Vector> sg_f =
202 Teuchos::RCP<Epetra_Vector> sg_dx =
204 Teuchos::RCP<Epetra_Operator> sg_J = sg_model->create_W();
205 Teuchos::RCP<Epetra_Operator> sg_M = sg_model->create_WPrec()->PrecOp;
209 EpetraExt::ModelEvaluator::InArgs sg_inArgs = sg_model->createInArgs();
210 EpetraExt::ModelEvaluator::OutArgs sg_outArgs = sg_model->createOutArgs();
211 sg_inArgs.set_p(1, sg_p);
212 sg_inArgs.set_x(sg_x);
213 sg_outArgs.set_f(sg_f);
214 sg_outArgs.set_W(sg_J);
215 sg_outArgs.set_WPrec(sg_M);
218 sg_model->evalModel(sg_inArgs, sg_outArgs);
222 sg_f->Norm2(&norm_f);
224 std::cout <<
"\nInitial residual norm = " << norm_f << std::endl;
229 aztec.SetAztecOption(AZ_solver, AZ_cg);
231 aztec.SetAztecOption(AZ_solver, AZ_gmres);
232 aztec.SetAztecOption(AZ_precond, AZ_none);
233 aztec.SetAztecOption(AZ_kspace, 20);
234 aztec.SetAztecOption(AZ_conv, AZ_r0);
235 aztec.SetAztecOption(AZ_output, 1);
236 aztec.SetUserOperator(sg_J.get());
237 aztec.SetPrecOperator(sg_M.get());
238 aztec.SetLHS(sg_dx.get());
239 aztec.SetRHS(sg_f.get());
242 aztec.Iterate(1000, 1e-12);
245 sg_x->Update(-1.0, *sg_dx, 1.0);
248 EpetraExt::VectorToMatrixMarketFile(
"stochastic_solution.mm", *sg_x);
251 EpetraExt::VectorToMatrixMarketFile(
"stochastic_RHS.mm", *sg_f);
254 Teuchos::RCP<Epetra_CrsMatrix> sg_J_crs =
255 Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(sg_J);
256 if (sg_J_crs != Teuchos::null)
257 EpetraExt::RowMatrixToMatrixMarketFile(
"stochastic_operator.mm",
261 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_x_poly =
262 sg_model->create_x_sg(
View, sg_x.get());
265 sg_x_poly->computeMean(mean);
266 sg_x_poly->computeStandardDeviation(std_dev);
267 EpetraExt::VectorToMatrixMarketFile(
"mean_gal.mm", mean);
268 EpetraExt::VectorToMatrixMarketFile(
"std_dev_gal.mm", std_dev);
271 EpetraExt::ModelEvaluator::OutArgs sg_outArgs2 = sg_model->createOutArgs();
272 Teuchos::RCP<Epetra_Vector> sg_g =
274 sg_f->PutScalar(0.0);
275 sg_outArgs2.set_f(sg_f);
276 sg_outArgs2.set_g(0, sg_g);
277 sg_model->evalModel(sg_inArgs, sg_outArgs2);
280 sg_f->Norm2(&norm_f);
282 std::cout <<
"\nFinal residual norm = " << norm_f << std::endl;
285 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_g_poly =
286 sg_model->create_g_sg(0,
View, sg_g.get());
289 sg_g_poly->computeMean(g_mean);
290 sg_g_poly->computeStandardDeviation(g_std_dev);
291 std::cout.precision(16);
295 std::cout << std::endl;
296 std::cout <<
"Response Mean = " << std::endl << g_mean << std::endl;
297 std::cout <<
"Response Std. Dev. = " << std::endl << g_std_dev << std::endl;
301 if (norm_f < 1.0e-10 &&
302 std::abs(g_mean[0]-g_mean_exp) < g_tol &&
303 std::abs(g_std_dev[0]-g_std_dev_exp) < g_tol)
307 std::cout <<
"Example Passed!" << std::endl;
309 std::cout <<
"Example Failed!" << std::endl;
314 Teuchos::TimeMonitor::summarize(std::cout);
315 Teuchos::TimeMonitor::zeroOutTimers();
319 catch (std::exception& e) {
320 std::cout << e.what() << std::endl;
322 catch (std::string& s) {
323 std::cout << s << std::endl;
326 std::cout << s << std::endl;
329 std::cout <<
"Caught unknown exception!" << std::endl;
Orthogonal polynomial expansions limited to algebraic operations.
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Legendre polynomial basis.
Nonlinear, stochastic Galerkin ModelEvaluator.
ModelEvaluator for a linear 2-D diffusion problem.
int main(int argc, char *argv[])