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