Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
pce_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// pce_example
43//
44// usage:
45// pce_example
46//
47// output:
48// prints the Hermite Polynomial Chaos Expansion of the simple function
49//
50// v = 1/(log(u)^2+1)
51//
52// where u = 1 + 0.1*x_1 + 0.05*x_2 + 0.01*x_3 x1,x2,x3 are zero-mean
53// unit-variance Gaussian random variables.
54
55#include "Stokhos_Sacado.hpp"
56
57// The function to compute the polynomial chaos expansion of,
58// written as a template function
59template <class ScalarType>
60ScalarType simple_function(const ScalarType& u) {
61 ScalarType z = std::log(u);
62 return 1.0/(z*z + 1.0);
63}
64
65int main(int argc, char **argv)
66{
67 // Typename of Polynomial Chaos scalar type
70
71 // Short-hand for several classes used below
72 using Teuchos::Array;
73 using Teuchos::RCP;
74 using Teuchos::rcp;
82
83 try {
84
85 // Basis of dimension 3, order 4
86 const int d = 3;
87 const int p = 4;
88 Array< RCP<const OneDOrthogPolyBasis<int,double> > > bases(d);
89 for (int i=0; i<d; i++) {
90 bases[i] = rcp(new HermiteBasis<int,double>(p));
91 }
92 RCP<const CompletePolynomialBasis<int,double> > basis =
93 rcp(new CompletePolynomialBasis<int,double>(bases));
94
95 // Quadrature method
96 RCP<const Quadrature<int,double> > quad =
97 rcp(new TensorProductQuadrature<int,double>(basis));
98
99 // Triple product tensor
100 RCP<Sparse3Tensor<int,double> > Cijk =
101 basis->computeTripleProductTensor();
102
103 // Expansion method
104 RCP<QuadOrthogPolyExpansion<int,double> > expn =
105 rcp(new QuadOrthogPolyExpansion<int,double>(basis, Cijk, quad));
106
107 // Polynomial expansion of u
108 pce_type u(expn);
109 u.term(0,0) = 1.0; // zeroth order term
110 u.term(0,1) = 0.1; // first order term for dimension 0
111 u.term(1,1) = 0.05; // first order term for dimension 1
112 u.term(2,1) = 0.01; // first order term for dimension 2
113
114 // Compute PCE expansion of function
116
117 // Print u and v
118 std::cout << "\tu = ";
119 u.print(std::cout);
120 std::cout << "\tv = ";
121 v.print(std::cout);
122
123 // Compute moments
124 double mean = v.mean();
125 double std_dev = v.standard_deviation();
126
127 // Evaluate PCE and function at a point = 0.25 in each dimension
128 Teuchos::Array<double> pt(d);
129 for (int i=0; i<d; i++)
130 pt[i] = 0.25;
131 double up = u.evaluate(pt);
132 double vp = simple_function(up);
133 double vp2 = v.evaluate(pt);
134
135 // Print results
136 std::cout << "\tv mean = " << mean << std::endl;
137 std::cout << "\tv std. dev. = " << std_dev << std::endl;
138 std::cout << "\tv(0.25) (true) = " << vp << std::endl;
139 std::cout << "\tv(0.25) (pce) = " << vp2 << std::endl;
140 }
141 catch (std::exception& e) {
142 std::cout << e.what() << std::endl;
143 }
144}
Stokhos::StandardStorage< int, double > storage_type
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Hermite polynomial basis.
Abstract base class for 1-D orthogonal polynomials.
Orthogonal polynomial expansions based on numerical quadrature.
Abstract base class for quadrature methods.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
Defines quadrature for a tensor product basis by tensor products of 1-D quadrature rules.
Sacado::ETPCE::OrthogPoly< double, Stokhos::StandardStorage< int, double > > pce_type
int main(int argc, char **argv)
Definition: pce_example.cpp:65
ScalarType simple_function(const ScalarType &u)
Definition: pce_example.cpp:60