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