Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
52 namespace 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  Stokhos::KL::ExponentialRandomField<double> rf(solverParams);
296 
297  int nqp = quad_weights.size();
298  success = true;
299 
300  // Evaluation point, scaled and shifted to proper domain
301  // Also map quadrature weights to right domain/density
302  Teuchos::Array<double> x(d);
303  Teuchos::Array<double> domain_center(d), domain_width(d);
304  double w_coeff = 1.0;
305  for (int i=0; i<d; i++) {
306  domain_center[i] = (domain_upper[i] + domain_lower[i])/2.0;
307  domain_width[i] = (domain_upper[i] - domain_lower[i])/2.0;
308  x[i] = domain_center[i] + 0.25*domain_width[i];
309  w_coeff *= 2.0*domain_width[i];
310  }
311 
312  out << std::endl;
313  // Loop over each eigenpair (lambda, b(x))
314  for (int i=0; i<M; i++) {
315 
316  // compute \int_D exp(-|x1-x1'|/L_1 - ... - |xd-xd'|/L_d)b(x') dx'
317  double integral = 0.0;
318  for (int qp=0; qp<nqp; qp++) {
319  Teuchos::Array<double> xp = quad_points[qp];
320  for (int j=0; j<d; j++)
321  xp[j] = domain_center[j] + xp[j]*domain_width[j];
322  double val = 0.0;
323  for (int j=0; j<d; j++)
324  val += std::abs(x[j] - xp[j])/correlation_length[j];
325  double w = w_coeff*quad_weights[qp];
326  integral +=
327  w*std::exp(-val)*rf.evaluate_eigenfunction(xp,i);
328  }
329 
330  // compute lambda*b(x)
331  double rhs = rf.eigenvalue(i)*rf.evaluate_eigenfunction(x,i);
332  out << "lambda = " << rf.eigenvalue(i) << ", integral = " << integral
333  << ", rhs = " << rhs << ", error = " << integral-rhs << std::endl;
334  success = success &&
335  Stokhos::compareValues(integral, "integral", rhs, "rhs",
336  1e-3, 1e-3, out);
337  }
338  }
339 
340 }
341 
342 int main( int argc, char* argv[] ) {
343  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
344  Kokkos::initialize();
345  int ret = Teuchos::UnitTestRepository::runUnitTestsFromMain(argc, argv);
346  Kokkos::finalize();
347  return ret;
348 }
argv
char * argv[]
Definition: Stokhos_HouseTriDiagUnitTest.cpp:286
Stokhos::TensorProductQuadrature::getQuadPoints
virtual const Teuchos::Array< Teuchos::Array< value_type > > & getQuadPoints() const
Get quadrature points.
Definition: Stokhos_TensorProductQuadratureImp.hpp:164
Teuchos::ParameterList::set
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
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::getQuadWeights
virtual const Teuchos::Array< value_type > & getQuadWeights() const
Get quadrature weights.
Definition: Stokhos_TensorProductQuadratureImp.hpp:172
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
Stokhos_KL_OneDExponentialCovarianceFunction.hpp
j
j
Definition: Sacado_Fad_Exp_MP_Vector.hpp:527
Teuchos_UnitTestRepository.hpp
Teuchos::UnitTestRepository::runUnitTestsFromMain
static int runUnitTestsFromMain(int argc, char *argv[])
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Stokhos::KL::ExponentialRandomField::evaluate_eigenfunction
KOKKOS_INLINE_FUNCTION Teuchos::PromotionTraits< typename point_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
Teuchos_UnitTestHarness.hpp
Sacado::UQ::abs
KOKKOS_INLINE_FUNCTION PCE< Storage > abs(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:1181
Stokhos::KL::OneDExponentialCovarianceFunction
Class representing an exponential covariance function and its KL eigevalues/eigenfunctions.
Definition: Stokhos_KL_OneDExponentialCovarianceFunction.hpp:119
Stokhos_KL_ExponentialRandomField.hpp
Stokhos_UnitTestHelpers.hpp
Sacado::UQ::atan
KOKKOS_INLINE_FUNCTION PCE< Storage > atan(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:1101
Sacado::UQ::exp
KOKKOS_INLINE_FUNCTION PCE< Storage > exp(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:827
Stokhos::LegendreBasis
Legendre polynomial basis.
Definition: Stokhos_LegendreBasis.hpp:67
main
int main(int argc, char **argv)
Definition: cijk_ltb_partition.cpp:57
Teuchos_TestingHelpers.hpp
val
expr val()
Teuchos_GlobalMPISession.hpp
ExponentialRandomFieldUnitTest::TEUCHOS_UNIT_TEST
TEUCHOS_UNIT_TEST(Stokhos_ExponentialRandomField, OneD_Eigenvalue_Ordering)
Definition: Stokhos_ExponentialRandomFieldUnitTest.cpp:54
Teuchos::Array::size
size_type size() const
Stokhos::KL::ExponentialRandomField
Class representing a KL expansion of an exponential random field.
Definition: Stokhos_KL_ExponentialRandomField.hpp:109
Stokhos::KL::ExponentialRandomField::eigenvalue
value_type KOKKOS_INLINE_FUNCTION eigenvalue(int i) const
Return eigenvalue.
Definition: Stokhos_KL_ExponentialRandomField.hpp:165
Stokhos::TensorProductQuadrature
Defines quadrature for a tensor product basis by tensor products of 1-D quadrature rules...
Definition: Stokhos_TensorProductQuadrature.hpp:58
Generated by   doxygen 1.8.5