Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
division_example.cpp
Go to the documentation of this file.
1// $Id$
2// $Source$
3// @HEADER
4// ***********************************************************************
5//
6// Stokhos Package
7// Copyright (2009) Sandia Corporation
8//
9// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
10// license for use of this work by or on behalf of the U.S. Government.
11//
12// Redistribution and use in source and binary forms, with or without
13// modification, are permitted provided that the following conditions are
14// met:
15//
16// 1. Redistributions of source code must retain the above copyright
17// notice, this list of conditions and the following disclaimer.
18//
19// 2. Redistributions in binary form must reproduce the above copyright
20// notice, this list of conditions and the following disclaimer in the
21// documentation and/or other materials provided with the distribution.
22//
23// 3. Neither the name of the Corporation nor the names of the
24// contributors may be used to endorse or promote products derived from
25// this software without specific prior written permission.
26//
27// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38//
39// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
40//
41// ***********************************************************************
42// @HEADER
43
44// hermite_example
45//
46// usage:
47// hermite_example
48//
49// output:
50// prints the Hermite Polynomial Chaos Expansion of the simple function
51//
52// v = 1/(log(u)^2+1)
53//
54// where u = 1 + 0.4*H_1(x) + 0.06*H_2(x) + 0.002*H_3(x), x is a zero-mean
55// and unit-variance Gaussian random variable, and H_i(x) is the i-th
56// Hermite polynomial.
57
58#include "Stokhos.hpp"
59#include "Stokhos_Sacado.hpp"
60#include "Teuchos_GlobalMPISession.hpp"
61#include "Teuchos_CommandLineProcessor.hpp"
62
63
64// Typename of PC expansion type
66
67// Linear solvers
69 const int num_division_solver = 3;
71 const char *division_solver_names[] = { "Dense_Direct", "GMRES", "CG" };
72
73// Preconditioners
75 const int num_division_prec = 5;
77 const char *division_prec_names[] = { "None", "Diag", "Jacobi","GS","Schur"};
78
79// Option for Schur complement precond: full or diag D
81 const int num_schur_option = 2;
83 const char *schur_option_names[] = { "full", "diag"};
84
85// Full matrix or linear matrix (pb = dim + 1 ) used for preconditioner
87 const int num_prec_option = 2;
89 const char *prec_option_names[] = { "full", "linear"};
90
91
92int main(int argc, char **argv)
93{
94 using Teuchos::Array;
95 using Teuchos::RCP;
96 using Teuchos::rcp;
97
98 // If applicable, set up MPI.
99 Teuchos::GlobalMPISession mpiSession (&argc, &argv);
100
101 try {
102
103 // Setup command line options
104 Teuchos::CommandLineProcessor CLP;
105 CLP.setDocString("This example tests the / operator.\n");
106 int d = 3;
107 CLP.setOption("dim", &d, "Stochastic dimension");
108 int p = 5;
109 CLP.setOption("order", &p, "Polynomial order");
110 int n = 100;
111 CLP.setOption("samples", &n, "Number of samples");
112 double shift = 5.0;
113 CLP.setOption("shift", &shift, "Shift point");
114 double tolerance = 1e-6;
115 CLP.setOption("tolerance", &tolerance, "Tolerance in Iterative Solver");
116 int prec_level = 1;
117 CLP.setOption("prec_level", &prec_level, "Level in Schur Complement Prec 0->Solve A0u0=g0 with division; 1->Form 1x1 Schur Complement");
118 int max_it_div = 50;
119 CLP.setOption("max_it_div", &max_it_div, "Maximum # of iterations for division iterative solver");
120 bool equilibrate = true;
121 CLP.setOption("equilibrate", "noequilibrate", &equilibrate,
122 "Equilibrate the linear system");
123
124 Division_Solver solve_method = Dense_Direct;
125 CLP.setOption("solver", &solve_method,
127 "Solver Method");
128 Division_Prec prec_method = None;
129 CLP.setOption("prec", &prec_method,
131 "Preconditioner Method");
132 Schur_option schur_option = diag;
133 CLP.setOption("schur_option", &schur_option,
135 "Schur option");
136 Prec_option prec_option = whole;
137 CLP.setOption("prec_option", &prec_option,
139 "Prec option");
140 CLP.parse( argc, argv );
141
142 // Basis
143 Array< RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
144 for (int i=0; i<d; i++) {
145 bases[i] = rcp(new Stokhos::LegendreBasis<int,double>(p));
146 }
147 RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
149
150 // Quadrature method
151 RCP<const Stokhos::Quadrature<int,double> > quad =
153
154 // Triple product tensor
155 RCP<Stokhos::Sparse3Tensor<int,double> > Cijk =
156 basis->computeTripleProductTensor();
157
158 // Expansion methods
159 Teuchos::RCP<Stokhos::QuadOrthogPolyExpansion<int,double> > quad_expn =
161 basis, Cijk, quad));
162 Teuchos::RCP<Teuchos::ParameterList> alg_params =
163 Teuchos::rcp(new Teuchos::ParameterList);
164
165 alg_params->set("Division Tolerance", tolerance);
166 alg_params->set("prec_iter", prec_level);
167 alg_params->set("max_it_div", max_it_div);
168 if (solve_method == Dense_Direct)
169 alg_params->set("Division Strategy", "Dense Direct");
170 else if (solve_method == GMRES)
171 alg_params->set("Division Strategy", "GMRES");
172 else if (solve_method == CG)
173 alg_params->set("Division Strategy", "CG");
174
175
176 if (prec_method == None)
177 alg_params->set("Prec Strategy", "None");
178 else if (prec_method == Diag)
179 alg_params->set("Prec Strategy", "Diag");
180 else if (prec_method == Jacobi)
181 alg_params->set("Prec Strategy", "Jacobi");
182 else if (prec_method == GS)
183 alg_params->set("Prec Strategy", "GS");
184 else if (prec_method == Schur)
185 alg_params->set("Prec Strategy", "Schur");
186
187 if (schur_option == diag)
188 alg_params->set("Schur option", "diag");
189 else
190 alg_params->set("Schur option", "full");
191 if (prec_option == linear)
192 alg_params->set("Prec option", "linear");
193
194 alg_params->set("Use Quadrature for Division", false);
195
196 if (equilibrate)
197 alg_params->set("Equilibrate", 1);
198 else
199 alg_params->set("Equilibrate", 0);
200
201
202
203 Teuchos::RCP<Stokhos::QuadOrthogPolyExpansion<int,double> > alg_expn =
205 basis, Cijk, quad, alg_params));
206
207 // Polynomial expansions
208 pce_type u_quad(quad_expn), v_quad(quad_expn);
209 u_quad.term(0,0) = 0.0;
210 for (int i=0; i<d; i++) {
211 u_quad.term(i,1) = 1.0;
212 }
213 pce_type u_alg(alg_expn), v_alg(alg_expn);
214 u_alg.term(0,0) = 0.0;
215 for (int i=0; i<d; i++) {
216 u_alg.term(i,1) = 1.0;
217 }
218
219 // Compute expansion
220 double scale = std::exp(shift);
221 pce_type b_alg = std::exp(shift + u_alg)/scale;
222 pce_type b_quad = std::exp(shift + u_quad)/scale;
223// v_alg = (1.0/scale) / b_alg;
224// v_quad = (1.0/scale) / b_quad;
225 v_alg = 1.0 / std::exp(shift + u_alg);
226 v_quad = 1.0 /std::exp(shift + u_quad);
227
228// std::cout << b_alg.getOrthogPolyApprox() << std::endl;
229
230 // Print u and v
231// std::cout << "quadrature: v = 1.0 / (shift + u) = ";
232// v_quad.print(std::cout);
233// std::cout << "dense solve: v = 1.0 / (shift + u) = ";
234// v_alg.print(std::cout);
235
236 double h = 2.0 / (n-1);
237 double err_quad = 0.0;
238 double err_alg = 0.0;
239 for (int i=0; i<n; i++) {
240
241 double x = -1.0 + h*i;
242 Array<double> pt(d);
243 for (int j=0; j<d; j++)
244 pt[j] = x;
245 double up = u_quad.evaluate(pt);
246 double vp = 1.0/(shift+up);
247 double vp_quad = v_quad.evaluate(pt);
248 double vp_alg = v_alg.evaluate(pt);
249 // std::cout << "vp = " << vp_quad << std::endl;
250 // std::cout << "vp_quad = " << vp_quad << std::endl;
251 // std::cout << "vp_alg = " << vp_alg << std::endl;
252 double point_err_quad = std::abs(vp-vp_quad);
253 double point_err_alg = std::abs(vp-vp_alg);
254 if (point_err_quad > err_quad) err_quad = point_err_quad;
255 if (point_err_alg > err_alg) err_alg = point_err_alg;
256 }
257 std::cout << "\tL_infty norm of quadrature error = " << err_quad
258 << std::endl;
259 std::cout << "\tL_infty norm of solver error = " << err_alg
260 << std::endl;
261
262 // Check the answer
263 //if (std::abs(err) < 1e-2)
264 std::cout << "\nExample Passed!" << std::endl;
265
266 Teuchos::TimeMonitor::summarize(std::cout);
267 }
268 catch (std::exception& e) {
269 std::cout << e.what() << std::endl;
270 }
271}
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Legendre polynomial basis.
Orthogonal polynomial expansions based on numerical quadrature.
Defines quadrature for a tensor product basis by tensor products of 1-D quadrature rules.
const char * prec_option_names[]
const int num_division_prec
const int num_prec_option
const char * division_prec_names[]
int main(int argc, char **argv)
const int num_schur_option
const Division_Prec Division_prec_values[]
const Prec_option Prec_option_values[]
const int num_division_solver
Sacado::PCE::OrthogPoly< double, Stokhos::StandardStorage< int, double > > pce_type
Division_Prec
const char * division_solver_names[]
const Schur_option Schur_option_values[]
Division_Solver
@ Dense_Direct
const Division_Solver Division_solver_values[]
const char * schur_option_names[]