54 #include "NOX_Epetra.H" 60 #include "Teuchos_TimeMonitor.hpp" 63 #include "EpetraExt_VectorOut.h" 71 bool nonlinear_expansion =
false;
73 bool matrix_free =
true;
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;
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 =
107 int sz = basis->size();
108 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk;
109 if (nonlinear_expansion)
110 Cijk = basis->computeTripleProductTensor();
112 Cijk = basis->computeLinearTripleProductTensor();
113 Teuchos::RCP<Stokhos::OrthogPolyExpansion<int,double> > expansion =
117 std::cout <<
"Stochastic Galerkin expansion size = " << sz << std::endl;
120 int num_spatial_procs = -1;
122 num_spatial_procs = std::atoi(
argv[1]);
123 Teuchos::ParameterList parallelParams;
124 parallelParams.set(
"Number of Spatial Processors", num_spatial_procs);
125 Teuchos::RCP<Stokhos::ParallelData> sg_parallel_data =
128 Teuchos::RCP<const EpetraExt::MultiComm> sg_comm =
129 sg_parallel_data->getMultiComm();
130 Teuchos::RCP<const Epetra_Comm> app_comm =
131 sg_parallel_data->getSpatialComm();
134 Teuchos::RCP<twoD_diffusion_ME> model =
136 nonlinear_expansion, symmetric));
139 Teuchos::RCP<Teuchos::ParameterList> sgParams =
140 Teuchos::rcp(
new Teuchos::ParameterList);
141 Teuchos::ParameterList& sgOpParams =
142 sgParams->sublist(
"SG Operator");
143 Teuchos::ParameterList& sgPrecParams =
144 sgParams->sublist(
"SG Preconditioner");
145 if (!nonlinear_expansion) {
146 sgParams->set(
"Parameter Expansion Type",
"Linear");
147 sgParams->set(
"Jacobian Expansion Type",
"Linear");
150 sgOpParams.set(
"Operator Method",
"Matrix Free");
151 sgPrecParams.set(
"Preconditioner Method",
"Approximate Gauss-Seidel");
152 sgPrecParams.set(
"Symmetric Gauss-Seidel", symmetric);
153 sgPrecParams.set(
"Mean Preconditioner Type",
"ML");
154 Teuchos::ParameterList& precParams =
155 sgPrecParams.sublist(
"Mean Preconditioner Parameters");
156 precParams.set(
"default values",
"SA");
157 precParams.set(
"ML output", 0);
158 precParams.set(
"max levels",5);
159 precParams.set(
"increasing or decreasing",
"increasing");
160 precParams.set(
"aggregation: type",
"Uncoupled");
161 precParams.set(
"smoother: type",
"ML symmetric Gauss-Seidel");
162 precParams.set(
"smoother: sweeps",2);
163 precParams.set(
"smoother: pre or post",
"both");
164 precParams.set(
"coarse: max size", 200);
165 #ifdef HAVE_ML_AMESOS 166 precParams.set(
"coarse: type",
"Amesos-KLU");
168 precParams.set(
"coarse: type",
"Jacobi");
172 sgOpParams.set(
"Operator Method",
"Fully Assembled");
173 sgPrecParams.set(
"Preconditioner Method",
"None");
177 Teuchos::RCP<Stokhos::SGModelEvaluator> sg_model =
179 expansion, sg_parallel_data,
185 Teuchos::Array<double> point(num_KL, 1.0);
186 Teuchos::Array<double> basis_vals(sz);
187 basis->evaluateBases(point, basis_vals);
188 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_p_init =
189 sg_model->create_p_sg(0);
190 for (
int i=0; i<num_KL; i++) {
191 sg_p_init->term(i,0)[i] = 0.0;
192 sg_p_init->term(i,1)[i] = 1.0 / basis_vals[i+1];
194 sg_model->set_p_sg_init(0, *sg_p_init);
197 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_x_init =
198 sg_model->create_x_sg();
199 sg_x_init->init(0.0);
200 sg_model->set_x_sg_init(*sg_x_init);
203 Teuchos::RCP<Teuchos::ParameterList> noxParams =
204 Teuchos::rcp(
new Teuchos::ParameterList);
207 noxParams->set(
"Nonlinear Solver",
"Line Search Based");
210 Teuchos::ParameterList& printParams = noxParams->sublist(
"Printing");
211 printParams.set(
"MyPID", MyPID);
212 printParams.set(
"Output Precision", 3);
213 printParams.set(
"Output Processor", 0);
214 printParams.set(
"Output Information",
215 NOX::Utils::OuterIteration +
216 NOX::Utils::OuterIterationStatusTest +
217 NOX::Utils::InnerIteration +
218 NOX::Utils::LinearSolverDetails +
219 NOX::Utils::Warning +
223 NOX::Utils utils(printParams);
226 Teuchos::ParameterList& searchParams = noxParams->sublist(
"Line Search");
227 searchParams.set(
"Method",
"Full Step");
230 Teuchos::ParameterList& dirParams = noxParams->sublist(
"Direction");
231 dirParams.set(
"Method",
"Newton");
232 Teuchos::ParameterList& newtonParams = dirParams.sublist(
"Newton");
233 newtonParams.set(
"Forcing Term Method",
"Constant");
236 Teuchos::ParameterList& lsParams = newtonParams.sublist(
"Linear Solver");
238 lsParams.set(
"Aztec Solver",
"CG");
240 lsParams.set(
"Aztec Solver",
"GMRES");
241 lsParams.set(
"Max Iterations", 1000);
242 lsParams.set(
"Size of Krylov Subspace", 100);
243 lsParams.set(
"Tolerance", 1e-12);
244 lsParams.set(
"Output Frequency", 1);
246 lsParams.set(
"Preconditioner",
"User Defined");
248 lsParams.set(
"Preconditioner",
"ML");
249 Teuchos::ParameterList& precParams =
250 lsParams.sublist(
"ML");
251 ML_Epetra::SetDefaults(
"DD", precParams);
252 lsParams.set(
"Write Linear System",
false);
256 Teuchos::ParameterList& statusParams = noxParams->sublist(
"Status Tests");
257 statusParams.set(
"Test Type",
"Combo");
258 statusParams.set(
"Number of Tests", 2);
259 statusParams.set(
"Combo Type",
"OR");
260 Teuchos::ParameterList& normF = statusParams.sublist(
"Test 0");
261 normF.set(
"Test Type",
"NormF");
262 normF.set(
"Tolerance", 1e-10);
263 normF.set(
"Scale Type",
"Scaled");
264 Teuchos::ParameterList& maxIters = statusParams.sublist(
"Test 1");
265 maxIters.set(
"Test Type",
"MaxIters");
266 maxIters.set(
"Maximum Iterations", 1);
269 Teuchos::RCP<NOX::Epetra::ModelEvaluatorInterface> nox_interface =
270 Teuchos::rcp(
new NOX::Epetra::ModelEvaluatorInterface(sg_model));
273 Teuchos::RCP<const Epetra_Vector> u = sg_model->get_x_init();
274 Teuchos::RCP<Epetra_Operator> A = sg_model->create_W();
275 Teuchos::RCP<NOX::Epetra::Interface::Required> iReq = nox_interface;
276 Teuchos::RCP<NOX::Epetra::Interface::Jacobian> iJac = nox_interface;
277 Teuchos::RCP<NOX::Epetra::LinearSystemAztecOO> linsys;
279 Teuchos::RCP<Epetra_Operator> M = sg_model->create_WPrec()->PrecOp;
280 Teuchos::RCP<NOX::Epetra::Interface::Preconditioner> iPrec = nox_interface;
282 Teuchos::rcp(
new NOX::Epetra::LinearSystemAztecOO(printParams, lsParams,
288 Teuchos::rcp(
new NOX::Epetra::LinearSystemAztecOO(printParams, lsParams,
294 Teuchos::RCP<NOX::Epetra::Group> grp =
295 Teuchos::rcp(
new NOX::Epetra::Group(printParams, iReq, *u, linsys));
298 Teuchos::RCP<NOX::StatusTest::Generic> statusTests =
299 NOX::StatusTest::buildStatusTests(statusParams, utils);
302 Teuchos::RCP<NOX::Solver::Generic> solver =
303 NOX::Solver::buildSolver(grp, statusTests, noxParams);
306 NOX::StatusTest::StatusType status = solver->solve();
309 const NOX::Epetra::Group& finalGroup =
310 dynamic_cast<const NOX::Epetra::Group&
>(solver->getSolutionGroup());
312 (
dynamic_cast<const NOX::Epetra::Vector&
>(finalGroup.getX())).getEpetraVector();
315 EpetraExt::VectorToMatrixMarketFile(
"nox_stochastic_solution.mm",
319 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_x_poly =
320 sg_model->create_x_sg(
View, &finalSolution);
323 sg_x_poly->computeMean(mean);
324 sg_x_poly->computeStandardDeviation(std_dev);
325 EpetraExt::VectorToMatrixMarketFile(
"mean_gal.mm", mean);
326 EpetraExt::VectorToMatrixMarketFile(
"std_dev_gal.mm", std_dev);
329 EpetraExt::ModelEvaluator::InArgs sg_inArgs = sg_model->createInArgs();
330 EpetraExt::ModelEvaluator::OutArgs sg_outArgs =
331 sg_model->createOutArgs();
332 Teuchos::RCP<const Epetra_Vector> sg_p = sg_model->get_p_init(1);
333 Teuchos::RCP<Epetra_Vector> sg_g =
335 sg_inArgs.set_p(1, sg_p);
336 sg_inArgs.set_x(Teuchos::rcp(&finalSolution,
false));
337 sg_outArgs.set_g(0, sg_g);
338 sg_model->evalModel(sg_inArgs, sg_outArgs);
341 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_g_poly =
342 sg_model->create_g_sg(0,
View, sg_g.get());
345 sg_g_poly->computeMean(g_mean);
346 sg_g_poly->computeStandardDeviation(g_std_dev);
347 std::cout.precision(16);
351 std::cout << std::endl;
352 std::cout <<
"Response Mean = " << std::endl << g_mean << std::endl;
353 std::cout <<
"Response Std. Dev. = " << std::endl << g_std_dev << std::endl;
357 if (status == NOX::StatusTest::Converged &&
358 std::abs(g_mean[0]-g_mean_exp) < g_tol &&
359 std::abs(g_std_dev[0]-g_std_dev_exp) < g_tol)
363 std::cout <<
"Example Passed!" << std::endl;
365 std::cout <<
"Example Failed!" << std::endl;
370 Teuchos::TimeMonitor::summarize(std::cout);
371 Teuchos::TimeMonitor::zeroOutTimers();
375 catch (std::exception& e) {
376 std::cout << e.what() << std::endl;
378 catch (std::string& s) {
379 std::cout << s << std::endl;
382 std::cout << s << std::endl;
385 std::cout <<
"Caught unknown exception!" << std::endl;
Orthogonal polynomial expansions limited to algebraic operations.
Nonlinear, stochastic Galerkin ModelEvaluator.
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.
int main(int argc, char *argv[])