Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
linear2d_diffusion_pce_nox.cpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Stokhos Package
5// Copyright (2009) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
38//
39// ***********************************************************************
40// @HEADER
41
42// ModelEvaluator implementing our problem
43#include "twoD_diffusion_ME.hpp"
44
45// Epetra communicator
46#ifdef HAVE_MPI
47#include "Epetra_MpiComm.h"
48#else
49#include "Epetra_SerialComm.h"
50#endif
51
52// NOX
53#include "NOX.H"
54#include "NOX_Epetra.H"
55
56// Stokhos Stochastic Galerkin
57#include "Stokhos_Epetra.hpp"
58
59// Timing utilities
60#include "Teuchos_TimeMonitor.hpp"
61
62// I/O utilities
63#include "EpetraExt_VectorOut.h"
64
65int main(int argc, char *argv[]) {
66 int n = 32; // spatial discretization (per dimension)
67 int num_KL = 2; // number of KL terms
68 int p = 3; // polynomial order
69 double mu = 0.1; // mean of exponential random field
70 double s = 0.2; // std. dev. of exponential r.f.
71 bool nonlinear_expansion = false; // nonlinear expansion of diffusion coeff
72 // (e.g., log-normal)
73 bool matrix_free = true; // use matrix-free stochastic operator
74 bool symmetric = false; // use symmetric formulation
75
76 double g_mean_exp = 0.172988; // expected response mean
77 double g_std_dev_exp = 0.0380007; // expected response std. dev.
78 double g_tol = 1e-6; // tolerance on determining success
79
80// Initialize MPI
81#ifdef HAVE_MPI
82 MPI_Init(&argc,&argv);
83#endif
84
85 int MyPID;
86
87 try {
88
89 {
90 TEUCHOS_FUNC_TIME_MONITOR("Total PCE Calculation Time");
91
92 // Create a communicator for Epetra objects
93 Teuchos::RCP<const Epetra_Comm> globalComm;
94#ifdef HAVE_MPI
95 globalComm = Teuchos::rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
96#else
97 globalComm = Teuchos::rcp(new Epetra_SerialComm);
98#endif
99 MyPID = globalComm->MyPID();
100
101 // Create Stochastic Galerkin basis and expansion
102 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(num_KL);
103 for (int i=0; i<num_KL; i++)
104 bases[i] = Teuchos::rcp(new Stokhos::LegendreBasis<int,double>(p, true));
105 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
106 Teuchos::rcp(new Stokhos::CompletePolynomialBasis<int,double>(bases));
107 int sz = basis->size();
108 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk;
109 if (nonlinear_expansion)
110 Cijk = basis->computeTripleProductTensor();
111 else
112 Cijk = basis->computeLinearTripleProductTensor();
113 Teuchos::RCP<Stokhos::OrthogPolyExpansion<int,double> > expansion =
115 Cijk));
116 if (MyPID == 0)
117 std::cout << "Stochastic Galerkin expansion size = " << sz << std::endl;
118
119 // Create stochastic parallel distribution
120 int num_spatial_procs = -1;
121 if (argc > 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 =
126 Teuchos::rcp(new Stokhos::ParallelData(basis, Cijk, globalComm,
127 parallelParams));
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();
132
133 // Create application
134 Teuchos::RCP<twoD_diffusion_ME> model =
135 Teuchos::rcp(new twoD_diffusion_ME(app_comm, n, num_KL, mu, s, basis,
136 nonlinear_expansion, symmetric));
137
138 // Setup stochastic Galerkin algorithmic parameters
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");
148 }
149 if (matrix_free) {
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");
167#else
168 precParams.set("coarse: type","Jacobi");
169#endif
170 }
171 else {
172 sgOpParams.set("Operator Method", "Fully Assembled");
173 sgPrecParams.set("Preconditioner Method", "None");
174 }
175
176 // Create stochastic Galerkin model evaluator
177 Teuchos::RCP<Stokhos::SGModelEvaluator> sg_model =
178 Teuchos::rcp(new Stokhos::SGModelEvaluator(model, basis, Teuchos::null,
179 expansion, sg_parallel_data,
180 sgParams));
181
182 // Set up stochastic parameters
183 // The current implementation of the model doesn't actually use these
184 // values, but is hard-coded to certain uncertainty models
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];
193 }
194 sg_model->set_p_sg_init(0, *sg_p_init);
195
196 // Setup stochastic initial guess
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);
201
202 // Set up NOX parameters
203 Teuchos::RCP<Teuchos::ParameterList> noxParams =
204 Teuchos::rcp(new Teuchos::ParameterList);
205
206 // Set the nonlinear solver method
207 noxParams->set("Nonlinear Solver", "Line Search Based");
208
209 // Set the printing parameters in the "Printing" sublist
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 +
220 NOX::Utils::Error);
221
222 // Create printing utilities
223 NOX::Utils utils(printParams);
224
225 // Sublist for line search
226 Teuchos::ParameterList& searchParams = noxParams->sublist("Line Search");
227 searchParams.set("Method", "Full Step");
228
229 // Sublist for direction
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");
234
235 // Sublist for linear solver for the Newton method
236 Teuchos::ParameterList& lsParams = newtonParams.sublist("Linear Solver");
237 if (symmetric)
238 lsParams.set("Aztec Solver", "CG");
239 else
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);
245 if (matrix_free)
246 lsParams.set("Preconditioner", "User Defined");
247 else {
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);
253 }
254
255 // Sublist for convergence tests
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);
267
268 // Create NOX interface
269 Teuchos::RCP<NOX::Epetra::ModelEvaluatorInterface> nox_interface =
270 Teuchos::rcp(new NOX::Epetra::ModelEvaluatorInterface(sg_model));
271
272 // Create NOX linear system object
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;
278 if (matrix_free) {
279 Teuchos::RCP<Epetra_Operator> M = sg_model->create_WPrec()->PrecOp;
280 Teuchos::RCP<NOX::Epetra::Interface::Preconditioner> iPrec = nox_interface;
281 linsys =
282 Teuchos::rcp(new NOX::Epetra::LinearSystemAztecOO(printParams, lsParams,
283 iJac, A, iPrec, M,
284 *u));
285 }
286 else {
287 linsys =
288 Teuchos::rcp(new NOX::Epetra::LinearSystemAztecOO(printParams, lsParams,
289 iReq, iJac, A,
290 *u));
291 }
292
293 // Build NOX group
294 Teuchos::RCP<NOX::Epetra::Group> grp =
295 Teuchos::rcp(new NOX::Epetra::Group(printParams, iReq, *u, linsys));
296
297 // Create the Solver convergence test
298 Teuchos::RCP<NOX::StatusTest::Generic> statusTests =
299 NOX::StatusTest::buildStatusTests(statusParams, utils);
300
301 // Create the solver
302 Teuchos::RCP<NOX::Solver::Generic> solver =
303 NOX::Solver::buildSolver(grp, statusTests, noxParams);
304
305 // Solve the system
306 NOX::StatusTest::StatusType status = solver->solve();
307
308 // Get final solution
309 const NOX::Epetra::Group& finalGroup =
310 dynamic_cast<const NOX::Epetra::Group&>(solver->getSolutionGroup());
311 const Epetra_Vector& finalSolution =
312 (dynamic_cast<const NOX::Epetra::Vector&>(finalGroup.getX())).getEpetraVector();
313
314 // Save final solution to file
315 EpetraExt::VectorToMatrixMarketFile("nox_stochastic_solution.mm",
316 finalSolution);
317
318 // Save mean and variance to file
319 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_x_poly =
320 sg_model->create_x_sg(View, &finalSolution);
321 Epetra_Vector mean(*(model->get_x_map()));
322 Epetra_Vector std_dev(*(model->get_x_map()));
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);
327
328 // Evaluate SG responses at SG parameters
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 =
334 Teuchos::rcp(new Epetra_Vector(*(sg_model->get_g_map(0))));
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);
339
340 // Print mean and standard deviation of response
341 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_g_poly =
342 sg_model->create_g_sg(0, View, sg_g.get());
343 Epetra_Vector g_mean(*(model->get_g_map(0)));
344 Epetra_Vector g_std_dev(*(model->get_g_map(0)));
345 sg_g_poly->computeMean(g_mean);
346 sg_g_poly->computeStandardDeviation(g_std_dev);
347 std::cout.precision(16);
348 // std::cout << "\nResponse Expansion = " << std::endl;
349 // std::cout.precision(12);
350 // sg_g_poly->print(std::cout);
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;
354
355 // Determine if example passed
356 bool passed = false;
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)
360 passed = true;
361 if (MyPID == 0) {
362 if (passed)
363 std::cout << "Example Passed!" << std::endl;
364 else
365 std::cout << "Example Failed!" << std::endl;
366 }
367
368 }
369
370 Teuchos::TimeMonitor::summarize(std::cout);
371 Teuchos::TimeMonitor::zeroOutTimers();
372
373 }
374
375 catch (std::exception& e) {
376 std::cout << e.what() << std::endl;
377 }
378 catch (std::string& s) {
379 std::cout << s << std::endl;
380 }
381 catch (char *s) {
382 std::cout << s << std::endl;
383 }
384 catch (...) {
385 std::cout << "Caught unknown exception!" << std::endl;
386 }
387
388#ifdef HAVE_MPI
389 MPI_Finalize() ;
390#endif
391
392}
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[])