Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
epetra/linear2d_diffusion_pce.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_Epetra.hpp"
57
58// Timing utilities
59#include "Teuchos_TimeMonitor.hpp"
60
61// I/O utilities
62#include "EpetraExt_VectorOut.h"
63
64int main(int argc, char *argv[]) {
65 int n = 32; // spatial discretization (per dimension)
66 int num_KL = 2; // number of KL terms
67 int p = 3; // polynomial order
68 double mu = 0.2; // mean of exponential random field
69 double s = 0.1; // std. dev. of exponential r.f.
70 bool nonlinear_expansion = false; // nonlinear expansion of diffusion coeff
71 // (e.g., log-normal)
72 bool symmetric = false; // use symmetric formulation
73
74 double g_mean_exp = 0.172988; // expected response mean
75 double g_std_dev_exp = 0.0380007; // expected response std. dev.
76 double g_tol = 1e-6; // tolerance on determining success
77
78// Initialize MPI
79#ifdef HAVE_MPI
80 MPI_Init(&argc,&argv);
81#endif
82
83 int MyPID;
84
85 try {
86
87 {
88 TEUCHOS_FUNC_TIME_MONITOR("Total PCE Calculation Time");
89
90 // Create a communicator for Epetra objects
91 Teuchos::RCP<const Epetra_Comm> globalComm;
92#ifdef HAVE_MPI
93 globalComm = Teuchos::rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
94#else
95 globalComm = Teuchos::rcp(new Epetra_SerialComm);
96#endif
97 MyPID = globalComm->MyPID();
98
99 // Create Stochastic Galerkin basis and expansion
100 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(num_KL);
101 for (int i=0; i<num_KL; i++)
102 bases[i] = Teuchos::rcp(new Stokhos::LegendreBasis<int,double>(p,true));
103 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
105 1e-12));
106 int sz = basis->size();
107 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk;
108 if (nonlinear_expansion)
109 Cijk = basis->computeTripleProductTensor();
110 else
111 Cijk = basis->computeLinearTripleProductTensor();
112 Teuchos::RCP<Stokhos::OrthogPolyExpansion<int,double> > expansion =
114 Cijk));
115 if (MyPID == 0)
116 std::cout << "Stochastic Galerkin expansion size = " << sz << std::endl;
117
118 // Create stochastic parallel distribution
119 int num_spatial_procs = -1;
120 if (argc > 1)
121 num_spatial_procs = std::atoi(argv[1]);
122 Teuchos::ParameterList parallelParams;
123 parallelParams.set("Number of Spatial Processors", num_spatial_procs);
124 // parallelParams.set("Rebalance Stochastic Graph", true);
125 // Teuchos::ParameterList& isorropia_params =
126 // parallelParams.sublist("Isorropia");
127 // isorropia_params.set("Balance objective", "nonzeros");
128 Teuchos::RCP<Stokhos::ParallelData> sg_parallel_data =
129 Teuchos::rcp(new Stokhos::ParallelData(basis, Cijk, globalComm,
130 parallelParams));
131 Teuchos::RCP<const EpetraExt::MultiComm> sg_comm =
132 sg_parallel_data->getMultiComm();
133 Teuchos::RCP<const Epetra_Comm> app_comm =
134 sg_parallel_data->getSpatialComm();
135
136 // Create application
137 Teuchos::RCP<twoD_diffusion_ME> model =
138 Teuchos::rcp(new twoD_diffusion_ME(app_comm, n, num_KL, s, mu, basis,
139 nonlinear_expansion, symmetric));
140
141 // Setup stochastic Galerkin algorithmic parameters
142 Teuchos::RCP<Teuchos::ParameterList> sgParams =
143 Teuchos::rcp(new Teuchos::ParameterList);
144 Teuchos::ParameterList& sgOpParams =
145 sgParams->sublist("SG Operator");
146 Teuchos::ParameterList& sgPrecParams =
147 sgParams->sublist("SG Preconditioner");
148 if (!nonlinear_expansion) {
149 sgParams->set("Parameter Expansion Type", "Linear");
150 sgParams->set("Jacobian Expansion Type", "Linear");
151 }
152 sgOpParams.set("Operator Method", "Matrix Free");
153 sgPrecParams.set("Preconditioner Method", "Mean-based");
154 //sgPrecParams.set("Preconditioner Method", "Approximate Gauss-Seidel");
155 //sgPrecParams.set("Preconditioner Method", "Approximate Schur Complement");
156 sgPrecParams.set("Symmetric Gauss-Seidel", symmetric);
157 sgPrecParams.set("Mean Preconditioner Type", "ML");
158 Teuchos::ParameterList& precParams =
159 sgPrecParams.sublist("Mean Preconditioner Parameters");
160 precParams.set("default values", "SA");
161 precParams.set("ML output", 0);
162 precParams.set("max levels",5);
163 precParams.set("increasing or decreasing","increasing");
164 precParams.set("aggregation: type", "Uncoupled");
165 precParams.set("smoother: type","ML symmetric Gauss-Seidel");
166 precParams.set("smoother: sweeps",2);
167 precParams.set("smoother: pre or post", "both");
168 precParams.set("coarse: max size", 200);
169#ifdef HAVE_ML_AMESOS
170 precParams.set("coarse: type","Amesos-KLU");
171#else
172 precParams.set("coarse: type","Jacobi");
173#endif
174
175 // Create stochastic Galerkin model evaluator
176 Teuchos::RCP<Stokhos::SGModelEvaluator> sg_model =
177 Teuchos::rcp(new Stokhos::SGModelEvaluator(model, basis, Teuchos::null,
178 expansion, sg_parallel_data,
179 sgParams));
180
181 // Set up stochastic parameters
182 // The current implementation of the model doesn't actually use these
183 // values, but is hard-coded to certain uncertainty models
184 Teuchos::Array<double> point(num_KL, 1.0);
185 Teuchos::Array<double> basis_vals(sz);
186 basis->evaluateBases(point, basis_vals);
187 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_p_poly =
188 sg_model->create_p_sg(0);
189 for (int i=0; i<num_KL; i++) {
190 sg_p_poly->term(i,0)[i] = 0.0;
191 sg_p_poly->term(i,1)[i] = 1.0 / basis_vals[i+1];
192 }
193 std::cout << *sg_p_poly << std::endl;
194
195 // Create vectors and operators
196 Teuchos::RCP<const Epetra_Vector> sg_p = sg_p_poly->getBlockVector();
197 Teuchos::RCP<Epetra_Vector> sg_x =
198 Teuchos::rcp(new Epetra_Vector(*(sg_model->get_x_map())));
199 sg_x->PutScalar(0.0);
200 Teuchos::RCP<Epetra_Vector> sg_f =
201 Teuchos::rcp(new Epetra_Vector(*(sg_model->get_f_map())));
202 Teuchos::RCP<Epetra_Vector> sg_dx =
203 Teuchos::rcp(new Epetra_Vector(*(sg_model->get_x_map())));
204 Teuchos::RCP<Epetra_Operator> sg_J = sg_model->create_W();
205 Teuchos::RCP<Epetra_Operator> sg_M = sg_model->create_WPrec()->PrecOp;
206
207
208 // Setup InArgs and OutArgs
209 EpetraExt::ModelEvaluator::InArgs sg_inArgs = sg_model->createInArgs();
210 EpetraExt::ModelEvaluator::OutArgs sg_outArgs = sg_model->createOutArgs();
211 sg_inArgs.set_p(1, sg_p);
212 sg_inArgs.set_x(sg_x);
213 sg_outArgs.set_f(sg_f);
214 sg_outArgs.set_W(sg_J);
215 sg_outArgs.set_WPrec(sg_M);
216
217 // Evaluate model
218 sg_model->evalModel(sg_inArgs, sg_outArgs);
219
220 // Print initial residual norm
221 double norm_f;
222 sg_f->Norm2(&norm_f);
223 if (MyPID == 0)
224 std::cout << "\nInitial residual norm = " << norm_f << std::endl;
225
226 // Setup AztecOO solver
227 AztecOO aztec;
228 if (symmetric)
229 aztec.SetAztecOption(AZ_solver, AZ_cg);
230 else
231 aztec.SetAztecOption(AZ_solver, AZ_gmres);
232 aztec.SetAztecOption(AZ_precond, AZ_none);
233 aztec.SetAztecOption(AZ_kspace, 20);
234 aztec.SetAztecOption(AZ_conv, AZ_r0);
235 aztec.SetAztecOption(AZ_output, 1);
236 aztec.SetUserOperator(sg_J.get());
237 aztec.SetPrecOperator(sg_M.get());
238 aztec.SetLHS(sg_dx.get());
239 aztec.SetRHS(sg_f.get());
240
241 // Solve linear system
242 aztec.Iterate(1000, 1e-12);
243
244 // Update x
245 sg_x->Update(-1.0, *sg_dx, 1.0);
246
247 // Save solution to file
248 EpetraExt::VectorToMatrixMarketFile("stochastic_solution.mm", *sg_x);
249
250 // Save RHS to file
251 EpetraExt::VectorToMatrixMarketFile("stochastic_RHS.mm", *sg_f);
252
253 // Save operator to file (set "Operator Method" to "Fully Assembled" above)
254 Teuchos::RCP<Epetra_CrsMatrix> sg_J_crs =
255 Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(sg_J);
256 if (sg_J_crs != Teuchos::null)
257 EpetraExt::RowMatrixToMatrixMarketFile("stochastic_operator.mm",
258 *sg_J_crs);
259
260 // Save mean and variance to file
261 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_x_poly =
262 sg_model->create_x_sg(View, sg_x.get());
263 Epetra_Vector mean(*(model->get_x_map()));
264 Epetra_Vector std_dev(*(model->get_x_map()));
265 sg_x_poly->computeMean(mean);
266 sg_x_poly->computeStandardDeviation(std_dev);
267 EpetraExt::VectorToMatrixMarketFile("mean_gal.mm", mean);
268 EpetraExt::VectorToMatrixMarketFile("std_dev_gal.mm", std_dev);
269
270 // Compute new residual & response function
271 EpetraExt::ModelEvaluator::OutArgs sg_outArgs2 = sg_model->createOutArgs();
272 Teuchos::RCP<Epetra_Vector> sg_g =
273 Teuchos::rcp(new Epetra_Vector(*(sg_model->get_g_map(0))));
274 sg_f->PutScalar(0.0);
275 sg_outArgs2.set_f(sg_f);
276 sg_outArgs2.set_g(0, sg_g);
277 sg_model->evalModel(sg_inArgs, sg_outArgs2);
278
279 // Print initial residual norm
280 sg_f->Norm2(&norm_f);
281 if (MyPID == 0)
282 std::cout << "\nFinal residual norm = " << norm_f << std::endl;
283
284 // Print mean and standard deviation of responses
285 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_g_poly =
286 sg_model->create_g_sg(0, View, sg_g.get());
287 Epetra_Vector g_mean(*(model->get_g_map(0)));
288 Epetra_Vector g_std_dev(*(model->get_g_map(0)));
289 sg_g_poly->computeMean(g_mean);
290 sg_g_poly->computeStandardDeviation(g_std_dev);
291 std::cout.precision(16);
292 // std::cout << "\nResponse Expansion = " << std::endl;
293 // std::cout.precision(12);
294 // sg_g_poly->print(std::cout);
295 std::cout << std::endl;
296 std::cout << "Response Mean = " << std::endl << g_mean << std::endl;
297 std::cout << "Response Std. Dev. = " << std::endl << g_std_dev << std::endl;
298
299 // Determine if example passed
300 bool passed = false;
301 if (norm_f < 1.0e-10 &&
302 std::abs(g_mean[0]-g_mean_exp) < g_tol &&
303 std::abs(g_std_dev[0]-g_std_dev_exp) < g_tol)
304 passed = true;
305 if (MyPID == 0) {
306 if (passed)
307 std::cout << "Example Passed!" << std::endl;
308 else
309 std::cout << "Example Failed!" << std::endl;
310 }
311
312 }
313
314 Teuchos::TimeMonitor::summarize(std::cout);
315 Teuchos::TimeMonitor::zeroOutTimers();
316
317 }
318
319 catch (std::exception& e) {
320 std::cout << e.what() << std::endl;
321 }
322 catch (std::string& s) {
323 std::cout << s << std::endl;
324 }
325 catch (char *s) {
326 std::cout << s << std::endl;
327 }
328 catch (...) {
329 std::cout << "Caught unknown exception!" << std::endl;
330 }
331
332#ifdef HAVE_MPI
333 MPI_Finalize() ;
334#endif
335
336}
Orthogonal polynomial expansions limited to algebraic operations.
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Legendre polynomial basis.
Nonlinear, stochastic Galerkin ModelEvaluator.
ModelEvaluator for a linear 2-D diffusion problem.
int main(int argc, char *argv[])