Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_NormalizedHermiteBasisUnitTest.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#include "Teuchos_UnitTestHarness.hpp"
45#include "Teuchos_TestingHelpers.hpp"
46#include "Teuchos_UnitTestRepository.hpp"
47#include "Teuchos_GlobalMPISession.hpp"
48
49#include "Stokhos.hpp"
51#ifdef HAVE_STOKHOS_FORUQTK
52#include "Stokhos_gaussq.h"
53#endif
54
55namespace HermiteBasisUnitTest {
56
57 // Common setup for unit tests
58 template <typename OrdinalType, typename ValueType>
59 struct UnitTestSetup {
60 ValueType rtol, atol;
61 OrdinalType p;
63
64 UnitTestSetup() : rtol(1e-12), atol(1e-7), p(10), basis(p,true) {}
65
66 };
67
69
70 TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, Order ) {
71 int order = setup.basis.order();
72 TEUCHOS_TEST_EQUALITY(order, setup.p, out, success);
73 }
74
75 TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, Size ) {
76 int size = setup.basis.size();
77 TEUCHOS_TEST_EQUALITY(size, setup.p+1, out, success);
78 }
79
80 TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, Norm ) {
81 const Teuchos::Array<double>& n1 = setup.basis.norm_squared();
82 Teuchos::Array<double> n2(setup.p+1);
83 n2[0] = 1.0;
84 for (int i=1; i<=setup.p; i++)
85 n2[i] = 1.0;
86 success = Stokhos::compareArrays(n1, "n1", n2, "n2", setup.rtol, setup.atol,
87 out);
88 }
89
90 TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, Norm2 ) {
91 Teuchos::Array<double> n1(setup.p+1);
92 Teuchos::Array<double> n2(setup.p+1);
93 n1[0] = setup.basis.norm_squared(0);
94 n2[0] = 1.0;
95 for (int i=1; i<=setup.p; i++) {
96 n1[i] = setup.basis.norm_squared(i);
97 n2[i] = 1.0;
98 }
99 success = Stokhos::compareArrays(n1, "n1", n2, "n2", setup.rtol, setup.atol,
100 out);
101 }
102
103 TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, QuadNorm ) {
104 const Teuchos::Array<double>& n1 = setup.basis.norm_squared();
105 Teuchos::Array<double> n2(setup.p+1);
106 Teuchos::Array<double> x, w;
107 Teuchos::Array< Teuchos::Array<double> > v;
108 setup.basis.getQuadPoints(2*setup.p, x, w, v);
109 int nqp = w.size();
110 for (int i=0; i<=setup.p; i++) {
111 n2[i] = 0;
112 for (int j=0; j<nqp; j++)
113 n2[i] += w[j]*v[j][i]*v[j][i];
114 }
115 success = Stokhos::compareArrays(n1, "n1", n2, "n2", setup.rtol, setup.atol,
116 out);
117 }
118
119 // Tests basis is orthogonal
120 TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, QuadOrthog ) {
121 const Teuchos::Array<double>& norms = setup.basis.norm_squared();
122 Teuchos::Array<double> x, w;
123 Teuchos::Array< Teuchos::Array<double> > v;
124 setup.basis.getQuadPoints(2*setup.p, x, w, v);
125 int nqp = w.size();
126 Teuchos::SerialDenseMatrix<int,double> mat(setup.p+1, setup.p+1);
127 for (int i=0; i<=setup.p; i++) {
128 for (int j=0; j<=setup.p; j++) {
129 for (int k=0; k<nqp; k++)
130 mat(i,j) += w[k]*v[k][i]*v[k][j];
131 mat(i,j) /= std::sqrt(norms[i]*norms[j]);
132 }
133 mat(i,i) -= 1.0;
134 }
135 success = mat.normInf() < setup.atol;
136 if (!success) {
137 out << "\n Error, mat.normInf() < atol = " << mat.normInf()
138 << " < " << setup.atol << ": failed!\n";
139 out << "mat = " << printMat(mat) << std::endl;
140 }
141 }
142
143 TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, TripleProduct ) {
144 Teuchos::RCP< Stokhos::Dense3Tensor<int, double> > Cijk =
145 setup.basis.computeTripleProductTensor();
146
147 Teuchos::Array<double> x, w;
148 Teuchos::Array< Teuchos::Array<double> > v;
149 setup.basis.getQuadPoints(3*setup.p, x, w, v);
150
151 success = true;
152 for (int i=0; i<=setup.p; i++) {
153 for (int j=0; j<=setup.p; j++) {
154 for (int k=0; k<=setup.p; k++) {
155 double c = 0.0;
156 int nqp = w.size();
157 for (int qp=0; qp<nqp; qp++)
158 c += w[qp]*v[qp][i]*v[qp][j]*v[qp][k];
159 std::stringstream ss;
160 ss << "Cijk(" << i << "," << j << "," << k << ")";
161 success = success &&
162 Stokhos::compareValues((*Cijk)(i,j,k), ss.str(), c, "c",
163 setup.rtol, setup.atol, out);
164 }
165 }
166 }
167 }
168
169 TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, DerivDoubleProduct ) {
170 Teuchos::RCP< Teuchos::SerialDenseMatrix<int, double> > Bij =
171 setup.basis.computeDerivDoubleProductTensor();
172
173 Teuchos::Array<double> x, w;
174 Teuchos::Array< Teuchos::Array<double> > v, val, deriv;
175 setup.basis.getQuadPoints(2*setup.p, x, w, v);
176 int nqp = w.size();
177 val.resize(nqp);
178 deriv.resize(nqp);
179 for (int i=0; i<nqp; i++) {
180 val[i].resize(setup.p+1);
181 deriv[i].resize(setup.p+1);
182 setup.basis.evaluateBasesAndDerivatives(x[i], val[i], deriv[i]);
183 }
184
185 success = true;
186 for (int i=0; i<=setup.p; i++) {
187 for (int j=0; j<=setup.p; j++) {
188 double b = 0.0;
189 for (int qp=0; qp<nqp; qp++)
190 b += w[qp]*deriv[qp][i]*val[qp][j];
191 std::stringstream ss;
192 ss << "Bij(" << i << "," << j << ")";
193 success = success &&
194 Stokhos::compareValues((*Bij)(i,j), ss.str(), b, "b",
195 setup.rtol, setup.atol, out);
196 }
197 }
198 }
199
200 // TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, DerivDoubleProduct2 ) {
201 // Teuchos::RCP< Teuchos::SerialDenseMatrix<int, double> > Bij =
202 // setup.basis.computeDerivDoubleProductTensor();
203 // const Teuchos::Array<double>& n = setup.basis.norm_squared();
204
205 // success = true;
206 // for (int i=0; i<=setup.p; i++) {
207 // for (int j=0; j<=setup.p; j++) {
208 // double b = 0.0;
209 // if (j == i-1)
210 // b = i*n[j];
211 // std::stringstream ss;
212 // ss << "Bij(" << i << "," << j << ")";
213 // success = success &&
214 // Stokhos::compareValues((*Bij)(i,j), ss.str(), b, "b",
215 // setup.rtol, setup.atol, out);
216 // }
217 // }
218 // }
219
220 TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, EvaluateBases ) {
221 double x = 0.1234;
222 Teuchos::Array<double> v1(setup.p+1), v2(setup.p+1);
223 setup.basis.evaluateBases(x, v1);
224
225 // evaluate bases using formula
226 // He_0(x) = 1
227 // He_1(x) = x
228 // He_i(x) = (x*He_{i-1}(x) - sqrt(i-1)*He_{i-2}(x))/sqrt(i), i=2,3,...
229 v2[0] = 1.0;
230 if (setup.p >= 1)
231 v2[1] = x;
232 for (int i=2; i<=setup.p; i++)
233 v2[i] = (x*v2[i-1] - std::sqrt(i-1.0)*v2[i-2])/std::sqrt(1.0*i);
234 success = Stokhos::compareArrays(v1, "v1", v2, "v2", setup.rtol, setup.atol,
235 out);
236 }
237
238 TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, EvaluateBasesAndDerivatives ) {
239 double x = 0.1234;
240 Teuchos::Array<double> val1(setup.p+1), deriv1(setup.p+1),
241 val2(setup.p+1), deriv2(setup.p+1);
242 setup.basis.evaluateBasesAndDerivatives(x, val1, deriv1);
243
244 // evaluate bases and derivatives using formula:
245 // He_0(x) = 1
246 // He_1(x) = x
247 // He_i(x) = (x*He_{i-1}(x) - sqrt(i-1)*He_{i-2}(x))/sqrt(i), i=2,3,...
248 val2[0] = 1.0;
249 if (setup.p >= 1)
250 val2[1] = x;
251 for (int i=2; i<=setup.p; i++)
252 val2[i] = (x*val2[i-1] - std::sqrt(i-1.0)*val2[i-2])/std::sqrt(1.0*i);
253
254 deriv2[0] = 0.0;
255 if (setup.p >= 1)
256 deriv2[1] = 1.0;
257 for (int i=2; i<=setup.p; i++)
258 deriv2[i] = (val2[i-1] + x*deriv2[i-1] - std::sqrt(i-1.0)*deriv2[i-2]) /
259 std::sqrt(1.0*i);
260 success = Stokhos::compareArrays(val1, "val1", val2, "val2",
261 setup.rtol, setup.atol, out);
262 success = success &&
263 Stokhos::compareArrays(deriv1, "deriv1", deriv2, "deriv2",
264 setup.rtol, setup.atol, out);
265 }
266
267 TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, Evaluate ) {
268 double x = 0.1234;
269 Teuchos::Array<double> v1(setup.p+1), v2(setup.p+1);
270 for (int i=0; i<=setup.p; i++)
271 v1[i] = setup.basis.evaluate(x, i);
272
273 // evaluate bases using formula
274 // He_0(x) = 1
275 // He_1(x) = x
276 // He_i(x) = (x*He_{i-1}(x) - sqrt(i-1)*He_{i-2}(x))/sqrt(i), i=2,3,...
277 v2[0] = 1.0;
278 if (setup.p >= 1)
279 v2[1] = x;
280 for (int i=2; i<=setup.p; i++)
281 v2[i] = (x*v2[i-1] - std::sqrt(i-1.0)*v2[i-2])/std::sqrt(1.0*i);
282 success = Stokhos::compareArrays(v1, "v1", v2, "v2", setup.rtol, setup.atol,
283 out);
284 }
285
286 TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, Recurrence ) {
287 Teuchos::Array<double> a1(setup.p+1), b1(setup.p+1), c1(setup.p+1),
288 d1(setup.p+1);
289 Teuchos::Array<double> a2(setup.p+1), b2(setup.p+1), c2(setup.p+1),
290 d2(setup.p+1);
291 setup.basis.getRecurrenceCoefficients(a1, b1, c1, d1);
292
293 // compute coefficients using formula
294 a2[0] = 0.0; b2[0] = 1.0; c2[0] = 1.0; d2[0] = 1.0;
295 for (int i=1; i<=setup.p; i++) {
296 a2[i] = 0.0;
297 b2[i] = std::sqrt(1.0*i);
298 c2[i] = 1.0;
299 d2[i] = std::sqrt(1.0*i);
300 }
301 success = true;
302 success = success &&
303 Stokhos::compareArrays(a1, "a1", a2, "a2", setup.rtol, setup.atol, out);
304 success = success &&
305 Stokhos::compareArrays(b1, "b1", b2, "b2", setup.rtol, setup.atol, out);
306 success = success &&
307 Stokhos::compareArrays(c1, "c1", c2, "c2", setup.rtol, setup.atol, out);
308 success = success &&
309 Stokhos::compareArrays(d1, "d1", d2, "d2", setup.rtol, setup.atol, out);
310 }
311
312 // Tests alpha coefficients satisfy
313 // alpha_k = delta_k * (t*psi_k,psi_k)/(psi_k,psi_k)
314 TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, RecurrenceAlpha ) {
315 Teuchos::Array<double> a1(setup.p+1), b1(setup.p+1), c1(setup.p+1),
316 d1(setup.p+1);
317 setup.basis.getRecurrenceCoefficients(a1, b1, c1, d1);
318
319 Teuchos::Array<double> a2(setup.p+1);
320 const Teuchos::Array<double>& n1 = setup.basis.norm_squared();
321 Teuchos::Array<double> x, w;
322 Teuchos::Array< Teuchos::Array<double> > v;
323 setup.basis.getQuadPoints(2*setup.p, x, w, v);
324 int nqp = w.size();
325 for (int i=0; i<=setup.p; i++) {
326 a2[i] = 0;
327 for (int j=0; j<nqp; j++)
328 a2[i] += w[j]*x[j]*v[j][i]*v[j][i];
329 a2[i] = a2[i]*c1[i]/n1[i];
330 }
331 success = Stokhos::compareArrays(a1, "a1", a2, "a2", setup.rtol, setup.atol,
332 out);
333 }
334
335 // Tests beta coefficients satisfy
336 // beta_k =
337 // gamma_k * delta_k/delta_{k-1} * (psi_k,psi_k)/(psi_{k-1},psi_{k-1})
338 TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, RecurrenceBeta ) {
339 Teuchos::Array<double> a1(setup.p+1), b1(setup.p+1), c1(setup.p+1),
340 d1(setup.p+1);
341 setup.basis.getRecurrenceCoefficients(a1, b1, c1, d1);
342
343 Teuchos::Array<double> b2(setup.p+1);
344 const Teuchos::Array<double>& n1 = setup.basis.norm_squared();
345 b2[0] = b1[0];
346 for (int i=1; i<=setup.p; i++) {
347 b2[i] = d1[i]*(c1[i]/c1[i-1])*(n1[i]/n1[i-1]);
348 }
349 success = Stokhos::compareArrays(b1, "b1", b2, "b2", setup.rtol, setup.atol,
350 out);
351 }
352
353#ifdef HAVE_STOKHOS_DAKOTA
354 TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, QuadPointsDakota ) {
355 int n = static_cast<int>(std::ceil((setup.p+1)/2.0));
356 Teuchos::Array<double> x1, w1;
357 Teuchos::Array< Teuchos::Array<double> > v1;
358 setup.basis.getQuadPoints(setup.p, x1, w1, v1);
359
360 Teuchos::Array<double> x2(n), w2(n);
361 Teuchos::Array< Teuchos::Array<double> > v2(n);
362 webbur::hermite_compute(n, &x2[0], &w2[0]);
363
364 for (int i=0; i<n; i++) {
365 w2[i] *= 0.5/std::sqrt(std::atan(1.0)); // 1/sqrt(pi)
366 x2[i] *= std::sqrt(2.0);
367 v2[i].resize(setup.p+1);
368 setup.basis.evaluateBases(x2[i], v2[i]);
369 }
370 success = true;
371 success = success &&
372 Stokhos::compareArrays(x1, "x1", x2, "x2", setup.rtol, setup.atol, out);
373 success = success &&
374 Stokhos::compareArrays(w1, "w1", w2, "w2", setup.rtol, setup.atol, out);
375 for (int i=0; i<n; i++) {
376 std::stringstream ss1, ss2;
377 ss1 << "v1[" << i << "]";
378 ss2 << "v2[" << i << "]";
379 success = success &&
380 Stokhos::compareArrays(v1[i], ss1.str(), v2[i], ss2.str(),
381 setup.rtol, setup.atol, out);
382 }
383 }
384#endif
385
386#ifdef HAVE_STOKHOS_FORUQTK
387 TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, QuadPointsForUQTK ) {
388 int n = static_cast<int>(std::ceil((setup.p+1)/2.0));
389 Teuchos::Array<double> x1, w1;
390 Teuchos::Array< Teuchos::Array<double> > v1;
391 setup.basis.getQuadPoints(setup.p, x1, w1, v1);
392
393 Teuchos::Array<double> x2(n), w2(n);
394 Teuchos::Array< Teuchos::Array<double> > v2(n);
395 int kind = 4;
396 int kpts = 0;
397 double endpts[2] = {0.0, 0.0};
398 Teuchos::Array<double> b(n);
399 double alpha = 0.0;
400 double beta = 0.0;
401 GAUSSQ_F77(&kind, &n, &alpha, &beta, &kpts, endpts, &b[0], &x2[0], &w2[0]);
402
403 for (int i=0; i<n; i++) {
404 w2[i] *= 0.5/std::sqrt(std::atan(1.0)); // 1/sqrt(pi)
405 x2[i] *= std::sqrt(2.0);
406 v2[i].resize(setup.p+1);
407 setup.basis.evaluateBases(x2[i], v2[i]);
408 }
409 success = true;
410 success = success &&
411 Stokhos::compareArrays(x1, "x1", x2, "x2", setup.rtol, setup.atol, out);
412 success = success &&
413 Stokhos::compareArrays(w1, "w1", w2, "w2", setup.rtol, setup.atol, out);
414 for (int i=0; i<n; i++) {
415 std::stringstream ss1, ss2;
416 ss1 << "v1[" << i << "]";
417 ss2 << "v2[" << i << "]";
418 success = success &&
419 Stokhos::compareArrays(v1[i], ss1.str(), v2[i], ss2.str(),
420 setup.rtol, setup.atol, out);
421 }
422 }
423#endif
424
425}
426
427int main( int argc, char* argv[] ) {
428 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
429 return Teuchos::UnitTestRepository::runUnitTestsFromMain(argc, argv);
430}
expr val()
int main(int argc, char *argv[])
#define GAUSSQ_F77
Hermite polynomial basis.
UnitTestSetup< int, double > setup
TEUCHOS_UNIT_TEST(Stokhos_HermiteBasis, Order)
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)
bool compareArrays(const Array1 &a1, const std::string &a1_name, const Array2 &a2, const std::string &a2_name, const ValueType &rel_tol, const ValueType &abs_tol, Teuchos::FancyOStream &out)
Stokhos::HermiteBasis< OrdinalType, ValueType > basis
void printMat(const char *name, Epetra_IntSerialDenseMatrix &matrix)