Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
linear2d_diffusion_multipoint_commuted.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// Block utilities
63#include "EpetraExt_BlockUtility.h"
64#include "EpetraExt_RowMatrixOut.h"
65
66// Krylov methods
68const int num_krylov_method = 2;
70const char *krylov_method_names[] = { "GMRES", "CG" };
71
72// Random field types
74const int num_sg_rf = 4;
76const char *sg_rf_names[] = { "Uniform", "CC-Uniform", "Rys", "Log-Normal" };
77
78// Quadrature types
80const int num_sg_quad = 2;
82const char *sg_quad_names[] = { "tensor", "sparse-grid" };
83
84// Sparse grid growth rules
86const int num_sg_growth = 3;
89const char *sg_growth_names[] = {
90 "Slow Restricted", "Moderate Restricted", "Unrestricted" };
91
92int main(int argc, char *argv[]) {
93 using Teuchos::RCP;
94 using Teuchos::rcp;
95 using Teuchos::rcp_dynamic_cast;
96 using Teuchos::Array;
97 using Teuchos::CommandLineProcessor;
98 using Teuchos::TimeMonitor;
99 using Teuchos::ParameterList;
100 using EpetraExt::BlockUtility;
101 using EpetraExt::ModelEvaluator;
102
103// Initialize MPI
104#ifdef HAVE_MPI
105 MPI_Init(&argc,&argv);
106#endif
107
108 // Create a communicator for Epetra objects
109 RCP<Epetra_Comm> Comm;
110#ifdef HAVE_MPI
111 Comm = rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
112#else
113 Comm = rcp(new Epetra_SerialComm);
114#endif
115
116 int MyPID = Comm->MyPID();
117
118 try {
119
120 // Setup command line options
121 CommandLineProcessor CLP;
122 CLP.setDocString(
123 "This example runs a stochastic collocation method.\n");
124
125 int n = 32;
126 CLP.setOption("num_mesh", &n, "Number of mesh points in each direction");
127
128 SG_RF randField = UNIFORM;
129 CLP.setOption("rand_field", &randField,
131 "Random field type");
132
133 double mean = 0.2;
134 CLP.setOption("mean", &mean, "Mean");
135
136 double sigma = 0.1;
137 CLP.setOption("std_dev", &sigma, "Standard deviation");
138
139 double weightCut = 1.0;
140 CLP.setOption("weight_cut", &weightCut, "Weight cut");
141
142 int num_KL = 2;
143 CLP.setOption("num_kl", &num_KL, "Number of KL terms");
144
145 int p = 3;
146 CLP.setOption("order", &p, "Polynomial order");
147
148 bool normalize_basis = true;
149 CLP.setOption("normalize", "unnormalize", &normalize_basis,
150 "Normalize PC basis");
151
152 Krylov_Method krylov_method = CG;
153 CLP.setOption("krylov_method", &krylov_method,
155 "Krylov method");
156
157 double tol = 1e-12;
158 CLP.setOption("tol", &tol, "Solver tolerance");
159
160#ifdef HAVE_STOKHOS_DAKOTA
161 SG_Quad quad_method = SPARSE_GRID;
162#else
163 SG_Quad quad_method = TENSOR;
164#endif
165 CLP.setOption("quadrature", &quad_method,
167 "Quadrature type");
168
169 SG_GROWTH sg_growth = MODERATE_RESTRICTED;
170 CLP.setOption("sg_growth", &sg_growth,
172 "Sparse grid growth rule");
173
174 CLP.parse( argc, argv );
175
176 if (MyPID == 0) {
177 std::cout << "Summary of command line options:" << std::endl
178 << "\tnum_mesh = " << n << std::endl
179 << "\trand_field = " << sg_rf_names[randField] << std::endl
180 << "\tmean = " << mean << std::endl
181 << "\tstd_dev = " << sigma << std::endl
182 << "\tnum_kl = " << num_KL << std::endl
183 << "\torder = " << p << std::endl
184 << "\tnormalize_basis = " << normalize_basis << std::endl
185 << "\tkrylov_method = " << krylov_method_names[krylov_method]
186 << std::endl
187 << "\ttol = " << tol << std::endl
188 << "\tquadrature = " << sg_quad_names[quad_method]
189 << std::endl
190 << "\tsg_growth = " << sg_growth_names[sg_growth]
191 << std::endl;
192 }
193
194 bool nonlinear_expansion = false;
195 if (randField == LOGNORMAL)
196 nonlinear_expansion = true;
197
198 {
199 TEUCHOS_FUNC_TIME_MONITOR("Total Collocation Calculation Time");
200
201 // Create Stochastic Galerkin basis
202 Array< RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(num_KL);
203 for (int i=0; i<num_KL; i++) {
204 RCP<Stokhos::OneDOrthogPolyBasis<int,double> > b;
205 if (randField == UNIFORM) {
206 b = rcp(new Stokhos::LegendreBasis<int,double>(p, normalize_basis));
207 }
208 else if (randField == CC_UNIFORM) {
209 b =
211 p, normalize_basis, true));
212 }
213 else if (randField == RYS) {
215 p, weightCut, normalize_basis));
216 }
217 else if (randField == LOGNORMAL) {
219 p, normalize_basis));
220 }
221 bases[i] = b;
222 }
223 RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
225
226 // Create sparse grid
227 RCP<Stokhos::Quadrature<int,double> > quad;
228 if (quad_method == SPARSE_GRID) {
229#ifdef HAVE_STOKHOS_DAKOTA
230 int sparse_grid_growth = Pecos::MODERATE_RESTRICTED_GROWTH;
231 if (sg_growth == SLOW_RESTRICTED)
232 sparse_grid_growth = Pecos::SLOW_RESTRICTED_GROWTH;
233 else if (sg_growth == MODERATE_RESTRICTED)
234 sparse_grid_growth = Pecos::MODERATE_RESTRICTED_GROWTH;
235 else if (sg_growth == UNRESTRICTED)
236 sparse_grid_growth = Pecos::UNRESTRICTED_GROWTH;
237 quad = rcp(new Stokhos::SparseGridQuadrature<int,double>(
238 basis,p,1e-12,sparse_grid_growth));
239#else
240 TEUCHOS_TEST_FOR_EXCEPTION(
241 true, std::logic_error,
242 "Sparse grids require building with Dakota support!");
243#endif
244 }
245 else if (quad_method == TENSOR) {
247 }
248 const Array< Array<double> >& quad_points = quad->getQuadPoints();
249 const Array<double>& quad_weights = quad->getQuadWeights();
250 int nqp = quad_weights.size();
251
252 // Create application
253 twoD_diffusion_ME model(Comm, n, num_KL, sigma, mean,
254 basis, nonlinear_expansion);
255 RCP<Epetra_Vector> p = rcp(new Epetra_Vector(*(model.get_p_map(0))));
256 RCP<Epetra_Vector> x = rcp(new Epetra_Vector(*(model.get_x_map())));
257 RCP<Epetra_Vector> f = rcp(new Epetra_Vector(*(model.get_f_map())));
258 RCP<Epetra_CrsMatrix> A =
259 rcp_dynamic_cast<Epetra_CrsMatrix>(model.create_W());
260
261 // Build commuted multipoint map and operator
262 RCP<const Epetra_Map> spatial_map = model.get_x_map();
263 Epetra_LocalMap stochastic_map(nqp, 0, *Comm);
264 RCP<const Epetra_Map> product_map =
265 rcp(BlockUtility::GenerateBlockMap(stochastic_map, *spatial_map, *Comm),
266 false);
267 const Epetra_CrsGraph& spatial_graph = A->Graph();
268 Epetra_CrsGraph stochastic_graph(Copy, stochastic_map, 1, true);
269 for (int i=0; i<nqp; ++i)
270 stochastic_graph.InsertGlobalIndices(i, 1, &i);
271 stochastic_graph.FillComplete();
272 RCP<const Epetra_CrsGraph> product_graph =
273 rcp(BlockUtility::GenerateBlockGraph(stochastic_graph, spatial_graph, *Comm));
274 Epetra_CrsMatrix A_mp(Copy, *product_graph);
275 Epetra_Vector f_mp(*product_map);
276 Epetra_Vector x_mp(*product_map);
277
278 if (MyPID == 0) {
279 std::cout << "spatial size = " << spatial_map->NumGlobalElements()
280 << " quadrature size = " << nqp
281 << " multipoint size = " << product_map->NumGlobalElements()
282 << std::endl;
283 }
284
285 // Loop over colloction points and assemble global operator, RHS
286 int max_num_entries = A->MaxNumEntries();
287 Array<int> block_indices(max_num_entries);
288 double *values;
289 int *indices;
290 int num_entries;
291 for (int qp=0; qp<nqp; qp++) {
292 TEUCHOS_FUNC_TIME_MONITOR("Compute MP Matrix, RHS");
293
294 // Set parameters
295 for (int i=0; i<num_KL; i++)
296 (*p)[i] = quad_points[qp][i];
297
298 // Set in/out args
299 ModelEvaluator::InArgs inArgs = model.createInArgs();
300 ModelEvaluator::OutArgs outArgs = model.createOutArgs();
301 inArgs.set_p(0, p);
302 inArgs.set_x(x);
303 outArgs.set_f(f);
304 outArgs.set_W(A);
305
306 // Evaluate model at collocation point
307 x->PutScalar(0.0);
308 model.evalModel(inArgs, outArgs);
309
310 // Put f, A into global objects
311 const int num_rows = spatial_map->NumMyElements();
312 for (int local_row=0; local_row<num_rows; ++local_row) {
313 const int global_row = spatial_map->GID(local_row);
314 f_mp[nqp*local_row+qp] = (*f)[local_row];
315 A->ExtractMyRowView(local_row, num_entries, values, indices);
316 for (int j=0; j<num_entries; ++j) {
317 int local_col = indices[j];
318 int global_col = A->GCID(local_col);
319 block_indices[j] = nqp*global_col+qp;
320 }
322 nqp*global_row+qp, num_entries, values, &block_indices[0]);
323 }
324
325 }
326 f_mp.Scale(-1.0);
327 A_mp.FillComplete();
328
329 //EpetraExt::RowMatrixToMatrixMarketFile("A_mp.mm", A_mp);
330
331 // Setup preconditioner
332 ParameterList precParams;
333 precParams.set("default values", "SA");
334 precParams.set("ML output", 10);
335 precParams.set("max levels",5);
336 precParams.set("increasing or decreasing","increasing");
337 precParams.set("aggregation: type", "Uncoupled");
338 precParams.set("smoother: type","ML symmetric Gauss-Seidel");
339 precParams.set("smoother: sweeps",2);
340 precParams.set("smoother: pre or post", "both");
341 precParams.set("coarse: max size", 200);
342 precParams.set("PDE equations",nqp);
343#ifdef HAVE_ML_AMESOS
344 precParams.set("coarse: type","Amesos-KLU");
345#else
346 precParams.set("coarse: type","Jacobi");
347#endif
348 RCP<ML_Epetra::MultiLevelPreconditioner> M_mp;
349 {
350 TEUCHOS_FUNC_TIME_MONITOR("Preconditioner Setup");
351 M_mp = rcp(new ML_Epetra::MultiLevelPreconditioner(A_mp, precParams));
352 }
353
354 // Setup AztecOO solver
355 AztecOO aztec;
356 if (krylov_method == GMRES)
357 aztec.SetAztecOption(AZ_solver, AZ_gmres);
358 else if (krylov_method == CG)
359 aztec.SetAztecOption(AZ_solver, AZ_cg);
360 aztec.SetAztecOption(AZ_precond, AZ_none);
361 aztec.SetAztecOption(AZ_kspace, 100);
362 aztec.SetAztecOption(AZ_conv, AZ_r0);
363 aztec.SetAztecOption(AZ_output, 10);
364 aztec.SetUserOperator(&A_mp);
365 aztec.SetPrecOperator(M_mp.get());
366 aztec.SetLHS(&x_mp);
367 aztec.SetRHS(&f_mp);
368
369 // Solve linear system
370 {
371 TEUCHOS_FUNC_TIME_MONITOR("System Solve");
372 aztec.Iterate(1000, tol);
373 }
374
375 // Compute responses
376 RCP<Epetra_Vector> g = rcp(new Epetra_Vector(*(model.get_g_map(0))));
377 Epetra_Vector g2(*(model.get_g_map(0)));
378 Epetra_Vector g_mean(*(model.get_g_map(0)));
379 Epetra_Vector g_var(*(model.get_g_map(0)));
380 for (int qp=0; qp<nqp; qp++) {
381 TEUCHOS_FUNC_TIME_MONITOR("Compute Responses");
382
383 // Set parameters
384 for (int i=0; i<num_KL; i++)
385 (*p)[i] = quad_points[qp][i];
386
387 // Extract x
388 const int num_rows = spatial_map->NumMyElements();
389 for (int local_row=0; local_row<num_rows; ++local_row)
390 (*x)[local_row] = x_mp[nqp*local_row+qp];
391
392 // Set in/out args
393 ModelEvaluator::InArgs inArgs = model.createInArgs();
394 ModelEvaluator::OutArgs outArgs = model.createOutArgs();
395 inArgs.set_p(0, p);
396 inArgs.set_x(x);
397 outArgs.set_g(0,g);
398
399 // Evaluate model at collocation point
400 model.evalModel(inArgs, outArgs);
401
402 // Sum contributions to mean and variance
403 g2.Multiply(1.0, *g, *g, 0.0);
404 g_mean.Update(quad_weights[qp], *g, 1.0);
405 g_var.Update(quad_weights[qp], g2, 1.0);
406 }
407 g2.Multiply(1.0, g_mean, g_mean, 0.0);
408 g_var.Update(-1.0, g2, 1.0);
409
410 // Compute standard deviations
411 for (int i=0; i<g_var.MyLength(); i++)
412 g_var[i] = std::sqrt(g_var[i]);
413
414 std::cout.precision(16);
415 std::cout << "\nResponse Mean = " << std::endl << g_mean << std::endl;
416 std::cout << "Response Std. Dev. = " << std::endl << g_var << std::endl;
417
418 }
419
420 TimeMonitor::summarize(std::cout);
421 TimeMonitor::zeroOutTimers();
422
423 }
424
425 catch (std::exception& e) {
426 std::cout << e.what() << std::endl;
427 }
428 catch (std::string& s) {
429 std::cout << s << std::endl;
430 }
431 catch (char *s) {
432 std::cout << s << std::endl;
433 }
434 catch (...) {
435 std::cout << "Caught unknown exception!" <<std:: endl;
436 }
437
438#ifdef HAVE_MPI
439 MPI_Finalize() ;
440#endif
441
442}
Copy
int InsertGlobalIndices(int_type GlobalRow, int NumIndices, int_type *Indices)
int FillComplete(bool OptimizeDataStorage=true)
virtual int ReplaceGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
int Scale(double ScalarValue)
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)
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< 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)
int main(int argc, char *argv[])
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 SG_GROWTH sg_growth_values[]
const SG_RF sg_rf_values[]