Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
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
68 const int num_krylov_method = 2;
70 const char *krylov_method_names[] = { "GMRES", "CG" };
71 
72 // Random field types
74 const int num_sg_rf = 4;
76 const char *sg_rf_names[] = { "Uniform", "CC-Uniform", "Rys", "Log-Normal" };
77 
78 // Quadrature types
80 const int num_sg_quad = 2;
82 const char *sg_quad_names[] = { "tensor", "sparse-grid" };
83 
84 // Sparse grid growth rules
86 const int num_sg_growth = 3;
89 const char *sg_growth_names[] = {
90  "Slow Restricted", "Moderate Restricted", "Unrestricted" };
91 
92 int 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) {
214  b = rcp(new Stokhos::RysBasis<int,double>(
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) {
246  quad = rcp(new Stokhos::TensorProductQuadrature<int,double>(basis));
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  }
321  A_mp.ReplaceGlobalValues(
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 }
KOKKOS_INLINE_FUNCTION PCE< Storage > sqrt(const PCE< Storage > &a)
Hermite polynomial basis.
int MyLength() const
Copy
const Krylov_Method krylov_method_values[]
Teuchos::RCP< const Epetra_Map > get_f_map() const
Return residual vector map.
const char * sg_quad_names[]
int InsertGlobalIndices(int_type GlobalRow, int NumIndices, int_type *Indices)
int FillComplete(bool OptimizeDataStorage=true)
const char * sg_rf_names[]
const char * krylov_method_names[]
int Scale(double ScalarValue)
Teuchos::RCP< const Epetra_Map > get_g_map(int j) const
Return response function map.
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
Definition: csr_vector.h:260
ModelEvaluator for a linear 2-D diffusion problem.
Rys polynomial basis.
int main(int argc, char *argv[])
Teuchos::RCP< const Epetra_Map > get_x_map() const
Return solution vector map.
const char * sg_growth_names[]
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
const SG_RF sg_rf_values[]
const SG_GROWTH sg_growth_values[]
Legendre polynomial basis.
const SG_Quad sg_quad_values[]
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
void evalModel(const InArgs &inArgs, const OutArgs &outArgs) const
Evaluate model on InArgs.
ScalarType f(const Teuchos::Array< ScalarType > &x, double a, double b)
Legendre polynomial basis using Clenshaw-Curtis quadrature points.
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_p_map(int l) const
Return parameter vector map.
OutArgs createOutArgs() const
Create OutArgs.
virtual int ReplaceGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
Defines quadrature for a tensor product basis by tensor products of 1-D quadrature rules...
ScalarType g(const Teuchos::Array< ScalarType > &x, const ScalarType &y)