47 #include "Epetra_MpiComm.h" 49 #include "Epetra_SerialComm.h" 60 #include "Teuchos_TimeMonitor.hpp" 63 #include "EpetraExt_VectorOut.h" 64 #include "EpetraExt_RowMatrixOut.h" 72 bool nonlinear_expansion =
false;
74 bool symmetric =
false;
76 double g_mean_exp = 0.172988;
77 double g_std_dev_exp = 0.0380007;
82 MPI_Init(&argc,&
argv);
90 TEUCHOS_FUNC_TIME_MONITOR(
"Total PCE Calculation Time");
93 Teuchos::RCP<const Epetra_Comm> globalComm;
95 globalComm = Teuchos::rcp(
new Epetra_MpiComm(MPI_COMM_WORLD));
97 globalComm = Teuchos::rcp(
new Epetra_SerialComm);
99 MyPID = globalComm->MyPID();
102 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(num_KL);
103 for (
int i=0; i<num_KL; i++)
105 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
108 int sz = basis->size();
109 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk;
110 if (nonlinear_expansion)
111 Cijk = basis->computeTripleProductTensor();
113 Cijk = basis->computeLinearTripleProductTensor();
114 Teuchos::RCP<Stokhos::OrthogPolyExpansion<int,double> > expansion =
118 std::cout <<
"Stochastic Galerkin expansion size = " << sz << std::endl;
121 int num_spatial_procs = -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 if (!nonlinear_expansion) {
145 sgParams->set(
"Parameter Expansion Type",
"Linear");
146 sgParams->set(
"Jacobian Expansion Type",
"Linear");
149 Teuchos::ParameterList precParams;
150 precParams.set(
"default values",
"SA");
151 precParams.set(
"ML output", 0);
152 precParams.set(
"max levels",5);
153 precParams.set(
"increasing or decreasing",
"increasing");
154 precParams.set(
"aggregation: type",
"Uncoupled");
155 precParams.set(
"smoother: type",
"ML symmetric Gauss-Seidel");
156 precParams.set(
"smoother: sweeps",2);
157 precParams.set(
"smoother: pre or post",
"both");
158 precParams.set(
"coarse: max size", 200);
159 precParams.set(
"PDE equations",sz);
160 #ifdef HAVE_ML_AMESOS 161 precParams.set(
"coarse: type",
"Amesos-KLU");
163 precParams.set(
"coarse: type",
"Jacobi");
167 Teuchos::RCP<Stokhos::SGModelEvaluator_Interlaced> sg_model =
169 model, basis, Teuchos::null,
170 expansion, sg_parallel_data,
176 Teuchos::Array<double> point(num_KL, 1.0);
177 Teuchos::Array<double> basis_vals(sz);
178 basis->evaluateBases(point, basis_vals);
179 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_p_poly =
180 sg_model->create_p_sg(0);
181 for (
int i=0; i<num_KL; i++) {
182 sg_p_poly->term(i,0)[i] = 0.0;
183 sg_p_poly->term(i,1)[i] = 1.0 / basis_vals[i+1];
187 Teuchos::RCP<const Epetra_Vector> sg_p = sg_p_poly->getBlockVector();
188 Teuchos::RCP<Epetra_Vector> sg_x =
189 Teuchos::rcp(
new Epetra_Vector(*(sg_model->get_x_map())));
190 sg_x->PutScalar(0.0);
191 Teuchos::RCP<Epetra_Vector> sg_f =
192 Teuchos::rcp(
new Epetra_Vector(*(sg_model->get_f_map())));
193 Teuchos::RCP<Epetra_Vector> sg_dx =
194 Teuchos::rcp(
new Epetra_Vector(*(sg_model->get_x_map())));
195 Teuchos::RCP<Epetra_CrsMatrix> sg_J =
196 Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(sg_model->create_W());
197 Teuchos::RCP<ML_Epetra::MultiLevelPreconditioner> sg_M =
198 Teuchos::rcp(
new ML_Epetra::MultiLevelPreconditioner(*sg_J, precParams,
202 EpetraExt::ModelEvaluator::InArgs sg_inArgs = sg_model->createInArgs();
203 EpetraExt::ModelEvaluator::OutArgs sg_outArgs = sg_model->createOutArgs();
204 sg_inArgs.set_p(1, sg_p);
205 sg_inArgs.set_x(sg_x);
206 sg_outArgs.set_f(sg_f);
207 sg_outArgs.set_W(sg_J);
210 sg_model->evalModel(sg_inArgs, sg_outArgs);
211 sg_M->ComputePreconditioner();
215 sg_f->Norm2(&norm_f);
217 std::cout <<
"\nInitial residual norm = " << norm_f << std::endl;
222 aztec.SetAztecOption(AZ_solver, AZ_cg);
224 aztec.SetAztecOption(AZ_solver, AZ_gmres);
225 aztec.SetAztecOption(AZ_precond, AZ_none);
226 aztec.SetAztecOption(AZ_kspace, 20);
227 aztec.SetAztecOption(AZ_conv, AZ_r0);
228 aztec.SetAztecOption(AZ_output, 1);
229 aztec.SetUserOperator(sg_J.get());
230 aztec.SetPrecOperator(sg_M.get());
231 aztec.SetLHS(sg_dx.get());
232 aztec.SetRHS(sg_f.get());
235 aztec.Iterate(1000, 1e-12);
238 sg_x->Update(-1.0, *sg_dx, 1.0);
241 EpetraExt::VectorToMatrixMarketFile(
"stochastic_solution_interlaced.mm",
245 EpetraExt::VectorToMatrixMarketFile(
"stochastic_RHS_interlaced.mm",
249 EpetraExt::RowMatrixToMatrixMarketFile(
"stochastic_operator_interlaced.mm",
253 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_x_poly =
254 sg_model->create_x_sg(View, sg_x.get());
255 Epetra_Vector mean(*(model->get_x_map()));
256 Epetra_Vector std_dev(*(model->get_x_map()));
257 sg_x_poly->computeMean(mean);
258 sg_x_poly->computeStandardDeviation(std_dev);
259 EpetraExt::VectorToMatrixMarketFile(
"mean_gal_interlaced.mm", mean);
260 EpetraExt::VectorToMatrixMarketFile(
"std_dev_gal_interlaced.mm", std_dev);
263 EpetraExt::ModelEvaluator::OutArgs sg_outArgs2 = sg_model->createOutArgs();
264 Teuchos::RCP<Epetra_Vector> sg_g =
265 Teuchos::rcp(
new Epetra_Vector(*(sg_model->get_g_map(0))));
266 sg_f->PutScalar(0.0);
267 sg_outArgs2.set_f(sg_f);
268 sg_outArgs2.set_g(0, sg_g);
269 sg_model->evalModel(sg_inArgs, sg_outArgs2);
272 sg_f->Norm2(&norm_f);
274 std::cout <<
"\nFinal residual norm = " << norm_f << std::endl;
277 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_g_poly =
278 sg_model->create_g_sg(0, View, sg_g.get());
279 Epetra_Vector g_mean(*(model->get_g_map(0)));
280 Epetra_Vector g_std_dev(*(model->get_g_map(0)));
281 sg_g_poly->computeMean(g_mean);
282 sg_g_poly->computeStandardDeviation(g_std_dev);
283 std::cout.precision(16);
287 std::cout << std::endl;
288 std::cout <<
"Response Mean = " << std::endl << g_mean << std::endl;
289 std::cout <<
"Response Std. Dev. = " << std::endl << g_std_dev << std::endl;
293 if (norm_f < 1.0e-10 &&
294 std::abs(g_mean[0]-g_mean_exp) < g_tol &&
295 std::abs(g_std_dev[0]-g_std_dev_exp) < g_tol)
299 std::cout <<
"Example Passed!" << std::endl;
301 std::cout <<
"Example Failed!" << std::endl;
306 Teuchos::TimeMonitor::summarize(std::cout);
307 Teuchos::TimeMonitor::zeroOutTimers();
311 catch (std::exception& e) {
312 std::cout << e.what() << std::endl;
314 catch (std::string& s) {
315 std::cout << s << std::endl;
318 std::cout << s << std::endl;
321 std::cout <<
"Caught unknown exception!" << std::endl;
Orthogonal polynomial expansions limited to algebraic operations.
int main(int argc, char *argv[])
ModelEvaluator for a linear 2-D diffusion problem.
KOKKOS_INLINE_FUNCTION PCE< Storage > abs(const PCE< Storage > &a)
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Legendre polynomial basis.
Nonlinear, stochastic Galerkin ModelEvaluator that constructs a interlaced Jacobian.