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