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 // $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"
48 #include "Stokhos_StieltjesPCEBasis.hpp"
49 #include "Stokhos_UserDefinedQuadrature.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 
58 typedef Stokhos::LegendreBasis<int,double> basis_type;
59 
60 double rel_error(double a, double b) {
61  return std::abs(a-b)/std::max(std::abs(a), std::abs(b));
62 }
63 
64 class RhoModel {
65 public:
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 
102  void computeRhoPCE(
103  const Stokhos::OrthogPolyBasis<int,double>& basis,
104  const Stokhos::Quadrature<int,double>& quad,
105  const Stokhos::OrthogPolyApprox<int,double>& xi2,
106  const Stokhos::OrthogPolyApprox<int,double>& gamma,
107  Stokhos::OrthogPolyApprox<int,double>& rho) {
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 =
115  quad.getBasisAtQuadPoints();
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 
138 private:
139  Teuchos::SerialDenseMatrix<int,double> A;
140  Teuchos::SerialDenseVector<int,double> b;
141  double h;
142  Teuchos::LAPACK<int,double> lapack;
143  Teuchos::Array<int> piv;
144 };
145 
146 class GammaModel {
147 public:
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 
188  void computeGammaPCE(
189  const Stokhos::OrthogPolyBasis<int,double>& basis,
190  const Stokhos::Quadrature<int,double>& quad,
191  const Stokhos::OrthogPolyApprox<int,double>& xi1,
192  const Stokhos::OrthogPolyApprox<int,double>& rho,
193  Stokhos::OrthogPolyApprox<int,double>& gamma) {
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 =
201  quad.getBasisAtQuadPoints();
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  }
223 private:
224  Teuchos::SerialDenseMatrix<int,double> A;
225  Teuchos::SerialDenseVector<int,double> b;
226  double h;
227  Teuchos::LAPACK<int,double> lapack;
228  Teuchos::Array<int> piv;
229 };
230 
231 class CoupledSolver {
232 public:
233  CoupledSolver(
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  }
274 private:
275  double tol;
276  int max_it;
277  RhoModel& rhoModel;
278  GammaModel& gammaModel;
279  Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> > basis;
280  Teuchos::RCP<const Stokhos::Quadrature<int,double> > quad;
281 };
282 
283 class NISPCoupledSolver {
284 public:
285  NISPCoupledSolver(
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 
293  void solve(const Stokhos::OrthogPolyApprox<int,double>& xi1,
294  const Stokhos::OrthogPolyApprox<int,double>& xi2,
295  Stokhos::OrthogPolyApprox<int,double>& rho_pce,
296  Stokhos::OrthogPolyApprox<int,double>& gamma_pce) {
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  }
327 private:
328  CoupledSolver coupledSolver;
329  Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> > basis;
330  Teuchos::RCP<const Stokhos::Quadrature<int,double> > quad;
331 };
332 
333 class SemiIntrusiveCoupledSolver {
334 public:
335  SemiIntrusiveCoupledSolver(
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 
343  void solve(const Stokhos::OrthogPolyApprox<int,double>& xi1,
344  const Stokhos::OrthogPolyApprox<int,double>& xi2,
345  Stokhos::OrthogPolyApprox<int,double>& rho,
346  Stokhos::OrthogPolyApprox<int,double>& gamma) {
347 
348  Stokhos::OrthogPolyApprox<int,double> rho0(basis);
349  Stokhos::OrthogPolyApprox<int,double> gamma0(basis);
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  }
389 private:
390  double tol;
391  int max_it;
392  RhoModel& rhoModel;
393  GammaModel& gammaModel;
394  Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> > basis;
395  Teuchos::RCP<const Stokhos::Quadrature<int,double> > quad;
396 };
397 
398 class StieltjesCoupledSolver {
399 public:
400  StieltjesCoupledSolver(
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 
408  void solve(const Stokhos::OrthogPolyApprox<int,double>& xi1,
409  const Stokhos::OrthogPolyApprox<int,double>& xi2,
410  Stokhos::OrthogPolyApprox<int,double>& rho,
411  Stokhos::OrthogPolyApprox<int,double>& gamma) {
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 
422  Stokhos::OrthogPolyApprox<int,double> rho0(basis);
423  Stokhos::OrthogPolyApprox<int,double> gamma0(basis);
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] =
451  Teuchos::rcp(new Stokhos::StieltjesPCEBasis<int,double>(
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] =
496  Teuchos::rcp(new Stokhos::StieltjesPCEBasis<int,double>(
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  }
560 private:
561  double tol;
562  int max_it;
563  RhoModel& rhoModel;
564  GammaModel& gammaModel;
565  Teuchos::RCP<const Stokhos::ProductBasis<int,double> > basis;
566  Teuchos::RCP<const Stokhos::Quadrature<int,double> > quad;
567 };
568 
569 int 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 }
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:67
argv
char * argv[]
Definition: Stokhos_HouseTriDiagUnitTest.cpp:286
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:408
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:248
RhoModel::lapack
Teuchos::LAPACK< int, double > lapack
Definition: stieltjes_coupled.cpp:142
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:148
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:400
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:395
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:335
RhoModel::A
Teuchos::SerialDenseMatrix< int, double > A
Definition: stieltjes_coupled.cpp:139
Stokhos::OrthogPolyApprox::init
void init(const value_type &v)
Initialize coefficients to value.
Definition: Stokhos_OrthogPolyApproxImp.hpp:106
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:343
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:60
CoupledSolver::quad
Teuchos::RCP< const Stokhos::Quadrature< int, double > > quad
Definition: stieltjes_coupled.cpp:280
GammaModel::computeGamma
void computeGamma(const Teuchos::Array< double > &s1, const Teuchos::Array< double > &rho, Teuchos::Array< double > &gamma)
Definition: stieltjes_coupled.cpp:149
GammaModel::lapack
Teuchos::LAPACK< int, double > lapack
Definition: stieltjes_coupled.cpp:227
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:188
basis_type
Stokhos::LegendreBasis< int, double > basis_type
Definition: gram_schmidt_example.cpp:52
j
j
Definition: Sacado_Fad_Exp_MP_Vector.hpp:527
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:64
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:143
CoupledSolver::solve
void solve(double xi1, double xi2, double &rho, double &gamma)
Definition: stieltjes_coupled.cpp:238
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:1211
Teuchos::SerialDenseVector::length
OrdinalType length() const
NISPCoupledSolver::coupledSolver
CoupledSolver coupledSolver
Definition: stieltjes_coupled.cpp:328
Sacado::UQ::abs
KOKKOS_INLINE_FUNCTION PCE< Storage > abs(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:1181
CoupledSolver::CoupledSolver
CoupledSolver(double tol_, int max_it_, RhoModel &rhoModel_, GammaModel &gammaModel_)
Definition: stieltjes_coupled.cpp:233
RhoModel::RhoModel
RhoModel(int n, double dx)
Definition: stieltjes_coupled.cpp:66
Stokhos::OrthogPolyApprox::mean
value_type mean() const
Compute mean of expansion.
Definition: Stokhos_OrthogPolyApproxImp.hpp:271
dx
expr expr expr dx(i, j)
Definition: Sacado_Fad_Ops_MP_Vector.hpp:144
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:330
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:102
Sacado::UQ::atan
KOKKOS_INLINE_FUNCTION PCE< Storage > atan(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:1101
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:67
main
int main(int argc, char **argv)
Definition: cijk_ltb_partition.cpp:57
GammaModel::b
Teuchos::SerialDenseVector< int, double > b
Definition: stieltjes_coupled.cpp:225
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:228
Sacado::UQ::sin
KOKKOS_INLINE_FUNCTION PCE< Storage > sin(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:965
StieltjesCoupledSolver::quad
Teuchos::RCP< const Stokhos::Quadrature< int, double > > quad
Definition: stieltjes_coupled.cpp:566
GammaModel::A
Teuchos::SerialDenseMatrix< int, double > A
Definition: stieltjes_coupled.cpp:224
StieltjesCoupledSolver::basis
Teuchos::RCP< const Stokhos::ProductBasis< int, double > > basis
Definition: stieltjes_coupled.cpp:565
RhoModel::b
Teuchos::SerialDenseVector< int, double > b
Definition: stieltjes_coupled.cpp:140
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:293
CoupledSolver::basis
Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > basis
Definition: stieltjes_coupled.cpp:279
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:285
SemiIntrusiveCoupledSolver::basis
Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > basis
Definition: stieltjes_coupled.cpp:394
CoupledSolver::gammaModel
GammaModel & gammaModel
Definition: stieltjes_coupled.cpp:278
n
int n
Stokhos::TensorProductQuadrature
Defines quadrature for a tensor product basis by tensor products of 1-D quadrature rules...
Definition: Stokhos_TensorProductQuadrature.hpp:58
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:109
NISPCoupledSolver::basis
Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > basis
Definition: stieltjes_coupled.cpp:329
Stokhos::OrthogPolyApprox::term
reference term(ordinal_type dimension, ordinal_type order)
Get coefficient term for given dimension and order.
Definition: Stokhos_OrthogPolyApproxImp.hpp:210
Stokhos_StieltjesPCEBasis.hpp
Generated by   doxygen 1.8.5