Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
linear2d_diffusion_collocation.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// AztecOO solver
53#include "AztecOO.h"
54
55// Stokhos Stochastic Galerkin
56#include "Stokhos.hpp"
57
58// Timing utilities
59#include "Teuchos_CommandLineProcessor.hpp"
60#include "Teuchos_TimeMonitor.hpp"
61
62// I/O utilities
63#include "EpetraExt_VectorOut.h"
64
65// Krylov methods
67const int num_krylov_method = 2;
69const char *krylov_method_names[] = { "GMRES", "CG" };
70
71// Collocation preconditioning strategies
73const int num_prec_strategy = 3;
75const char *prec_strategy_names[] = { "Mean", "Reuse", "Rebuild" };
76
77// Random field types
79const int num_sg_rf = 4;
81const char *sg_rf_names[] = { "Uniform", "CC-Uniform", "Rys", "Log-Normal" };
82
83// Quadrature types
85const int num_sg_quad = 2;
87const char *sg_quad_names[] = { "tensor", "sparse-grid" };
88
89// Sparse grid growth rules
91const int num_sg_growth = 3;
94const char *sg_growth_names[] = {
95 "Slow Restricted", "Moderate Restricted", "Unrestricted" };
96
97int main(int argc, char *argv[]) {
98
99// Initialize MPI
100#ifdef HAVE_MPI
101 MPI_Init(&argc,&argv);
102#endif
103
104 // Create a communicator for Epetra objects
105 Teuchos::RCP<Epetra_Comm> Comm;
106#ifdef HAVE_MPI
107 Comm = Teuchos::rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
108#else
109 Comm = Teuchos::rcp(new Epetra_SerialComm);
110#endif
111
112 int MyPID = Comm->MyPID();
113
114 try {
115
116 // Setup command line options
117 Teuchos::CommandLineProcessor CLP;
118 CLP.setDocString(
119 "This example runs a stochastic collocation method.\n");
120
121 int n = 32;
122 CLP.setOption("num_mesh", &n, "Number of mesh points in each direction");
123
124 SG_RF randField = UNIFORM;
125 CLP.setOption("rand_field", &randField,
127 "Random field type");
128
129 double mean = 0.2;
130 CLP.setOption("mean", &mean, "Mean");
131
132 double sigma = 0.1;
133 CLP.setOption("std_dev", &sigma, "Standard deviation");
134
135 double weightCut = 1.0;
136 CLP.setOption("weight_cut", &weightCut, "Weight cut");
137
138 int num_KL = 2;
139 CLP.setOption("num_kl", &num_KL, "Number of KL terms");
140
141 int p = 3;
142 CLP.setOption("order", &p, "Polynomial order");
143
144 bool normalize_basis = true;
145 CLP.setOption("normalize", "unnormalize", &normalize_basis,
146 "Normalize PC basis");
147
148 Krylov_Method krylov_method = CG;
149 CLP.setOption("krylov_method", &krylov_method,
151 "Krylov method");
152
153 PrecStrategy precStrategy = MEAN;
154 CLP.setOption("prec_strategy", &precStrategy,
156 "Preconditioner strategy");
157
158 double tol = 1e-12;
159 CLP.setOption("tol", &tol, "Solver tolerance");
160
161#ifdef HAVE_STOKHOS_DAKOTA
162 SG_Quad quad_method = SPARSE_GRID;
163#else
164 SG_Quad quad_method = TENSOR;
165#endif
166 CLP.setOption("quadrature", &quad_method,
168 "Quadrature type");
169
170 SG_GROWTH sg_growth = MODERATE_RESTRICTED;
171 CLP.setOption("sg_growth", &sg_growth,
173 "Sparse grid growth rule");
174
175 CLP.parse( argc, argv );
176
177 if (MyPID == 0) {
178 std::cout << "Summary of command line options:" << std::endl
179 << "\tnum_mesh = " << n << std::endl
180 << "\trand_field = " << sg_rf_names[randField] << std::endl
181 << "\tmean = " << mean << std::endl
182 << "\tstd_dev = " << sigma << std::endl
183 << "\tnum_kl = " << num_KL << std::endl
184 << "\torder = " << p << std::endl
185 << "\tnormalize_basis = " << normalize_basis << std::endl
186 << "\tkrylov_method = " << krylov_method_names[krylov_method]
187 << std::endl
188 << "\tprec_strategy = " << prec_strategy_names[precStrategy]
189 << std::endl
190 << "\ttol = " << tol << std::endl
191 << "\tquadrature = " << sg_quad_names[quad_method]
192 << std::endl
193 << "\tsg_growth = " << sg_growth_names[sg_growth]
194 << std::endl;
195 }
196
197 bool nonlinear_expansion = false;
198 if (randField == LOGNORMAL)
199 nonlinear_expansion = true;
200
201 {
202 TEUCHOS_FUNC_TIME_MONITOR("Total Collocation Calculation Time");
203
204 // Create Stochastic Galerkin basis
205 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(num_KL);
206 for (int i=0; i<num_KL; i++) {
207 Teuchos::RCP<Stokhos::OneDOrthogPolyBasis<int,double> > b;
208 if (randField == UNIFORM) {
209 b = Teuchos::rcp(new Stokhos::LegendreBasis<int,double>(
210 p, normalize_basis));
211 }
212 else if (randField == CC_UNIFORM) {
213 b =
215 p, normalize_basis, true));
216 }
217 else if (randField == RYS) {
218 b = Teuchos::rcp(new Stokhos::RysBasis<int,double>(
219 p, weightCut, normalize_basis));
220 }
221 else if (randField == LOGNORMAL) {
222 b = Teuchos::rcp(new Stokhos::HermiteBasis<int,double>(
223 p, normalize_basis));
224 }
225 bases[i] = b;
226 }
227 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
228 Teuchos::rcp(new Stokhos::CompletePolynomialBasis<int,double>(bases));
229
230 // Create sparse grid
231 Teuchos::RCP<Stokhos::Quadrature<int,double> > quad;
232 if (quad_method == SPARSE_GRID) {
233#ifdef HAVE_STOKHOS_DAKOTA
234 int sparse_grid_growth = Pecos::MODERATE_RESTRICTED_GROWTH;
235 if (sg_growth == SLOW_RESTRICTED)
236 sparse_grid_growth = Pecos::SLOW_RESTRICTED_GROWTH;
237 else if (sg_growth == MODERATE_RESTRICTED)
238 sparse_grid_growth = Pecos::MODERATE_RESTRICTED_GROWTH;
239 else if (sg_growth == UNRESTRICTED)
240 sparse_grid_growth = Pecos::UNRESTRICTED_GROWTH;
241 quad = Teuchos::rcp(new Stokhos::SparseGridQuadrature<int,double>(
242 basis,p,1e-12,sparse_grid_growth));
243#else
244 TEUCHOS_TEST_FOR_EXCEPTION(
245 true, std::logic_error,
246 "Sparse grids require building with Dakota support!");
247#endif
248 }
249 else if (quad_method == TENSOR) {
250 quad =
251 Teuchos::rcp(new Stokhos::TensorProductQuadrature<int,double>(basis));
252 }
253 const Teuchos::Array< Teuchos::Array<double> >& quad_points =
254 quad->getQuadPoints();
255 const Teuchos::Array<double>& quad_weights = quad->getQuadWeights();
256 int nqp = quad_weights.size();
257
258 // Create application
259 twoD_diffusion_ME model(Comm, n, num_KL, sigma, mean,
260 basis, nonlinear_expansion);
261
262 // Model data
263 Teuchos::RCP<Epetra_Vector> p =
264 Teuchos::rcp(new Epetra_Vector(*(model.get_p_map(0))));
265 Teuchos::RCP<Epetra_Vector> x =
266 Teuchos::rcp(new Epetra_Vector(*(model.get_x_map())));
267 Teuchos::RCP<Epetra_Vector> f =
268 Teuchos::rcp(new Epetra_Vector(*(model.get_f_map())));
269 Teuchos::RCP<Epetra_Vector> g =
270 Teuchos::rcp(new Epetra_Vector(*(model.get_g_map(0))));
271 Teuchos::RCP<Epetra_CrsMatrix> A =
272 Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(model.create_W());
273 EpetraExt::ModelEvaluator::InArgs inArgs = model.createInArgs();
274 EpetraExt::ModelEvaluator::OutArgs outArgs = model.createOutArgs();
275 EpetraExt::ModelEvaluator::OutArgs outArgs2 = model.createOutArgs();
276
277 // Data to compute/store mean & variance
278 Epetra_Vector x2(*(model.get_x_map()));
279 Epetra_Vector x_mean(*(model.get_x_map()));
280 Epetra_Vector x_var(*(model.get_x_map()));
281 Epetra_Vector g2(*(model.get_g_map(0)));
282 Epetra_Vector g_mean(*(model.get_g_map(0)));
283 Epetra_Vector g_var(*(model.get_g_map(0)));
284
285 // Setup preconditioner
286 Teuchos::ParameterList precParams;
287 precParams.set("default values", "SA");
288 precParams.set("ML output", 0);
289 precParams.set("max levels",5);
290 precParams.set("increasing or decreasing","increasing");
291 precParams.set("aggregation: type", "Uncoupled");
292 precParams.set("smoother: type","ML symmetric Gauss-Seidel");
293 precParams.set("smoother: sweeps",2);
294 precParams.set("smoother: pre or post", "both");
295 precParams.set("coarse: max size", 200);
296#ifdef HAVE_ML_AMESOS
297 precParams.set("coarse: type","Amesos-KLU");
298#else
299 precParams.set("coarse: type","Jacobi");
300#endif
301 bool checkFiltering = false;
302 if (precStrategy == REUSE) {
303 checkFiltering = true;
304 precParams.set("reuse: enable", true);
305 }
306 Teuchos::RCP<ML_Epetra::MultiLevelPreconditioner> M =
307 Teuchos::rcp(new ML_Epetra::MultiLevelPreconditioner(*A, precParams,
308 false));
309 if (precStrategy == MEAN) {
310 TEUCHOS_FUNC_TIME_MONITOR("Deterministic Preconditioner Calculation");
311 *A = *(model.get_mean());
312 M->ComputePreconditioner();
313 }
314
315 // Setup AztecOO solver
316 AztecOO aztec;
317 if (krylov_method == GMRES)
318 aztec.SetAztecOption(AZ_solver, AZ_gmres);
319 else if (krylov_method == CG)
320 aztec.SetAztecOption(AZ_solver, AZ_cg);
321 aztec.SetAztecOption(AZ_precond, AZ_none);
322 aztec.SetAztecOption(AZ_kspace, 100);
323 aztec.SetAztecOption(AZ_conv, AZ_r0);
324 aztec.SetAztecOption(AZ_output, 0);
325 aztec.SetUserOperator(A.get());
326 aztec.SetPrecOperator(M.get());
327 aztec.SetLHS(x.get());
328 aztec.SetRHS(f.get());
329
330 x_mean.PutScalar(0.0);
331 x_var.PutScalar(0.0);
332 // Loop over colloction points
333 for (int qp=0; qp<nqp; qp++) {
334 TEUCHOS_FUNC_TIME_MONITOR("Collocation Loop");
335
336 if (qp%100 == 0 || qp == nqp-1)
337 std::cout << "Collocation point " << qp+1 <<'/' << nqp << "\n";
338
339 // Set parameters
340 for (int i=0; i<num_KL; i++)
341 (*p)[i] = quad_points[qp][i];
342
343 // Set in/out args
344 inArgs.set_p(0, p);
345 inArgs.set_x(x);
346 outArgs.set_f(f);
347 outArgs.set_W(A);
348
349 // Evaluate model at collocation point
350 x->PutScalar(0.0);
351 model.evalModel(inArgs, outArgs);
352 f->Scale(-1.0);
353
354 // Compute preconditioner
355 if (precStrategy != MEAN) {
356 TEUCHOS_FUNC_TIME_MONITOR("Deterministic Preconditioner Calculation");
357 M->ComputePreconditioner(checkFiltering);
358 }
359
360 // Solve linear system
361 {
362 TEUCHOS_FUNC_TIME_MONITOR("Deterministic Solve");
363 aztec.Iterate(1000, tol);
364 }
365
366 // Compute responses
367 outArgs2.set_g(0, g);
368 model.evalModel(inArgs, outArgs2);
369
370 // Sum contributions to mean and variance
371 x2.Multiply(1.0, *x, *x, 0.0);
372 g2.Multiply(1.0, *g, *g, 0.0);
373 x_mean.Update(quad_weights[qp], *x, 1.0);
374 x_var.Update(quad_weights[qp], x2, 1.0);
375 g_mean.Update(quad_weights[qp], *g, 1.0);
376 g_var.Update(quad_weights[qp], g2, 1.0);
377
378 }
379 x2.Multiply(1.0, x_mean, x_mean, 0.0);
380 g2.Multiply(1.0, g_mean, g_mean, 0.0);
381 x_var.Update(-1.0, x2, 1.0);
382 g_var.Update(-1.0, g2, 1.0);
383
384 // Compute standard deviations
385 for (int i=0; i<x_var.MyLength(); i++)
386 x_var[i] = std::sqrt(x_var[i]);
387 for (int i=0; i<g_var.MyLength(); i++)
388 g_var[i] = std::sqrt(g_var[i]);
389
390 std::cout.precision(16);
391 std::cout << "\nResponse Mean = " << std::endl << g_mean << std::endl;
392 std::cout << "Response Std. Dev. = " << std::endl << g_var << std::endl;
393
394 // Save mean and variance to file
395 EpetraExt::VectorToMatrixMarketFile("mean_col.mm", x_mean);
396 EpetraExt::VectorToMatrixMarketFile("std_dev_col.mm", x_var);
397
398 }
399
400 Teuchos::TimeMonitor::summarize(std::cout);
401 Teuchos::TimeMonitor::zeroOutTimers();
402
403 }
404
405 catch (std::exception& e) {
406 std::cout << e.what() << std::endl;
407 }
408 catch (std::string& s) {
409 std::cout << s << std::endl;
410 }
411 catch (char *s) {
412 std::cout << s << std::endl;
413 }
414 catch (...) {
415 std::cout << "Caught unknown exception!" <<std:: endl;
416 }
417
418#ifdef HAVE_MPI
419 MPI_Finalize() ;
420#endif
421
422}
int MyLength() const
int Multiply(char TransA, char TransB, double ScalarAB, const Epetra_MultiVector &A, const Epetra_MultiVector &B, double ScalarThis)
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
int PutScalar(double ScalarConstant)
Legendre polynomial basis using Clenshaw-Curtis quadrature points.
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Hermite polynomial basis.
Legendre polynomial basis.
Rys polynomial basis.
Defines quadrature for a tensor product basis by tensor products of 1-D quadrature rules.
ModelEvaluator for a linear 2-D diffusion problem.
Teuchos::RCP< Epetra_Operator > create_W() const
Create W = alpha*M + beta*J matrix.
InArgs createInArgs() const
Create InArgs.
Teuchos::RCP< const Epetra_Map > get_f_map() const
Return residual vector map.
OutArgs createOutArgs() const
Create OutArgs.
Teuchos::RCP< Epetra_CrsMatrix > get_mean() const
Get mean matrix.
Teuchos::RCP< const Epetra_Map > get_g_map(int j) const
Return response function map.
void evalModel(const InArgs &inArgs, const OutArgs &outArgs) const
Evaluate model on InArgs.
Teuchos::RCP< const Epetra_Map > get_p_map(int l) const
Return parameter vector map.
Teuchos::RCP< const Epetra_Map > get_x_map() const
Return solution vector map.
ScalarType g(const Teuchos::Array< ScalarType > &x, const ScalarType &y)
ScalarType f(const Teuchos::Array< ScalarType > &x, double a, double b)
const int num_krylov_method
const int num_sg_quad
int main(int argc, char *argv[])
const PrecStrategy prec_strategy_values[]
const int num_sg_growth
const int num_prec_strategy
const SG_Quad sg_quad_values[]
const char * sg_quad_names[]
const char * krylov_method_names[]
const Krylov_Method krylov_method_values[]
const char * sg_rf_names[]
const char * sg_growth_names[]
const int num_sg_rf
const SG_GROWTH sg_growth_values[]
const SG_RF sg_rf_values[]
const char * prec_strategy_names[]