Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
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 
65 int 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.
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[])