Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_LogNormalUnitTest.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// Test the exponential of a linear Hermite expansion using an analytic
45// closed form solution
46
47#include "Teuchos_UnitTestHarness.hpp"
48#include "Teuchos_TestingHelpers.hpp"
49#include "Teuchos_UnitTestRepository.hpp"
50#include "Teuchos_GlobalMPISession.hpp"
51
52#include "Stokhos.hpp"
54
55TEUCHOS_UNIT_TEST( Stokhos_LogNormal_TP, UnitTest ) {
56 // Basis of dimension 3, order 5
57 const int d = 3;
58 const int p = 5;
59 const double mean = 0.2;
60 const double std_dev = 0.1;
61 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
62 for (int i=0; i<d; i++) {
63 bases[i] = Teuchos::rcp(new Stokhos::HermiteBasis<int,double>(p));
64 }
65 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
66 Teuchos::rcp(new Stokhos::CompletePolynomialBasis<int,double>(bases));
67
68 // Quadrature method
69 Teuchos::RCP<const Stokhos::Quadrature<int,double> > quad =
70 Teuchos::rcp(new Stokhos::TensorProductQuadrature<int,double>(basis));
71
72 // Triple product tensor
73 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk =
74 basis->computeTripleProductTensor();
75
76 // Expansion method
78
79 // Polynomial expansions
80 Stokhos::OrthogPolyApprox<int,double> u(basis), v(basis), w(basis);
81 u.term(0,0) = mean;
82 for (int i=0; i<d; i++)
83 u.term(i,1) = std_dev / (i+1);
84
85 // Compute expansion
86 expn.exp(v,u);
87
88 // Compute expansion analytically
89 double w_mean = mean;
90 for (int j=0; j<d; j++)
91 w_mean += u.term(j,1)*u.term(j,1)/2.0;
92 w_mean = std::exp(w_mean);
93 for (int i=0; i<basis->size(); i++) {
94 Stokhos::MultiIndex<int> multiIndex = basis->term(i);
95 double s = 1.0;
96 for (int j=0; j<d; j++)
97 s *= std::pow(u.term(j,1), multiIndex[j]);
98 w[i] = w_mean*s/basis->norm_squared(i);
99 }
100
101 success = Stokhos::comparePCEs(v, "quad_expansion", w, "analytic",
102 1e-12, 1e-9, out);
103}
104
105#ifdef HAVE_STOKHOS_DAKOTA
106TEUCHOS_UNIT_TEST( Stokhos_LogNormal_SG, UnitTest ) {
107 // Basis of dimension 3, order 5
108 const int d = 3;
109 const int p = 5;
110 const double mean = 0.2;
111 const double std_dev = 0.1;
112 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
113 for (int i=0; i<d; i++) {
114 bases[i] = Teuchos::rcp(new Stokhos::HermiteBasis<int,double>(p));
115 }
116 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
117 Teuchos::rcp(new Stokhos::CompletePolynomialBasis<int,double>(bases));
118
119 // Quadrature method
120 Teuchos::RCP<const Stokhos::Quadrature<int,double> > quad =
121 Teuchos::rcp(new Stokhos::SparseGridQuadrature<int,double>(basis, p));
122
123 // Triple product tensor
124 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk =
125 basis->computeTripleProductTensor();
126
127 // Expansion method
128 Stokhos::QuadOrthogPolyExpansion<int,double> expn(basis, Cijk, quad);
129
130 // Polynomial expansions
131 Stokhos::OrthogPolyApprox<int,double> u(basis), v(basis), w(basis);
132 u.term(0,0) = mean;
133 for (int i=0; i<d; i++)
134 u.term(i,1) = std_dev / (i+1);
135
136 // Compute expansion
137 expn.exp(v,u);
138
139 // Compute expansion analytically
140 double w_mean = mean;
141 for (int j=0; j<d; j++)
142 w_mean += u.term(j,1)*u.term(j,1)/2.0;
143 w_mean = std::exp(w_mean);
144 for (int i=0; i<basis->size(); i++) {
145 Stokhos::MultiIndex<int> multiIndex = basis->term(i);
146 double s = 1.0;
147 for (int j=0; j<d; j++)
148 s *= std::pow(u.term(j,1), multiIndex[j]);
149 w[i] = w_mean*s/basis->norm_squared(i);
150 }
151
152 success = Stokhos::comparePCEs(v, "quad_expansion", w, "analytic",
153 1e-12, 1e-9, out);
154}
155#endif
156
157int main( int argc, char* argv[] ) {
158 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
159 return Teuchos::UnitTestRepository::runUnitTestsFromMain(argc, argv);
160}
int main(int argc, char *argv[])
TEUCHOS_UNIT_TEST(Stokhos_LogNormal_TP, UnitTest)
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Hermite polynomial basis.
A multidimensional index.
Class to store coefficients of a projection onto an orthogonal polynomial basis.
Orthogonal polynomial expansions based on numerical quadrature.
Defines quadrature for a tensor product basis by tensor products of 1-D quadrature rules.
bool comparePCEs(const PCEType &a1, const std::string &a1_name, const Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &a2, const std::string &a2_name, const ValueType &rel_tol, const ValueType &abs_tol, Teuchos::FancyOStream &out)