Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
stieltjes_coupled.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 <iostream>
45#include <iomanip>
46
47#include "Stokhos.hpp"
50
51#include "Teuchos_SerialDenseMatrix.hpp"
52#include "Teuchos_SerialDenseVector.hpp"
53#include "Teuchos_LAPACK.hpp"
54
55#include <fstream>
56#include <iostream>
57
59
60double rel_error(double a, double b) {
61 return std::abs(a-b)/std::max(std::abs(a), std::abs(b));
62}
63
64class RhoModel {
65public:
66 RhoModel(int n, double dx) : A(n,n), b(n), h(dx), piv(n) {}
67 void computeRho(const Teuchos::Array<double>& s2,
68 const Teuchos::Array<double>& gamma,
69 Teuchos::Array<double>& rho) {
70 // Loop over quadrature points
71 int n = b.length();
72 int m = s2.size();
73 for (int qp=0; qp<m; qp++) {
74
75 // Fill matrix
76 A.putScalar(0.0);
77 A(0,0) = -2.0/(h*h);
78 A(0,1) = s2[qp]/(2.*h) + 1.0/(h*h);
79 for (int i=1; i<n-1; i++) {
80 A(i,i) = -2.0/(h*h);
81 A(i,i-1) = -s2[qp]/(2.*h) + 1.0/(h*h);
82 A(i,i+1) = s2[qp]/(2.*h) + 1.0/(h*h);
83 }
84 A(n-1,n-2) = -s2[qp]/(2.*h) + 1.0/(h*h);
85 A(n-1,n-1) = -2.0/(h*h);
86
87 // Fill rhs
88 b.putScalar(gamma[qp]);
89
90 // Solve system
91 int info;
92 lapack.GESV(n, 1, A.values(), n, &piv[0], b.values(), n, &info);
93
94 // Compute rho
95 rho[qp] = 0.0;
96 for (int i=0; i<n; i++)
97 rho[qp] += b(i);
98 rho[qp] *= h;
99 }
100 }
101
108
109 // Get quadrature data
110 const Teuchos::Array<double>& weights =
111 quad.getQuadWeights();
112 const Teuchos::Array<Teuchos::Array<double> >& points =
113 quad.getQuadPoints();
114 const Teuchos::Array<Teuchos::Array<double> >& basis_vals =
116 const Teuchos::Array<double>& norms = basis.norm_squared();
117 int nqp = weights.size();
118 int sz = basis.size();
119
120 // Evaluate input PCEs at quad points
121 Teuchos::Array<double> s2(nqp), r(nqp), g(nqp);
122 for (int qp=0; qp<nqp; qp++) {
123 s2[qp] = xi2.evaluate(points[qp], basis_vals[qp]);
124 g[qp] = gamma.evaluate(points[qp], basis_vals[qp]);
125 }
126
127 // Compute rho at quad points
128 computeRho(s2, g, r);
129
130 // Compute rho pce
131 rho.init(0.0);
132 for (int qp=0; qp<nqp; qp++) {
133 for (int i=0; i<sz; i++)
134 rho[i] += r[qp]*weights[qp]*basis_vals[qp][i]/norms[i];
135 }
136 }
137
138private:
139 Teuchos::SerialDenseMatrix<int,double> A;
140 Teuchos::SerialDenseVector<int,double> b;
141 double h;
143 Teuchos::Array<int> piv;
144};
145
147public:
148 GammaModel(int n, double dx) : A(n,n), b(n), h(dx), piv(n) {}
149 void computeGamma(const Teuchos::Array<double>& s1,
150 const Teuchos::Array<double>& rho,
151 Teuchos::Array<double>& gamma) {
152 // Loop over quadrature points
153 int n = b.length();
154 int m = s1.size();
155 for (int qp=0; qp<m; qp++) {
156
157 // Fill matrix
158 A.putScalar(0.0);
159 A(0,0) = -s1[qp]*2.0/(h*h) + rho[qp];
160 A(0,1) = s1[qp]/(h*h);
161 for (int i=1; i<n-1; i++) {
162 A(i,i) = -s1[qp]*2.0/(h*h) + rho[qp];
163 A(i,i-1) = s1[qp]/(h*h);
164 A(i,i+1) = s1[qp]/(h*h);
165 }
166 A(n-1,n-2) = s1[qp]/(h*h);
167 A(n-1,n-1) = -s1[qp]*2.0/(h*h) + rho[qp];
168
169 // Fill rhs
170 double pi = 4.0*std::atan(1.0);
171 for (int i=0; i<n; i++) {
172 double x = h + i*h;
173 b(i) = std::sin(2.0*pi*x);
174 }
175
176 // Solve system
177 int info;
178 lapack.GESV(n, 1, A.values(), n, &piv[0], b.values(), n, &info);
179
180 // Compute gamma
181 gamma[qp] = 0.0;
182 for (int i=0; i<n; i++)
183 gamma[qp] += b(i)*b(i);
184 gamma[qp] *= 100.0*h;
185 }
186 }
187
194
195 // Get quadrature data
196 const Teuchos::Array<double>& weights =
197 quad.getQuadWeights();
198 const Teuchos::Array<Teuchos::Array<double> >& points =
199 quad.getQuadPoints();
200 const Teuchos::Array<Teuchos::Array<double> >& basis_vals =
202 const Teuchos::Array<double>& norms = basis.norm_squared();
203 int nqp = weights.size();
204 int sz = basis.size();
205
206 // Evaluate input PCEs at quad points
207 Teuchos::Array<double> s1(nqp), r(nqp), g(nqp);
208 for (int qp=0; qp<nqp; qp++) {
209 s1[qp] = xi1.evaluate(points[qp], basis_vals[qp]);
210 r[qp] = rho.evaluate(points[qp], basis_vals[qp]);
211 }
212
213 // Compute gamma at quad points
214 computeGamma(s1, r, g);
215
216 // Compute gamma pce
217 gamma.init(0.0);
218 for (int qp=0; qp<nqp; qp++) {
219 for (int i=0; i<sz; i++)
220 gamma[i] += g[qp]*weights[qp]*basis_vals[qp][i]/norms[i];
221 }
222 }
223private:
224 Teuchos::SerialDenseMatrix<int,double> A;
225 Teuchos::SerialDenseVector<int,double> b;
226 double h;
228 Teuchos::Array<int> piv;
229};
230
232public:
234 double tol_, int max_it_,
235 RhoModel& rhoModel_, GammaModel& gammaModel_) :
236 tol(tol_), max_it(max_it_), rhoModel(rhoModel_), gammaModel(gammaModel_) {}
237
238 void solve(double xi1, double xi2, double& rho, double& gamma) {
239 double rho0 = 0.0;
240 double gamma0 = 0.0;
241 double rho_error = 1.0;
242 double gamma_error = 1.0;
243 int it = 0;
244 Teuchos::Array<double> s1(1), s2(1), r(1), g(1);
245 s1[0] = xi1;
246 s2[0] = xi2;
247 r[0] = 1.0;
248 g[0] = 1.0;
249 while( (rho_error > tol || gamma_error > tol) && it < max_it) {
250
251 // Compute rho
252 rhoModel.computeRho(s2, g, r);
253
254 // Compute gamma
255 gammaModel.computeGamma(s1, r, g);
256
257 // Compute errors
258 rho_error = rel_error(r[0], rho0);
259 gamma_error = rel_error(g[0],gamma0);
260
261 // Update
262 rho0 = r[0];
263 gamma0 = g[0];
264 ++it;
265 }
266
267 // Check if we converged
268 if (it >= max_it)
269 throw std::string("model didn't converge!");
270
271 rho = r[0];
272 gamma = g[0];
273 }
274private:
275 double tol;
279 Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> > basis;
280 Teuchos::RCP<const Stokhos::Quadrature<int,double> > quad;
281};
282
284public:
286 double tol, int max_it,
287 RhoModel& rhoModel, GammaModel& gammaModel,
288 const Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> >& basis_,
289 const Teuchos::RCP<const Stokhos::Quadrature<int,double> >& quad_) :
290 coupledSolver(tol, max_it, rhoModel, gammaModel),
291 basis(basis_), quad(quad_) {}
292
297
298 rho_pce.init(0.0);
299 gamma_pce.init(0.0);
300
301 // Get quadrature data
302 const Teuchos::Array<double>& weights =
303 quad->getQuadWeights();
304 const Teuchos::Array<Teuchos::Array<double> >& points =
305 quad->getQuadPoints();
306 const Teuchos::Array<Teuchos::Array<double> >& basis_vals =
307 quad->getBasisAtQuadPoints();
308 const Teuchos::Array<double>& norms = basis->norm_squared();
309
310 // Solve coupled system at each quadrature point
311 double rho, gamma;
312 int nqp = weights.size();
313 for (int qp=0; qp<nqp; qp++) {
314 double s1 = xi1.evaluate(points[qp], basis_vals[qp]);
315 double s2 = xi2.evaluate(points[qp], basis_vals[qp]);
316
317 coupledSolver.solve(s1, s2, rho, gamma);
318
319 // Sum rho and gamma into their respective PCEs
320 int sz = basis->size();
321 for (int i=0; i<sz; i++) {
322 rho_pce[i] += rho*weights[qp]*basis_vals[qp][i]/norms[i];
323 gamma_pce[i] += gamma*weights[qp]*basis_vals[qp][i]/norms[i];
324 }
325 }
326 }
327private:
329 Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> > basis;
330 Teuchos::RCP<const Stokhos::Quadrature<int,double> > quad;
331};
332
334public:
336 double tol_, int max_it_,
337 RhoModel& rhoModel_, GammaModel& gammaModel_,
338 const Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> >& basis_,
339 const Teuchos::RCP<const Stokhos::Quadrature<int,double> >& quad_) :
340 tol(tol_), max_it(max_it_), rhoModel(rhoModel_), gammaModel(gammaModel_),
341 basis(basis_), quad(quad_) {}
342
347
350 double rho_error = 1.0;
351 double gamma_error = 1.0;
352 int it = 0;
353 int sz = basis->size();
354
355 while( (rho_error > tol || gamma_error > tol) && it < max_it) {
356
357 // Compute rho pce
358 rhoModel.computeRhoPCE(*basis, *quad, xi2, gamma, rho);
359
360 // Compute gamma pce
361 gammaModel.computeGammaPCE(*basis, *quad, xi1, rho, gamma);
362
363 // Compute errors
364 rho_error = 0.0;
365 gamma_error = 0.0;
366 for (int i=0; i<sz; i++) {
367 double re = rel_error(rho[i],rho0[i]);
368 double ge = rel_error(gamma[i],gamma0[i]);
369 if (re > rho_error)
370 rho_error = re;
371 if (ge > gamma_error)
372 gamma_error = ge;
373 }
374
375 std::cout << "it = " << it
376 << " rho_error = " << rho_error
377 << " gamma_error = " << gamma_error << std::endl;
378
379 // Update
380 rho0 = rho;
381 gamma0 = gamma;
382 ++it;
383 }
384
385 // Check if we converged
386 if (it >= max_it)
387 throw std::string("model didn't converge!");
388 }
389private:
390 double tol;
394 Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> > basis;
395 Teuchos::RCP<const Stokhos::Quadrature<int,double> > quad;
396};
397
399public:
401 double tol_, int max_it_,
402 RhoModel& rhoModel_, GammaModel& gammaModel_,
403 const Teuchos::RCP<const Stokhos::ProductBasis<int,double> >& basis_,
404 const Teuchos::RCP<const Stokhos::Quadrature<int,double> >& quad_) :
405 tol(tol_), max_it(max_it_), rhoModel(rhoModel_), gammaModel(gammaModel_),
406 basis(basis_), quad(quad_) {}
407
412
413 // Get quadrature data
414 const Teuchos::Array<double>& weights =
415 quad->getQuadWeights();
416 const Teuchos::Array<Teuchos::Array<double> >& points =
417 quad->getQuadPoints();
418 const Teuchos::Array<Teuchos::Array<double> >& basis_vals =
419 quad->getBasisAtQuadPoints();
420 const Teuchos::Array<double>& norms = basis->norm_squared();
421
424 double rho_error = 1.0;
425 double gamma_error = 1.0;
426 int it = 0;
427 int nqp = weights.size();
428 int sz = basis->size();
429 int p = basis->order();
430 bool use_pce_quad_points = false;
431 bool normalize = false;
432 bool project_integrals = false;
433
434 // Triple product tensor
435 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk;
436 if (project_integrals)
437 Cijk = basis->computeTripleProductTensor();
438
439 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > >
440 coordinate_bases = basis->getCoordinateBases();
441 Teuchos::RCP<const Stokhos::Quadrature<int,double> > st_quad = quad;
442 // Teuchos::rcp(new Stokhos::TensorProductQuadrature<int,double>(basis,
443 // 2*(p+1)));
444
445 while( (rho_error > tol || gamma_error > tol) && it < max_it) {
446
447 // Compute basis for rho
448 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > rho_bases(2);
449 rho_bases[0] = coordinate_bases[1];
450 rho_bases[1] =
452 p, Teuchos::rcp(&gamma,false), st_quad,
453 use_pce_quad_points,
454 normalize, project_integrals, Cijk));
455 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> >
456 rho_basis =
457 Teuchos::rcp(new Stokhos::CompletePolynomialBasis<int,double>(rho_bases));
458
459 // Compute quadrature for rho
460 Teuchos::RCP<const Stokhos::Quadrature<int,double> > rho_quad =
461 Teuchos::rcp(new Stokhos::TensorProductQuadrature<int,double>(rho_basis));
462
463 // Write xi2 and gamma in rho basis
464 Stokhos::OrthogPolyApprox<int,double> xi2_rho(rho_basis),
465 gamma_rho(rho_basis), rho_rho(rho_basis);
466 xi2_rho.term(0, 0) = xi2.term(1,0); // this assumes linear expansion
467 xi2_rho.term(0, 1) = xi2.term(1,1);
468 gamma_rho.term(1, 0) = gamma.mean();
469 //gamma_rho.term(1, 1) = gamma.standard_deviation();
470 gamma_rho.term(1, 1) = 1.0;
471
472 // Compute rho pce in transformed basis
473 rhoModel.computeRhoPCE(*rho_basis, *rho_quad, xi2_rho, gamma_rho,
474 rho_rho);
475
476 // if (it == 0)
477 // std::cout << rho_rho << std::endl;
478
479 // Project rho back to original basis
480 Teuchos::Array<double> rho_point(2);
481 Teuchos::Array<double> rho_basis_vals(rho_basis->size());
482 rho.init(0.0);
483 for (int qp=0; qp<nqp; qp++) {
484 rho_point[0] = points[qp][1];
485 rho_point[1] = gamma.evaluate(points[qp], basis_vals[qp]);
486 rho_basis->evaluateBases(rho_point, rho_basis_vals);
487 double r = rho_rho.evaluate(rho_point ,rho_basis_vals);
488 for (int i=0; i<sz; i++)
489 rho[i] += r*weights[qp]*basis_vals[qp][i]/norms[i];
490 }
491
492 // Compute basis for gamma
493 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > gamma_bases(2);
494 gamma_bases[0] = coordinate_bases[0];
495 gamma_bases[1] =
497 p, Teuchos::rcp(&rho,false), st_quad,
498 use_pce_quad_points,
499 normalize, project_integrals, Cijk));
500 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> >
501 gamma_basis =
502 Teuchos::rcp(new Stokhos::CompletePolynomialBasis<int,double>(gamma_bases));
503
504 // Compute quadrature for gamma
505 Teuchos::RCP<const Stokhos::Quadrature<int,double> > gamma_quad =
506 Teuchos::rcp(new Stokhos::TensorProductQuadrature<int,double>(gamma_basis));
507
508 // Write xi1 and rho in gamma basis
509 Stokhos::OrthogPolyApprox<int,double> xi1_gamma(gamma_basis),
510 gamma_gamma(gamma_basis), rho_gamma(gamma_basis);
511 xi1_gamma.term(0, 0) = xi1.term(0,0); // this assumes linear expansion
512 xi1_gamma.term(0, 1) = xi1.term(0,1);
513 rho_gamma.term(1, 0) = rho.mean();
514 //rho_gamma.term(1, 1) = rho.standard_deviation();
515 rho_gamma.term(1, 1) = 1.0;
516
517 // Compute gamma pce in transformed basis
518 gammaModel.computeGammaPCE(*gamma_basis, *gamma_quad, xi1_gamma,
519 rho_gamma, gamma_gamma);
520
521 // Project gamma back to original basis
522 Teuchos::Array<double> gamma_point(2);
523 Teuchos::Array<double> gamma_basis_vals(gamma_basis->size());
524 gamma.init(0.0);
525 for (int qp=0; qp<nqp; qp++) {
526 gamma_point[0] = points[qp][0];
527 gamma_point[1] = rho.evaluate(points[qp], basis_vals[qp]);
528 gamma_basis->evaluateBases(gamma_point, gamma_basis_vals);
529 double g = gamma_gamma.evaluate(gamma_point, gamma_basis_vals);
530 for (int i=0; i<sz; i++)
531 gamma[i] += g*weights[qp]*basis_vals[qp][i]/norms[i];
532 }
533
534 // Compute errors
535 rho_error = 0.0;
536 gamma_error = 0.0;
537 for (int i=0; i<sz; i++) {
538 double re = rel_error(rho[i],rho0[i]);
539 double ge = rel_error(gamma[i],gamma0[i]);
540 if (re > rho_error)
541 rho_error = re;
542 if (ge > gamma_error)
543 gamma_error = ge;
544 }
545
546 std::cout << "it = " << it
547 << " rho_error = " << rho_error
548 << " gamma_error = " << gamma_error << std::endl;
549
550 // Update
551 rho0 = rho;
552 gamma0 = gamma;
553 ++it;
554 }
555
556 // Check if we converged
557 if (it >= max_it)
558 throw std::string("model didn't converge!");
559 }
560private:
561 double tol;
565 Teuchos::RCP<const Stokhos::ProductBasis<int,double> > basis;
566 Teuchos::RCP<const Stokhos::Quadrature<int,double> > quad;
567};
568
569int main(int argc, char **argv)
570{
571 try {
572
573 const int n = 102;
574 const double h = 1.0/(n+1);
575 RhoModel rhoModel(n,h);
576 GammaModel gammaModel(n,h);
577 const double tol = 1.0e-7;
578 const int max_it = 20;
579 CoupledSolver coupledSolver(1e-12, max_it, rhoModel, gammaModel);
580
581 // Create product basis
582 const int d = 2;
583 const int p = 10;
584 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
585 for (int i=0; i<d; i++)
586 bases[i] = Teuchos::rcp(new basis_type(p));
587 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
588 Teuchos::rcp(new Stokhos::CompletePolynomialBasis<int,double>(bases));
589
590 // Quadrature
591 Teuchos::RCP<const Stokhos::Quadrature<int,double> > quad =
592 Teuchos::rcp(new Stokhos::TensorProductQuadrature<int,double>(basis));
593
594 // Create expansion for xi1 and xi2
595 Stokhos::OrthogPolyApprox<int,double> x1(basis), x2(basis), rho_pce(basis), gamma_pce(basis);
596
597 // Uniform on [0.2,0.5]
598 x1.term(0,0) = (0.5+0.2)/2.0;
599 x1.term(0,1) = (0.5-0.2)/2.0;
600
601 // Uniform on [0,40]
602 x2.term(1,0) = (40.0+0)/2.0;
603 x2.term(1,1) = (40.0-0)/2.0;
604
605 // Evaluate coupled model at a point
606 Teuchos::Array<double> point(2);
607 point[0] = 0.1234;
608 point[1] = -0.23465;
609 double s1 = x1.evaluate(point);
610 double s2 = x2.evaluate(point);
611 double rho0, gamma0;
612 coupledSolver.solve(s1, s2, rho0, gamma0);
613 std::cout << "s1 = " << s1 << " s2 = " << s2 << std::endl;
614
615 NISPCoupledSolver nispSolver(1e-12, max_it, rhoModel, gammaModel,
616 basis, quad);
617 Stokhos::OrthogPolyApprox<int,double> rho_nisp(basis), gamma_nisp(basis);
618 nispSolver.solve(x1, x2, rho_nisp, gamma_nisp);
619 double rho1 = rho_nisp.evaluate(point);
620 double gamma1 = gamma_nisp.evaluate(point);
621
622 std::cout << "rho_nisp = " << rho1
623 << " rho0 = " << rho0 << " error = "
624 << std::abs(rho1-rho0)/std::abs(rho0) << std::endl;
625 std::cout << "gamma_nisp = " << gamma1
626 << " gamma0 = " << gamma0 << " error = "
627 << std::abs(gamma1-gamma0)/std::abs(gamma0) << std::endl;
628
629 SemiIntrusiveCoupledSolver siSolver(tol, max_it, rhoModel, gammaModel,
630 basis, quad);
631 Stokhos::OrthogPolyApprox<int,double> rho_si(basis), gamma_si(basis);
632 for (int i=0; i<basis->size(); i++) {
633 rho_si[i] = x2[i];
634 gamma_si[i] = x1[i];
635 }
636 siSolver.solve(x1, x2, rho_si, gamma_si);
637 double rho2 = rho_si.evaluate(point);
638 double gamma2 = gamma_si.evaluate(point);
639
640 std::cout << "rho_si = " << rho2
641 << " rho0 = " << rho0 << " error = "
642 << std::abs(rho2-rho0)/std::abs(rho0) << std::endl;
643 std::cout << "gamma_si = " << gamma2
644 << " gamma0 = " << gamma0 << " error = "
645 << std::abs(gamma2-gamma0)/std::abs(gamma0) << std::endl;
646
647 StieltjesCoupledSolver stSolver(tol, max_it, rhoModel, gammaModel,
648 basis, quad);
649 Stokhos::OrthogPolyApprox<int,double> rho_st(basis), gamma_st(basis);
650 for (int i=0; i<basis->size(); i++) {
651 rho_st[i] = x2[i];
652 gamma_st[i] = x1[i];
653 // rho_st[i] = rho_si[i];
654 // gamma_st[i] = gamma_si[i];
655 }
656 stSolver.solve(x1, x2, rho_st, gamma_st);
657 double rho3 = rho_st.evaluate(point);
658 double gamma3 = gamma_st.evaluate(point);
659
660 std::cout << "rho_st = " << rho3
661 << " rho0 = " << rho0 << " error = "
662 << std::abs(rho3-rho0)/std::abs(rho0) << std::endl;
663 std::cout << "gamma_st = " << gamma3
664 << " gamma0 = " << gamma0 << " error = "
665 << std::abs(gamma3-gamma0)/std::abs(gamma0) << std::endl;
666
667 int num_samples = 100;
668 double h_grid = 2.0/(num_samples - 1);
669 std::ofstream coupled("coupled.txt");
670 std::ofstream nisp("nisp.txt");
671 std::ofstream si("si.txt");
672 std::ofstream st("st.txt");
673 for (int i=0; i<num_samples; i++) {
674 point[0] = -1.0 + h_grid*i;
675 for (int j=0; j<num_samples; j++) {
676 point[1] = -1.0 + h_grid*j;
677
678 std::cout << "i = " << i << "j = " << j << std::endl;
679
680 double s1 = x1.evaluate(point);
681 double s2 = x2.evaluate(point);
682 coupledSolver.solve(s1, s2, rho0, gamma0);
683 coupled << s1 << " " << s2 << " " << rho0 << " " << gamma0 << std::endl;
684
685 rho1 = rho_nisp.evaluate(point);
686 gamma1 = gamma_nisp.evaluate(point);
687 nisp << s1 << " " << s2 << " " << rho1 << " " << gamma1 << std::endl;
688
689 rho2 = rho_si.evaluate(point);
690 gamma2 = gamma_si.evaluate(point);
691 si << s1 << " " << s2 << " " << rho2 << " " << gamma2 << std::endl;
692
693 rho3 = rho_st.evaluate(point);
694 gamma3 = gamma_st.evaluate(point);
695 st << s1 << " " << s2 << " " << rho3 << " " << gamma3 << std::endl;
696 }
697 }
698 coupled.close();
699 nisp.close();
700 si.close();
701 st.close();
702
703 }
704 catch (std::string& s) {
705 std::cout << "caught exception: " << s << std::endl;
706 }
707 catch (std::exception& e) {
708 std::cout << e.what() << std::endl;
709 }
710}
expr expr expr dx(i, j)
Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > basis
void solve(double xi1, double xi2, double &rho, double &gamma)
Teuchos::RCP< const Stokhos::Quadrature< int, double > > quad
GammaModel & gammaModel
CoupledSolver(double tol_, int max_it_, RhoModel &rhoModel_, GammaModel &gammaModel_)
Teuchos::Array< int > piv
Teuchos::SerialDenseMatrix< int, double > A
Teuchos::SerialDenseVector< int, double > b
Teuchos::LAPACK< int, double > lapack
void computeGammaPCE(const Stokhos::OrthogPolyBasis< int, double > &basis, const Stokhos::Quadrature< int, double > &quad, const Stokhos::OrthogPolyApprox< int, double > &xi1, const Stokhos::OrthogPolyApprox< int, double > &rho, Stokhos::OrthogPolyApprox< int, double > &gamma)
GammaModel(int n, double dx)
void computeGamma(const Teuchos::Array< double > &s1, const Teuchos::Array< double > &rho, Teuchos::Array< double > &gamma)
Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > basis
void solve(const Stokhos::OrthogPolyApprox< int, double > &xi1, const Stokhos::OrthogPolyApprox< int, double > &xi2, Stokhos::OrthogPolyApprox< int, double > &rho_pce, Stokhos::OrthogPolyApprox< int, double > &gamma_pce)
CoupledSolver coupledSolver
NISPCoupledSolver(double tol, int max_it, RhoModel &rhoModel, GammaModel &gammaModel, const Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > &basis_, const Teuchos::RCP< const Stokhos::Quadrature< int, double > > &quad_)
Teuchos::RCP< const Stokhos::Quadrature< int, double > > quad
Teuchos::Array< int > piv
Teuchos::LAPACK< int, double > lapack
void computeRhoPCE(const Stokhos::OrthogPolyBasis< int, double > &basis, const Stokhos::Quadrature< int, double > &quad, const Stokhos::OrthogPolyApprox< int, double > &xi2, const Stokhos::OrthogPolyApprox< int, double > &gamma, Stokhos::OrthogPolyApprox< int, double > &rho)
void computeRho(const Teuchos::Array< double > &s2, const Teuchos::Array< double > &gamma, Teuchos::Array< double > &rho)
Teuchos::SerialDenseMatrix< int, double > A
Teuchos::SerialDenseVector< int, double > b
RhoModel(int n, double dx)
SemiIntrusiveCoupledSolver(double tol_, int max_it_, RhoModel &rhoModel_, GammaModel &gammaModel_, const Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > &basis_, const Teuchos::RCP< const Stokhos::Quadrature< int, double > > &quad_)
void solve(const Stokhos::OrthogPolyApprox< int, double > &xi1, const Stokhos::OrthogPolyApprox< int, double > &xi2, Stokhos::OrthogPolyApprox< int, double > &rho, Stokhos::OrthogPolyApprox< int, double > &gamma)
Teuchos::RCP< const Stokhos::Quadrature< int, double > > quad
Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > basis
Teuchos::RCP< const Stokhos::Quadrature< int, double > > quad
Teuchos::RCP< const Stokhos::ProductBasis< int, double > > basis
StieltjesCoupledSolver(double tol_, int max_it_, RhoModel &rhoModel_, GammaModel &gammaModel_, const Teuchos::RCP< const Stokhos::ProductBasis< int, double > > &basis_, const Teuchos::RCP< const Stokhos::Quadrature< int, double > > &quad_)
void solve(const Stokhos::OrthogPolyApprox< int, double > &xi1, const Stokhos::OrthogPolyApprox< int, double > &xi2, Stokhos::OrthogPolyApprox< int, double > &rho, Stokhos::OrthogPolyApprox< int, double > &gamma)
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Legendre polynomial basis.
Class to store coefficients of a projection onto an orthogonal polynomial basis.
value_type mean() const
Compute mean of expansion.
void init(const value_type &v)
Initialize coefficients to value.
value_type evaluate(const Teuchos::Array< value_type > &point) const
Evaluate polynomial approximation at a point.
Abstract base class for multivariate orthogonal polynomials.
virtual const Teuchos::Array< value_type > & norm_squared() const =0
Return array storing norm-squared of each basis polynomial.
virtual ordinal_type size() const =0
Return total size of basis.
Abstract base class for multivariate orthogonal polynomials generated from tensor products of univari...
Abstract base class for quadrature methods.
virtual const Teuchos::Array< value_type > & getQuadWeights() const =0
Get quadrature weights.
virtual const Teuchos::Array< Teuchos::Array< value_type > > & getBasisAtQuadPoints() const =0
Get values of basis at quadrature points.
virtual const Teuchos::Array< Teuchos::Array< value_type > > & getQuadPoints() const =0
Get quadrature points.
Generates three-term recurrence using the Discretized Stieltjes procedure applied to a polynomial cha...
Defines quadrature for a tensor product basis by tensor products of 1-D quadrature rules.
Specialization for Sacado::UQ::PCE< Storage<...> >
ScalarType g(const Teuchos::Array< ScalarType > &x, const ScalarType &y)
int main(int argc, char **argv)
double rel_error(double a, double b)
Stokhos::LegendreBasis< int, double > basis_type