Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
sacado_ensemble_example.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#include "Stokhos_Sacado.hpp"
44#include "Teuchos_CommandLineProcessor.hpp"
45
46// The function to compute the polynomial chaos expansion of,
47// written as a template function
48template <class ScalarType>
49ScalarType simple_function(const ScalarType& u) {
50 ScalarType z = std::log(u);
51 return 1.0/(z*z + 1.0);
52}
53
54int main(int argc, char **argv)
55{
56 // Typename of Polynomial Chaos scalar type
57 typedef Stokhos::StandardStorage<int,double> pce_storage_type;
59
60 // Typename of ensemble scalar type
61 const int EnsembleSize = 8;
64
65 // Short-hand for several classes used below
66 using Teuchos::Array;
67 using Teuchos::RCP;
68 using Teuchos::rcp;
79
80 try {
81
82 // Setup command line options
83 Teuchos::CommandLineProcessor CLP;
84 CLP.setDocString(
85 "This example computes the PC expansion of a simple function.\n");
86 int p = 4;
87 CLP.setOption("order", &p, "Polynomial order");
88 bool sparse = false;
89 CLP.setOption("sparse", "tensor", &sparse,
90 "Use sparse grid or tensor product quadrature");
91
92 // Parse arguments
93 CLP.parse( argc, argv );
94
95 // Basis of dimension 3, order given by command-line option
96 const int d = 3;
97 Array< RCP<const OneDOrthogPolyBasis<int,double> > > bases(d);
98 for (int i=0; i<d; i++) {
99 bases[i] = rcp(new HermiteBasis<int,double>(p, true));
100 }
101 RCP<const CompletePolynomialBasis<int,double> > basis =
102 rcp(new CompletePolynomialBasis<int,double>(bases));
103 const int pce_size = basis->size();
104 std::cout << "basis size = " << pce_size << std::endl;
105
106 // Quadrature method
107 RCP<const Quadrature<int,double> > quad;
108 if (sparse) {
109 const TotalOrderIndexSet<int> index_set(d, p);
110 quad = rcp(new SmolyakSparseGridQuadrature<int,double>(basis, index_set));
111 }
112 else {
113 quad = rcp(new TensorProductQuadrature<int,double>(basis));
114 }
115 std::cout << "quadrature size = " << quad->size() << std::endl;
116
117 // Triple product tensor
118 RCP<Sparse3Tensor<int,double> > Cijk =
119 basis->computeTripleProductTensor();
120
121 // Expansion method
122 RCP<QuadOrthogPolyExpansion<int,double> > expn =
123 rcp(new QuadOrthogPolyExpansion<int,double>(basis, Cijk, quad));
124
125 // Polynomial expansion of u (note: these are coefficients in the
126 // normalized basis)
127 pce_type u(expn);
128 u.term(0,0) = 1.0; // zeroth order term
129 u.term(0,1) = 0.1; // first order term for dimension 0
130 u.term(1,1) = 0.05; // first order term for dimension 1
131 u.term(2,1) = 0.01; // first order term for dimension 2
132
133 //
134 // Compute PCE expansion of function using NISP with ensemble propagation
135 //
136
137 // Extract quadrature data
138 const int num_quad_points = quad->size();
139 const Array<double>& quad_weights = quad->getQuadWeights();
140 const Array< Array<double> >& quad_points = quad->getQuadPoints();
141 const Array< Array<double> >& quad_values = quad->getBasisAtQuadPoints();
142
143 // Loop over quadrature points in blocks of size EnsembleSize
144 pce_type v(expn);
145 ensemble_type u_ensemble;
146 for (int qp_block=0; qp_block<num_quad_points; qp_block+=EnsembleSize) {
147 const int qp_sz = qp_block+EnsembleSize <= num_quad_points ?
148 EnsembleSize : num_quad_points-qp_block;
149
150 // Evaluate u at each quadrature point
151 for (int qp=0; qp<qp_sz; ++qp)
152 u_ensemble.fastAccessCoeff(qp) =
153 u.evaluate(quad_points[qp_block+qp], quad_values[qp_block+qp]);
154 for (int qp=qp_sz; qp<EnsembleSize; ++qp)
155 u_ensemble.fastAccessCoeff(qp) = u_ensemble.fastAccessCoeff(qp_sz-1);
156
157 // Evaluate function at each quadrature point
158 ensemble_type v_ensemble = simple_function(u_ensemble);
159
160 // Sum results into PCE integral
161 for (int pc=0; pc<pce_size; ++pc)
162 for (int qp=0; qp<qp_sz; ++qp)
163 v.fastAccessCoeff(pc) += v_ensemble.fastAccessCoeff(qp)*quad_weights[qp_block+qp]*quad_values[qp_block+qp][pc];
164 }
165
166 /*
167 for (int qp=0; qp<num_quad_points; ++qp) {
168 double u_qp = u.evaluate(quad_points[qp]);
169 double v_qp = simple_function(u_qp);
170 double w = quad_weights[qp];
171 for (int pc=0; pc<pce_size; ++pc)
172 v.fastAccessCoeff(pc) += v_qp*w*quad_values[qp][pc];
173 }
174 */
175
176 // Print u and v
177 std::cout << "\tu = ";
178 u.print(std::cout);
179 std::cout << "\tv = ";
180 v.print(std::cout);
181
182 // Compute moments
183 double mean = v.mean();
184 double std_dev = v.standard_deviation();
185
186 // Evaluate PCE and function at a point = 0.25 in each dimension
187 Teuchos::Array<double> pt(d);
188 for (int i=0; i<d; i++)
189 pt[i] = 0.25;
190 double up = u.evaluate(pt);
191 double vp = simple_function(up);
192 double vp2 = v.evaluate(pt);
193
194 // Print results
195 std::cout << "\tv mean = " << mean << std::endl;
196 std::cout << "\tv std. dev. = " << std_dev << std::endl;
197 std::cout << "\tv(0.25) (true) = " << vp << std::endl;
198 std::cout << "\tv(0.25) (pce) = " << vp2 << std::endl;
199
200 // Check the answer
201 if (std::abs(vp - vp2) < 1e-2)
202 std::cout << "\nExample Passed!" << std::endl;
203 }
204 catch (std::exception& e) {
205 std::cout << e.what() << std::endl;
206 }
207}
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Hermite polynomial basis.
Legendre polynomial basis.
Abstract base class for 1-D orthogonal polynomials.
Orthogonal polynomial expansions based on numerical quadrature.
Abstract base class for quadrature methods.
Defines quadrature for a tensor product basis by Smolyak sparse grids.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
Statically allocated storage class.
Defines quadrature for a tensor product basis by tensor products of 1-D quadrature rules.
An isotropic total order index set.
Sacado::ETPCE::OrthogPoly< double, Stokhos::StandardStorage< int, double > > pce_type
int main(int argc, char **argv)
ScalarType simple_function(const ScalarType &u)