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 // Stokhos Package
4 //
5 // Copyright 2009 NTESS and the Stokhos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #include "Teuchos_UnitTestHarness.hpp"
11 #include "Teuchos_TestingHelpers.hpp"
12 #include "Teuchos_UnitTestRepository.hpp"
13 #include "Teuchos_GlobalMPISession.hpp"
14 
15 #include "Stokhos.hpp"
16 #include "Stokhos_KL_OneDExponentialCovarianceFunction.hpp"
17 #include "Stokhos_KL_ExponentialRandomField.hpp"
18 #include "Stokhos_UnitTestHelpers.hpp"
19 
20 namespace ExponentialRandomFieldUnitTest {
21 
22  TEUCHOS_UNIT_TEST( Stokhos_ExponentialRandomField, OneD_Eigenvalue_Ordering ) {
23  int M = 100;
24  double a = -1.0;
25  double b = 1.5;
26  double L = 1.1;
27 
28  // Setup covariance function
29  Teuchos::ParameterList solverParams;
30  solverParams.set("Nonlinear Solver Tolerance", 1e-8);
31  typedef Stokhos::KL::OneDExponentialCovarianceFunction<double> cov_type;
32  typedef cov_type::eigen_pair_type eigen_pair_type;
33  cov_type cov(M, a, b, L, 0, solverParams);
34 
35  // Get eigenpairs
36  const Teuchos::Array<eigen_pair_type>& eigs = cov.getEigenPairs();
37 
38  success = true;
39  out << std::endl;
40  for (int i=0; i<M-1; i++) {
41  out << "eigs[" << i << "] = " << eigs[i].eig_val << std::endl;
42  if (eigs[i].eig_val < eigs[i+1].eig_val)
43  success = false;
44  }
45  }
46 
47  TEUCHOS_UNIT_TEST( Stokhos_ExponentialRandomField, OneD_Frequency_Spacing ) {
48  int M = 100;
49  double a = -1.0;
50  double b = 1.5;
51  double L = 1.1;
52 
53  // Setup covariance function
54  Teuchos::ParameterList solverParams;
55  solverParams.set("Nonlinear Solver Tolerance", 1e-8);
56  typedef Stokhos::KL::OneDExponentialCovarianceFunction<double> cov_type;
57  typedef cov_type::eigen_pair_type eigen_pair_type;
58  cov_type cov(M, a, b, L, 0, solverParams);
59 
60  // Get eigenpairs
61  const Teuchos::Array<eigen_pair_type>& eigs = cov.getEigenPairs();
62 
63  success = true;
64  out << std::endl;
65  double pi = 4.0*std::atan(1.0);
66  for (int i=0; i<M-1; i++) {
67  double omega1 = eigs[i].eig_func.getFrequency();
68  double omega2 = eigs[i].eig_func.getFrequency();
69 
70  out << "eigs[" << i << "].frequency = " << omega1 << std::endl;
71  if (omega2 - omega1 > pi)
72  success = false;
73  }
74  }
75 
76  TEUCHOS_UNIT_TEST( Stokhos_ExponentialRandomField, OneD_Eigenfunction_Norm ) {
77  int p = 200;
78  int M = 10;
79  double a = -1.0;
80  double b = 1.5;
81  double L = 1.1;
82  double center = (b+a)/2.0;
83  double width = (b-a)/2.0;
84  double w_coeff = 2.0*width;
85 
86  // Create basis for getting quadrature points
87  Stokhos::LegendreBasis<int,double> basis(p);
88  Teuchos::Array<double> quad_points, quad_weights;
89  Teuchos::Array< Teuchos::Array<double> > quad_values;
90  basis.getQuadPoints(p, quad_points, quad_weights, quad_values);
91 
92  // Setup covariance function
93  Teuchos::ParameterList solverParams;
94  typedef Stokhos::KL::OneDExponentialCovarianceFunction<double> cov_type;
95  typedef cov_type::eigen_pair_type eigen_pair_type;
96  cov_type cov(M, a, b, L, 0, solverParams);
97 
98  // Get eigenpairs
99  const Teuchos::Array<eigen_pair_type>& eigs = cov.getEigenPairs();
100 
101  int nqp = quad_weights.size();
102  success = true;
103 
104  out << std::endl;
105  // Loop over each eigenpair (lambda, b(x))
106  for (int i=0; i<M; i++) {
107 
108  // compute \int_D b(x)*b(x) dx
109  double integral = 0.0;
110  double rhs = 1.0;
111  for (int qp=0; qp<nqp; qp++) {
112  double xp = center + quad_points[qp]*width;
113  double w = w_coeff*quad_weights[qp];
114  double val = eigs[i].eig_func.evaluate(xp);
115  integral += w*val*val;
116  }
117 
118  out << "lambda = " << eigs[i].eig_val << ", integral = " << integral
119  << ", rhs = " << rhs << ", error = " << integral-rhs << std::endl;
120  success = success &&
121  Stokhos::compareValues(integral, "integral", rhs, "rhs",
122  1e-3, 1e-3, out);
123  }
124  }
125 
126  TEUCHOS_UNIT_TEST( Stokhos_ExponentialRandomField, OneD_Eigenfunction_Orthogonality ) {
127  int p = 200;
128  int M = 10;
129  double a = -1.0;
130  double b = 1.5;
131  double L = 1.1;
132  double center = (b+a)/2.0;
133  double width = (b-a)/2.0;
134  double w_coeff = 2.0*width;
135 
136  // Create basis for getting quadrature points
137  Stokhos::LegendreBasis<int,double> basis(p);
138  Teuchos::Array<double> quad_points, quad_weights;
139  Teuchos::Array< Teuchos::Array<double> > quad_values;
140  basis.getQuadPoints(p, quad_points, quad_weights, quad_values);
141 
142  // Setup covariance function
143  Teuchos::ParameterList solverParams;
144  typedef Stokhos::KL::OneDExponentialCovarianceFunction<double> cov_type;
145  typedef cov_type::eigen_pair_type eigen_pair_type;
146  cov_type cov(M, a, b, L, 0, solverParams);
147 
148  // Get eigenpairs
149  const Teuchos::Array<eigen_pair_type>& eigs = cov.getEigenPairs();
150 
151  int nqp = quad_weights.size();
152  success = true;
153 
154  out << std::endl;
155  for (int i=0; i<M; i++) {
156  for (int j=0; j<i; j++) {
157 
158  // compute \int_D b_i(x)*b_j(x) dx
159  double integral = 0.0;
160  double rhs = 0.0;
161  for (int qp=0; qp<nqp; qp++) {
162  double xp = center + quad_points[qp]*width;
163  double w = w_coeff*quad_weights[qp];
164  double val1 = eigs[i].eig_func.evaluate(xp);
165  double val2 = eigs[j].eig_func.evaluate(xp);
166  integral += w*val1*val2;
167  }
168 
169  out << "lambda = " << eigs[i].eig_val << ", integral = " << integral
170  << ", rhs = " << rhs << ", error = " << integral-rhs << std::endl;
171  success = success &&
172  Stokhos::compareValues(integral, "integral", rhs, "rhs",
173  1e-3, 1e-3, out);
174  }
175  }
176  }
177 
178  TEUCHOS_UNIT_TEST( Stokhos_ExponentialRandomField, OneD_Eigen_Solution ) {
179  int p = 200;
180  int M = 10;
181  double a = -1.0;
182  double b = 1.5;
183  double L = 1.1;
184  double center = (b+a)/2.0;
185  double width = (b-a)/2.0;
186  double x = center + 0.25*width;
187  double w_coeff = 2.0*width;
188 
189  // Create basis for getting quadrature points
190  Stokhos::LegendreBasis<int,double> basis(p);
191  Teuchos::Array<double> quad_points, quad_weights;
192  Teuchos::Array< Teuchos::Array<double> > quad_values;
193  basis.getQuadPoints(p, quad_points, quad_weights, quad_values);
194 
195  // Setup covariance function
196  Teuchos::ParameterList solverParams;
197  typedef Stokhos::KL::OneDExponentialCovarianceFunction<double> cov_type;
198  typedef cov_type::eigen_pair_type eigen_pair_type;
199  cov_type cov(M, a, b, L, 0, solverParams);
200 
201  // Get eigenpairs
202  const Teuchos::Array<eigen_pair_type>& eigs = cov.getEigenPairs();
203 
204  int nqp = quad_weights.size();
205  success = true;
206 
207  out << std::endl;
208  // Loop over each eigenpair (lambda, b(x))
209  for (int i=0; i<M; i++) {
210 
211  // compute \int_D exp(-|x-x'|/L)b(x') dx'
212  double integral = 0.0;
213  for (int qp=0; qp<nqp; qp++) {
214  double xp = center + quad_points[qp]*width;
215  double w = w_coeff*quad_weights[qp];
216  integral +=
217  w*cov.evaluateCovariance(x,xp)*eigs[i].eig_func.evaluate(xp);
218  }
219 
220  // compute lambda*b(x)
221  double rhs = eigs[i].eig_val*eigs[i].eig_func.evaluate(x);
222  out << "lambda = " << eigs[i].eig_val << ", integral = " << integral
223  << ", rhs = " << rhs << ", error = " << integral-rhs << std::endl;
224  success = success &&
225  Stokhos::compareValues(integral, "integral", rhs, "rhs",
226  1e-3, 1e-3, out);
227  }
228  }
229 
230  TEUCHOS_UNIT_TEST( Stokhos_ExponentialRandomField, Product_Eigen_Solution ) {
231  // Create product basis
232  int p = 20;
233  int d = 2;
234  int M = 10;
235  Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
236  for (int i=0; i<d; i++)
237  bases[i] = Teuchos::rcp(new Stokhos::LegendreBasis<int,double>(p));
238  Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
239  Teuchos::rcp(new Stokhos::CompletePolynomialBasis<int,double>(bases));
240 
241  // Tensor product quadrature
242  Stokhos::TensorProductQuadrature<int,double> quad(basis);
243  const Teuchos::Array< Teuchos::Array<double> >& quad_points =
244  quad.getQuadPoints();
245  const Teuchos::Array<double>& quad_weights = quad.getQuadWeights();
246 
247  // Setup random field
248  Teuchos::ParameterList solverParams;
249  solverParams.set("Number of KL Terms", M);
250  solverParams.set("Mean", 0.5);
251  solverParams.set("Standard Deviation", 1.25);
252  Teuchos::Array<double> domain_upper(d);
253  Teuchos::Array<double> domain_lower(d);
254  Teuchos::Array<double>correlation_length(d);
255  for (int i=0; i<d; i++) {
256  domain_upper[i] = 1.5;
257  domain_lower[i] = -1.0;
258  correlation_length[i] = 10.0;
259  }
260  solverParams.set("Domain Upper Bounds", domain_upper);
261  solverParams.set("Domain Lower Bounds", domain_lower);
262  solverParams.set("Correlation Lengths", correlation_length);
263 
264  // Since this test only runs on the host, instantiate the random field
265  // on the host
266  Stokhos::KL::ExponentialRandomField<double,Kokkos::DefaultHostExecutionSpace> rf(solverParams);
267 
268  int nqp = quad_weights.size();
269  success = true;
270 
271  // Evaluation point, scaled and shifted to proper domain
272  // Also map quadrature weights to right domain/density
273  Teuchos::Array<double> x(d);
274  Teuchos::Array<double> domain_center(d), domain_width(d);
275  double w_coeff = 1.0;
276  for (int i=0; i<d; i++) {
277  domain_center[i] = (domain_upper[i] + domain_lower[i])/2.0;
278  domain_width[i] = (domain_upper[i] - domain_lower[i])/2.0;
279  x[i] = domain_center[i] + 0.25*domain_width[i];
280  w_coeff *= 2.0*domain_width[i];
281  }
282 
283  out << std::endl;
284  // Loop over each eigenpair (lambda, b(x))
285  for (int i=0; i<M; i++) {
286 
287  // compute \int_D exp(-|x1-x1'|/L_1 - ... - |xd-xd'|/L_d)b(x') dx'
288  double integral = 0.0;
289  for (int qp=0; qp<nqp; qp++) {
290  Teuchos::Array<double> xp = quad_points[qp];
291  for (int j=0; j<d; j++)
292  xp[j] = domain_center[j] + xp[j]*domain_width[j];
293  double val = 0.0;
294  for (int j=0; j<d; j++)
295  val += std::abs(x[j] - xp[j])/correlation_length[j];
296  double w = w_coeff*quad_weights[qp];
297  integral +=
298  w*std::exp(-val)*rf.evaluate_eigenfunction(xp,i);
299  }
300 
301  // compute lambda*b(x)
302  double rhs = rf.eigenvalue(i)*rf.evaluate_eigenfunction(x,i);
303  out << "lambda = " << rf.eigenvalue(i) << ", integral = " << integral
304  << ", rhs = " << rhs << ", error = " << integral-rhs << std::endl;
305  success = success &&
306  Stokhos::compareValues(integral, "integral", rhs, "rhs",
307  1e-3, 1e-3, out);
308  }
309  }
310 
311 }
Stokhos::TensorProductQuadrature::getQuadPoints
virtual const Teuchos::Array< Teuchos::Array< value_type > > & getQuadPoints() const
Get quadrature points.
Definition: Stokhos_TensorProductQuadratureImp.hpp:131
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:310
Stokhos::TensorProductQuadrature::getQuadWeights
virtual const Teuchos::Array< value_type > & getQuadWeights() const
Get quadrature weights.
Definition: Stokhos_TensorProductQuadratureImp.hpp:139
Teuchos::ParameterList::set
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
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:69
Stokhos_KL_OneDExponentialCovarianceFunction.hpp
j
j
Definition: Sacado_Fad_Exp_MP_Vector.hpp:495
Teuchos_UnitTestRepository.hpp
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:171
Teuchos_UnitTestHarness.hpp
Sacado::UQ::abs
KOKKOS_INLINE_FUNCTION PCE< Storage > abs(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:1150
Stokhos::KL::OneDExponentialCovarianceFunction
Class representing an exponential covariance function and its KL eigevalues/eigenfunctions.
Definition: Stokhos_KL_OneDExponentialCovarianceFunction.hpp:87
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:1070
Sacado::UQ::exp
KOKKOS_INLINE_FUNCTION PCE< Storage > exp(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:796
Stokhos::LegendreBasis
Legendre polynomial basis.
Definition: Stokhos_LegendreBasis.hpp:33
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:22
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:76
Stokhos::KL::ExponentialRandomField::eigenvalue
value_type KOKKOS_INLINE_FUNCTION eigenvalue(int i) const
Return eigenvalue.
Definition: Stokhos_KL_ExponentialRandomField.hpp:132
Stokhos::TensorProductQuadrature
Defines quadrature for a tensor product basis by tensor products of 1-D quadrature rules...
Definition: Stokhos_TensorProductQuadrature.hpp:24
Generated on Fri Nov 22 2024 09:18:11 for Stokhos Package Browser (Single Doxygen Collection) by   doxygen 1.8.5