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