Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_ExponentialRandomFieldUnitTest.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 "Teuchos_UnitTestHarness.hpp"
43#include "Teuchos_TestingHelpers.hpp"
44#include "Teuchos_UnitTestRepository.hpp"
45#include "Teuchos_GlobalMPISession.hpp"
46
47#include "Stokhos.hpp"
48#include "Stokhos_KL_OneDExponentialCovarianceFunction.hpp"
49#include "Stokhos_KL_ExponentialRandomField.hpp"
50#include "Stokhos_UnitTestHelpers.hpp"
51
52namespace ExponentialRandomFieldUnitTest {
53
54 TEUCHOS_UNIT_TEST( Stokhos_ExponentialRandomField, OneD_Eigenvalue_Ordering ) {
55 int M = 100;
56 double a = -1.0;
57 double b = 1.5;
58 double L = 1.1;
59
60 // Setup covariance function
61 Teuchos::ParameterList solverParams;
62 solverParams.set("Nonlinear Solver Tolerance", 1e-8);
63 typedef Stokhos::KL::OneDExponentialCovarianceFunction<double> cov_type;
64 typedef cov_type::eigen_pair_type eigen_pair_type;
65 cov_type cov(M, a, b, L, 0, solverParams);
66
67 // Get eigenpairs
68 const Teuchos::Array<eigen_pair_type>& eigs = cov.getEigenPairs();
69
70 success = true;
71 out << std::endl;
72 for (int i=0; i<M-1; i++) {
73 out << "eigs[" << i << "] = " << eigs[i].eig_val << std::endl;
74 if (eigs[i].eig_val < eigs[i+1].eig_val)
75 success = false;
76 }
77 }
78
79 TEUCHOS_UNIT_TEST( Stokhos_ExponentialRandomField, OneD_Frequency_Spacing ) {
80 int M = 100;
81 double a = -1.0;
82 double b = 1.5;
83 double L = 1.1;
84
85 // Setup covariance function
86 Teuchos::ParameterList solverParams;
87 solverParams.set("Nonlinear Solver Tolerance", 1e-8);
88 typedef Stokhos::KL::OneDExponentialCovarianceFunction<double> cov_type;
89 typedef cov_type::eigen_pair_type eigen_pair_type;
90 cov_type cov(M, a, b, L, 0, solverParams);
91
92 // Get eigenpairs
93 const Teuchos::Array<eigen_pair_type>& eigs = cov.getEigenPairs();
94
95 success = true;
96 out << std::endl;
97 double pi = 4.0*std::atan(1.0);
98 for (int i=0; i<M-1; i++) {
99 double omega1 = eigs[i].eig_func.getFrequency();
100 double omega2 = eigs[i].eig_func.getFrequency();
101
102 out << "eigs[" << i << "].frequency = " << omega1 << std::endl;
103 if (omega2 - omega1 > pi)
104 success = false;
105 }
106 }
107
108 TEUCHOS_UNIT_TEST( Stokhos_ExponentialRandomField, OneD_Eigenfunction_Norm ) {
109 int p = 200;
110 int M = 10;
111 double a = -1.0;
112 double b = 1.5;
113 double L = 1.1;
114 double center = (b+a)/2.0;
115 double width = (b-a)/2.0;
116 double w_coeff = 2.0*width;
117
118 // Create basis for getting quadrature points
119 Stokhos::LegendreBasis<int,double> basis(p);
120 Teuchos::Array<double> quad_points, quad_weights;
121 Teuchos::Array< Teuchos::Array<double> > quad_values;
122 basis.getQuadPoints(p, quad_points, quad_weights, quad_values);
123
124 // Setup covariance function
125 Teuchos::ParameterList solverParams;
126 typedef Stokhos::KL::OneDExponentialCovarianceFunction<double> cov_type;
127 typedef cov_type::eigen_pair_type eigen_pair_type;
128 cov_type cov(M, a, b, L, 0, solverParams);
129
130 // Get eigenpairs
131 const Teuchos::Array<eigen_pair_type>& eigs = cov.getEigenPairs();
132
133 int nqp = quad_weights.size();
134 success = true;
135
136 out << std::endl;
137 // Loop over each eigenpair (lambda, b(x))
138 for (int i=0; i<M; i++) {
139
140 // compute \int_D b(x)*b(x) dx
141 double integral = 0.0;
142 double rhs = 1.0;
143 for (int qp=0; qp<nqp; qp++) {
144 double xp = center + quad_points[qp]*width;
145 double w = w_coeff*quad_weights[qp];
146 double val = eigs[i].eig_func.evaluate(xp);
147 integral += w*val*val;
148 }
149
150 out << "lambda = " << eigs[i].eig_val << ", integral = " << integral
151 << ", rhs = " << rhs << ", error = " << integral-rhs << std::endl;
152 success = success &&
153 Stokhos::compareValues(integral, "integral", rhs, "rhs",
154 1e-3, 1e-3, out);
155 }
156 }
157
158 TEUCHOS_UNIT_TEST( Stokhos_ExponentialRandomField, OneD_Eigenfunction_Orthogonality ) {
159 int p = 200;
160 int M = 10;
161 double a = -1.0;
162 double b = 1.5;
163 double L = 1.1;
164 double center = (b+a)/2.0;
165 double width = (b-a)/2.0;
166 double w_coeff = 2.0*width;
167
168 // Create basis for getting quadrature points
169 Stokhos::LegendreBasis<int,double> basis(p);
170 Teuchos::Array<double> quad_points, quad_weights;
171 Teuchos::Array< Teuchos::Array<double> > quad_values;
172 basis.getQuadPoints(p, quad_points, quad_weights, quad_values);
173
174 // Setup covariance function
175 Teuchos::ParameterList solverParams;
176 typedef Stokhos::KL::OneDExponentialCovarianceFunction<double> cov_type;
177 typedef cov_type::eigen_pair_type eigen_pair_type;
178 cov_type cov(M, a, b, L, 0, solverParams);
179
180 // Get eigenpairs
181 const Teuchos::Array<eigen_pair_type>& eigs = cov.getEigenPairs();
182
183 int nqp = quad_weights.size();
184 success = true;
185
186 out << std::endl;
187 for (int i=0; i<M; i++) {
188 for (int j=0; j<i; j++) {
189
190 // compute \int_D b_i(x)*b_j(x) dx
191 double integral = 0.0;
192 double rhs = 0.0;
193 for (int qp=0; qp<nqp; qp++) {
194 double xp = center + quad_points[qp]*width;
195 double w = w_coeff*quad_weights[qp];
196 double val1 = eigs[i].eig_func.evaluate(xp);
197 double val2 = eigs[j].eig_func.evaluate(xp);
198 integral += w*val1*val2;
199 }
200
201 out << "lambda = " << eigs[i].eig_val << ", integral = " << integral
202 << ", rhs = " << rhs << ", error = " << integral-rhs << std::endl;
203 success = success &&
204 Stokhos::compareValues(integral, "integral", rhs, "rhs",
205 1e-3, 1e-3, out);
206 }
207 }
208 }
209
210 TEUCHOS_UNIT_TEST( Stokhos_ExponentialRandomField, OneD_Eigen_Solution ) {
211 int p = 200;
212 int M = 10;
213 double a = -1.0;
214 double b = 1.5;
215 double L = 1.1;
216 double center = (b+a)/2.0;
217 double width = (b-a)/2.0;
218 double x = center + 0.25*width;
219 double w_coeff = 2.0*width;
220
221 // Create basis for getting quadrature points
222 Stokhos::LegendreBasis<int,double> basis(p);
223 Teuchos::Array<double> quad_points, quad_weights;
224 Teuchos::Array< Teuchos::Array<double> > quad_values;
225 basis.getQuadPoints(p, quad_points, quad_weights, quad_values);
226
227 // Setup covariance function
228 Teuchos::ParameterList solverParams;
229 typedef Stokhos::KL::OneDExponentialCovarianceFunction<double> cov_type;
230 typedef cov_type::eigen_pair_type eigen_pair_type;
231 cov_type cov(M, a, b, L, 0, solverParams);
232
233 // Get eigenpairs
234 const Teuchos::Array<eigen_pair_type>& eigs = cov.getEigenPairs();
235
236 int nqp = quad_weights.size();
237 success = true;
238
239 out << std::endl;
240 // Loop over each eigenpair (lambda, b(x))
241 for (int i=0; i<M; i++) {
242
243 // compute \int_D exp(-|x-x'|/L)b(x') dx'
244 double integral = 0.0;
245 for (int qp=0; qp<nqp; qp++) {
246 double xp = center + quad_points[qp]*width;
247 double w = w_coeff*quad_weights[qp];
248 integral +=
249 w*cov.evaluateCovariance(x,xp)*eigs[i].eig_func.evaluate(xp);
250 }
251
252 // compute lambda*b(x)
253 double rhs = eigs[i].eig_val*eigs[i].eig_func.evaluate(x);
254 out << "lambda = " << eigs[i].eig_val << ", integral = " << integral
255 << ", rhs = " << rhs << ", error = " << integral-rhs << std::endl;
256 success = success &&
257 Stokhos::compareValues(integral, "integral", rhs, "rhs",
258 1e-3, 1e-3, out);
259 }
260 }
261
262 TEUCHOS_UNIT_TEST( Stokhos_ExponentialRandomField, Product_Eigen_Solution ) {
263 // Create product basis
264 int p = 20;
265 int d = 2;
266 int M = 10;
267 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
268 for (int i=0; i<d; i++)
269 bases[i] = Teuchos::rcp(new Stokhos::LegendreBasis<int,double>(p));
270 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
271 Teuchos::rcp(new Stokhos::CompletePolynomialBasis<int,double>(bases));
272
273 // Tensor product quadrature
274 Stokhos::TensorProductQuadrature<int,double> quad(basis);
275 const Teuchos::Array< Teuchos::Array<double> >& quad_points =
276 quad.getQuadPoints();
277 const Teuchos::Array<double>& quad_weights = quad.getQuadWeights();
278
279 // Setup random field
280 Teuchos::ParameterList solverParams;
281 solverParams.set("Number of KL Terms", M);
282 solverParams.set("Mean", 0.5);
283 solverParams.set("Standard Deviation", 1.25);
284 Teuchos::Array<double> domain_upper(d);
285 Teuchos::Array<double> domain_lower(d);
286 Teuchos::Array<double>correlation_length(d);
287 for (int i=0; i<d; i++) {
288 domain_upper[i] = 1.5;
289 domain_lower[i] = -1.0;
290 correlation_length[i] = 10.0;
291 }
292 solverParams.set("Domain Upper Bounds", domain_upper);
293 solverParams.set("Domain Lower Bounds", domain_lower);
294 solverParams.set("Correlation Lengths", correlation_length);
295
296 // Since this test only runs on the host, instantiate the random field
297 // on the host
298 Stokhos::KL::ExponentialRandomField<double,Kokkos::DefaultHostExecutionSpace> rf(solverParams);
299
300 int nqp = quad_weights.size();
301 success = true;
302
303 // Evaluation point, scaled and shifted to proper domain
304 // Also map quadrature weights to right domain/density
305 Teuchos::Array<double> x(d);
306 Teuchos::Array<double> domain_center(d), domain_width(d);
307 double w_coeff = 1.0;
308 for (int i=0; i<d; i++) {
309 domain_center[i] = (domain_upper[i] + domain_lower[i])/2.0;
310 domain_width[i] = (domain_upper[i] - domain_lower[i])/2.0;
311 x[i] = domain_center[i] + 0.25*domain_width[i];
312 w_coeff *= 2.0*domain_width[i];
313 }
314
315 out << std::endl;
316 // Loop over each eigenpair (lambda, b(x))
317 for (int i=0; i<M; i++) {
318
319 // compute \int_D exp(-|x1-x1'|/L_1 - ... - |xd-xd'|/L_d)b(x') dx'
320 double integral = 0.0;
321 for (int qp=0; qp<nqp; qp++) {
322 Teuchos::Array<double> xp = quad_points[qp];
323 for (int j=0; j<d; j++)
324 xp[j] = domain_center[j] + xp[j]*domain_width[j];
325 double val = 0.0;
326 for (int j=0; j<d; j++)
327 val += std::abs(x[j] - xp[j])/correlation_length[j];
328 double w = w_coeff*quad_weights[qp];
329 integral +=
330 w*std::exp(-val)*rf.evaluate_eigenfunction(xp,i);
331 }
332
333 // compute lambda*b(x)
334 double rhs = rf.eigenvalue(i)*rf.evaluate_eigenfunction(x,i);
335 out << "lambda = " << rf.eigenvalue(i) << ", integral = " << integral
336 << ", rhs = " << rhs << ", error = " << integral-rhs << std::endl;
337 success = success &&
338 Stokhos::compareValues(integral, "integral", rhs, "rhs",
339 1e-3, 1e-3, out);
340 }
341 }
342
343}
j
j
Definition Sacado_Fad_Exp_MP_Vector.hpp:527
val
expr val()
Stokhos_KL_ExponentialRandomField.hpp
Stokhos_KL_OneDExponentialCovarianceFunction.hpp
Stokhos_UnitTestHelpers.hpp
Stokhos::CompletePolynomialBasis
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Definition Stokhos_CompletePolynomialBasis.hpp:74
Stokhos::KL::ExponentialRandomField
Class representing a KL expansion of an exponential random field.
Definition Stokhos_KL_ExponentialRandomField.hpp:108
Stokhos::KL::ExponentialRandomField::evaluate_eigenfunction
KOKKOS_INLINE_FUNCTION Teuchos::PromotionTraits< typenamepoint_type::value_type, value_type >::promote evaluate_eigenfunction(const point_type &point, int i) const
Evaluate given eigenfunction at a point.
Definition Stokhos_KL_ExponentialRandomFieldImp.hpp:203
Stokhos::KL::ExponentialRandomField::eigenvalue
value_type KOKKOS_INLINE_FUNCTION eigenvalue(int i) const
Return eigenvalue.
Definition Stokhos_KL_ExponentialRandomField.hpp:164
Stokhos::KL::OneDExponentialCovarianceFunction
Class representing an exponential covariance function and its KL eigevalues/eigenfunctions.
Definition Stokhos_KL_OneDExponentialCovarianceFunction.hpp:119
Stokhos::LegendreBasis
Legendre polynomial basis.
Definition Stokhos_LegendreBasis.hpp:68
Stokhos::RecurrenceBasis::getQuadPoints
virtual void getQuadPoints(ordinal_type quad_order, Teuchos::Array< value_type > &points, Teuchos::Array< value_type > &weights, Teuchos::Array< Teuchos::Array< value_type > > &values) const
Compute quadrature points, weights, and values of basis polynomials at given set of points points.
Definition Stokhos_RecurrenceBasisImp.hpp:344
Stokhos::TensorProductQuadrature
Defines quadrature for a tensor product basis by tensor products of 1-D quadrature rules.
Definition Stokhos_TensorProductQuadrature.hpp:58
Stokhos::TensorProductQuadrature::getQuadPoints
virtual const Teuchos::Array< Teuchos::Array< value_type > > & getQuadPoints() const
Get quadrature points.
Definition Stokhos_TensorProductQuadratureImp.hpp:164
Stokhos::TensorProductQuadrature::getQuadWeights
virtual const Teuchos::Array< value_type > & getQuadWeights() const
Get quadrature weights.
Definition Stokhos_TensorProductQuadratureImp.hpp:172
ExponentialRandomFieldUnitTest::TEUCHOS_UNIT_TEST
TEUCHOS_UNIT_TEST(Stokhos_ExponentialRandomField, OneD_Eigenvalue_Ordering)
Definition Stokhos_ExponentialRandomFieldUnitTest.cpp:54
Stokhos::compareValues
bool compareValues(const ValueType &a1, const std::string &a1_name, const ValueType &a2, const std::string &a2_name, const ValueType &rel_tol, const ValueType &abs_tol, Teuchos::FancyOStream &out)
Definition Stokhos_UnitTestHelpers.hpp:101
Generated by doxygen 1.17.0