Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
linear2d_diffusion.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 // Class implementing our problem
44 
45 // Epetra communicator
46 #ifdef HAVE_MPI
47 #include "Epetra_MpiComm.h"
48 #else
49 #include "Epetra_SerialComm.h"
50 #endif
51 
52 // solver
53 #include "Ifpack2_Factory.hpp"
54 #include "BelosLinearProblem.hpp"
55 #include "BelosTpetraAdapter.hpp"
56 #include "BelosPseudoBlockCGSolMgr.hpp"
57 #include "BelosPseudoBlockGmresSolMgr.hpp"
58 #include "MatrixMarket_Tpetra.hpp"
59 
60 // Stokhos Stochastic Galerkin
61 #include "Stokhos_Epetra.hpp"
62 
63 // Timing utilities
64 #include "Teuchos_TimeMonitor.hpp"
65 
66 int main(int argc, char *argv[]) {
67  typedef double Scalar;
68  typedef double MeshScalar;
69  typedef double BasisScalar;
70  typedef int LocalOrdinal;
71  typedef int GlobalOrdinal;
72  typedef Tpetra::DefaultPlatform::DefaultPlatformType::NodeType Node;
73 
74  using Teuchos::RCP;
75  using Teuchos::rcp;
76  using Teuchos::ParameterList;
77 
78  LocalOrdinal n = 32; // spatial discretization (per dimension)
79  LocalOrdinal num_KL = 2; // number of KL terms
80  LocalOrdinal p = 3; // polynomial order
81  BasisScalar mu = 0.2; // mean of exponential random field
82  BasisScalar s = 0.1; // std. dev. of exponential r.f.
83  bool nonlinear_expansion = false; // nonlinear expansion of diffusion coeff
84  // (e.g., log-normal)
85  bool symmetric = false; // use symmetric formulation
86 
87 // Initialize MPI
88 #ifdef HAVE_MPI
89  MPI_Init(&argc,&argv);
90 #endif
91 
92  LocalOrdinal MyPID;
93 
94  try {
95 
96  {
97  TEUCHOS_FUNC_TIME_MONITOR("Total PCE Calculation Time");
98 
99  // Create a communicator for Epetra objects
100  RCP<const Epetra_Comm> globalComm;
101 #ifdef HAVE_MPI
102  globalComm = rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
103 #else
104  globalComm = rcp(new Epetra_SerialComm);
105 #endif
106  MyPID = globalComm->MyPID();
107 
108  // Create Stochastic Galerkin basis and expansion
109  Teuchos::Array< RCP<const Stokhos::OneDOrthogPolyBasis<LocalOrdinal,BasisScalar> > > bases(num_KL);
110  for (LocalOrdinal i=0; i<num_KL; i++)
111  bases[i] = rcp(new Stokhos::LegendreBasis<LocalOrdinal,BasisScalar>(p,true));
112  RCP<const Stokhos::CompletePolynomialBasis<LocalOrdinal,BasisScalar> > basis =
114  1e-12));
115  LocalOrdinal sz = basis->size();
116  RCP<Stokhos::Sparse3Tensor<LocalOrdinal,BasisScalar> > Cijk;
117  if (nonlinear_expansion)
118  Cijk = basis->computeTripleProductTensor();
119  else
120  Cijk = basis->computeLinearTripleProductTensor();
121  RCP<Stokhos::OrthogPolyExpansion<LocalOrdinal,BasisScalar> > expansion =
123  Cijk));
124  if (MyPID == 0)
125  std::cout << "Stochastic Galerkin expansion size = " << sz << std::endl;
126 
127  // Create stochastic parallel distribution
128  LocalOrdinal num_spatial_procs = -1;
129  if (argc > 1)
130  num_spatial_procs = std::atoi(argv[1]);
131  ParameterList parallelParams;
132  parallelParams.set("Number of Spatial Processors", num_spatial_procs);
133  // parallelParams.set("Rebalance Stochastic Graph", true);
134  // Teuchos::ParameterList& isorropia_params =
135  // parallelParams.sublist("Isorropia");
136  // isorropia_params.set("Balance objective", "nonzeros");
137  RCP<Stokhos::ParallelData> sg_parallel_data =
138  rcp(new Stokhos::ParallelData(basis, Cijk, globalComm,
139  parallelParams));
140  RCP<const EpetraExt::MultiComm> sg_comm =
141  sg_parallel_data->getMultiComm();
142  RCP<const Epetra_Comm> app_comm =
143  sg_parallel_data->getSpatialComm();
144 
145  // Create Teuchos::Comm from Epetra_Comm
146  RCP< Teuchos::Comm<int> > teuchos_app_comm;
147 #ifdef HAVE_MPI
148  RCP<const Epetra_MpiComm> app_mpi_comm =
149  Teuchos::rcp_dynamic_cast<const Epetra_MpiComm>(app_comm);
150  RCP<const Teuchos::OpaqueWrapper<MPI_Comm> > raw_mpi_comm =
151  Teuchos::opaqueWrapper(app_mpi_comm->Comm());
152  teuchos_app_comm = rcp(new Teuchos::MpiComm<int>(raw_mpi_comm));
153 #else
154  teuchos_app_comm = rcp(new Teuchos::SerialComm<int>());
155 #endif
156 
157  // Create application
159  RCP<problem_type> model =
160  rcp(new problem_type(teuchos_app_comm, n, num_KL, s, mu,
161  nonlinear_expansion, symmetric));
162 
163 
164  // Create vectors and operators
165  typedef problem_type::Tpetra_Vector Tpetra_Vector;
166  typedef problem_type::Tpetra_CrsMatrix Tpetra_CrsMatrix;
167  typedef Tpetra::MatrixMarket::Writer<Tpetra_CrsMatrix> Writer;
168  RCP<const Tpetra_Vector> p = model->get_p_init(0);
169  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(model->get_x_map());
170  x->putScalar(0.0);
171  RCP<Tpetra_Vector> f = Tpetra::createVector<Scalar>(model->get_f_map());
172  RCP<Tpetra_Vector> dx = Tpetra::createVector<Scalar>(model->get_x_map());
173  RCP<Tpetra_CrsMatrix> J = model->create_W();
174 
175  // Create preconditioner
176  Teuchos::ParameterList precParams;
177  std::string prec_name = "RILUK";
178  precParams.set("fact: iluk level-of-fill", 1);
179  precParams.set("fact: iluk level-of-overlap", 0);
180  typedef Ifpack2::Preconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tprec;
181  Teuchos::RCP<Tprec> M;
182  Ifpack2::Factory factory;
183  M = factory.create<Tpetra_CrsMatrix>(prec_name, J);
184  M->setParameters(precParams);
185 
186  // Evaluate model
187  model->computeResidual(*x, *p, *f);
188  model->computeJacobian(*x, *p, *J);
189  M->initialize();
190  M->compute();
191 
192  // Print initial residual norm
193  Scalar norm_f;
194  f->norm2(Teuchos::arrayView(&norm_f,1));
195  if (MyPID == 0)
196  std::cout << "\nInitial residual norm = " << norm_f << std::endl;
197 
198  // Setup Belos solver
199  RCP<ParameterList> belosParams = rcp(new ParameterList);
200  belosParams->set("Num Blocks", 20);
201  belosParams->set("Convergence Tolerance", 1e-12);
202  belosParams->set("Maximum Iterations", 1000);
203  belosParams->set("Verbosity", 33);
204  belosParams->set("Output Style", 1);
205  belosParams->set("Output Frequency", 1);
206  typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
208  typedef Belos::OperatorTraits<Scalar,MV,OP> BOPT;
209  typedef Belos::MultiVecTraits<Scalar,MV> BMVT;
210  typedef Belos::LinearProblem<Scalar,MV,OP> BLinProb;
211  RCP< BLinProb > problem = rcp(new BLinProb(J, dx, f));
212  problem->setRightPrec(M);
213  problem->setProblem();
214  RCP<Belos::SolverManager<Scalar,MV,OP> > solver;
215  if (symmetric)
216  solver = rcp(new Belos::PseudoBlockCGSolMgr<Scalar,MV,OP>(problem,
217  belosParams));
218  else
219  solver = rcp(new Belos::PseudoBlockGmresSolMgr<Scalar,MV,OP>(problem,
220  belosParams));
221 
222  // Solve linear system
223  Belos::ReturnType ret = solver->solve();
224 
225  // Update x
226  x->update(-1.0, *dx, 1.0);
227 
228  // Compute new residual & response function
229  RCP<Tpetra_Vector> g = Tpetra::createVector<Scalar>(model->get_g_map(0));
230  f->putScalar(0.0);
231  model->computeResidual(*x, *p, *f);
232  model->computeResponse(*x, *p, *g);
233 
234  // Print final residual norm
235  f->norm2(Teuchos::arrayView(&norm_f,1));
236  if (MyPID == 0)
237  std::cout << "\nFinal residual norm = " << norm_f << std::endl;
238 
239  // Print response
240  std::cout << "\nResponse = " << std::endl;
241  Writer::writeDense(std::cout, g);
242 
243  }
244 
245  Teuchos::TimeMonitor::summarize(std::cout);
246  Teuchos::TimeMonitor::zeroOutTimers();
247 
248  }
249 
250  catch (std::exception& e) {
251  std::cout << e.what() << std::endl;
252  }
253  catch (string& s) {
254  std::cout << s << std::endl;
255  }
256  catch (char *s) {
257  std::cout << s << std::endl;
258  }
259  catch (...) {
260  std::cout << "Caught unknown exception!" <<std:: endl;
261  }
262 
263 #ifdef HAVE_MPI
264  MPI_Finalize() ;
265 #endif
266 
267 }
int main(int argc, char *argv[])
Orthogonal polynomial expansions limited to algebraic operations.
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
Definition: csr_vector.h:260
KokkosClassic::DefaultNode::DefaultNodeType Node
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Legendre polynomial basis.
ScalarType f(const Teuchos::Array< ScalarType > &x, double a, double b)
pce_type Scalar
expr expr expr dx(i, j)
ScalarType g(const Teuchos::Array< ScalarType > &x, const ScalarType &y)
A linear 2-D diffusion problem.