51#include "Teuchos_SerialDenseMatrix.hpp"
52#include "Teuchos_SerialDenseVector.hpp"
53#include "Teuchos_LAPACK.hpp"
61 return std::abs(a-b)/std::max(std::abs(a), std::abs(b));
68 const Teuchos::Array<double>& gamma,
69 Teuchos::Array<double>& rho) {
73 for (
int qp=0; qp<m; qp++) {
78 A(0,1) = s2[qp]/(2.*
h) + 1.0/(
h*
h);
79 for (
int i=1; i<n-1; i++) {
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);
84 A(n-1,n-2) = -s2[qp]/(2.*
h) + 1.0/(
h*
h);
85 A(n-1,n-1) = -2.0/(
h*
h);
88 b.putScalar(gamma[qp]);
92 lapack.GESV(n, 1,
A.values(), n, &
piv[0],
b.values(), n, &info);
96 for (
int i=0; i<n; i++)
110 const Teuchos::Array<double>& weights =
112 const Teuchos::Array<Teuchos::Array<double> >& points =
114 const Teuchos::Array<Teuchos::Array<double> >& basis_vals =
116 const Teuchos::Array<double>& norms = basis.
norm_squared();
117 int nqp = weights.size();
118 int sz = basis.
size();
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]);
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];
139 Teuchos::SerialDenseMatrix<int,double>
A;
140 Teuchos::SerialDenseVector<int,double>
b;
150 const Teuchos::Array<double>& rho,
151 Teuchos::Array<double>& gamma) {
155 for (
int qp=0; qp<m; qp++) {
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);
166 A(n-1,n-2) = s1[qp]/(
h*
h);
167 A(n-1,n-1) = -s1[qp]*2.0/(
h*
h) + rho[qp];
170 double pi = 4.0*std::atan(1.0);
171 for (
int i=0; i<n; i++) {
173 b(i) = std::sin(2.0*pi*x);
178 lapack.GESV(n, 1,
A.values(), n, &
piv[0],
b.values(), n, &info);
182 for (
int i=0; i<n; i++)
183 gamma[qp] +=
b(i)*
b(i);
184 gamma[qp] *= 100.0*
h;
196 const Teuchos::Array<double>& weights =
198 const Teuchos::Array<Teuchos::Array<double> >& points =
200 const Teuchos::Array<Teuchos::Array<double> >& basis_vals =
202 const Teuchos::Array<double>& norms = basis.
norm_squared();
203 int nqp = weights.size();
204 int sz = basis.
size();
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]);
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];
224 Teuchos::SerialDenseMatrix<int,double>
A;
225 Teuchos::SerialDenseVector<int,double>
b;
234 double tol_,
int max_it_,
238 void solve(
double xi1,
double xi2,
double& rho,
double& gamma) {
241 double rho_error = 1.0;
242 double gamma_error = 1.0;
244 Teuchos::Array<double> s1(1), s2(1), r(1),
g(1);
249 while( (rho_error >
tol || gamma_error >
tol) && it <
max_it) {
269 throw std::string(
"model didn't converge!");
279 Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> >
basis;
280 Teuchos::RCP<const Stokhos::Quadrature<int,double> >
quad;
286 double tol,
int max_it,
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();
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]);
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];
329 Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> >
basis;
330 Teuchos::RCP<const Stokhos::Quadrature<int,double> >
quad;
336 double tol_,
int max_it_,
350 double rho_error = 1.0;
351 double gamma_error = 1.0;
353 int sz =
basis->size();
355 while( (rho_error >
tol || gamma_error >
tol) && it <
max_it) {
366 for (
int i=0; i<sz; i++) {
368 double ge =
rel_error(gamma[i],gamma0[i]);
371 if (ge > gamma_error)
375 std::cout <<
"it = " << it
376 <<
" rho_error = " << rho_error
377 <<
" gamma_error = " << gamma_error << std::endl;
387 throw std::string(
"model didn't converge!");
394 Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> >
basis;
395 Teuchos::RCP<const Stokhos::Quadrature<int,double> >
quad;
401 double tol_,
int max_it_,
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();
424 double rho_error = 1.0;
425 double gamma_error = 1.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;
435 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk;
436 if (project_integrals)
437 Cijk =
basis->computeTripleProductTensor();
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;
445 while( (rho_error >
tol || gamma_error >
tol) && it <
max_it) {
448 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > rho_bases(2);
449 rho_bases[0] = coordinate_bases[1];
452 p, Teuchos::rcp(&gamma,
false), st_quad,
454 normalize, project_integrals, Cijk));
455 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> >
460 Teuchos::RCP<const Stokhos::Quadrature<int,double> > rho_quad =
465 gamma_rho(rho_basis), rho_rho(rho_basis);
466 xi2_rho.term(0, 0) = xi2.term(1,0);
467 xi2_rho.term(0, 1) = xi2.term(1,1);
468 gamma_rho.term(1, 0) = gamma.
mean();
470 gamma_rho.term(1, 1) = 1.0;
480 Teuchos::Array<double> rho_point(2);
481 Teuchos::Array<double> rho_basis_vals(rho_basis->size());
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];
493 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > gamma_bases(2);
494 gamma_bases[0] = coordinate_bases[0];
497 p, Teuchos::rcp(&rho,
false), st_quad,
499 normalize, project_integrals, Cijk));
500 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> >
505 Teuchos::RCP<const Stokhos::Quadrature<int,double> > gamma_quad =
510 gamma_gamma(gamma_basis), rho_gamma(gamma_basis);
511 xi1_gamma.term(0, 0) = xi1.term(0,0);
512 xi1_gamma.term(0, 1) = xi1.term(0,1);
513 rho_gamma.term(1, 0) = rho.mean();
515 rho_gamma.term(1, 1) = 1.0;
519 rho_gamma, gamma_gamma);
522 Teuchos::Array<double> gamma_point(2);
523 Teuchos::Array<double> gamma_basis_vals(gamma_basis->size());
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];
537 for (
int i=0; i<sz; i++) {
539 double ge =
rel_error(gamma[i],gamma0[i]);
542 if (ge > gamma_error)
546 std::cout <<
"it = " << it
547 <<
" rho_error = " << rho_error
548 <<
" gamma_error = " << gamma_error << std::endl;
558 throw std::string(
"model didn't converge!");
565 Teuchos::RCP<const Stokhos::ProductBasis<int,double> >
basis;
566 Teuchos::RCP<const Stokhos::Quadrature<int,double> >
quad;
574 const double h = 1.0/(n+1);
577 const double tol = 1.0e-7;
578 const int max_it = 20;
579 CoupledSolver coupledSolver(1e-12, max_it, rhoModel, gammaModel);
584 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
585 for (
int i=0; i<d; i++)
587 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
591 Teuchos::RCP<const Stokhos::Quadrature<int,double> > quad =
598 x1.term(0,0) = (0.5+0.2)/2.0;
599 x1.term(0,1) = (0.5-0.2)/2.0;
602 x2.term(1,0) = (40.0+0)/2.0;
603 x2.term(1,1) = (40.0-0)/2.0;
606 Teuchos::Array<double> point(2);
609 double s1 = x1.evaluate(point);
610 double s2 = x2.evaluate(point);
612 coupledSolver.
solve(s1, s2, rho0, gamma0);
613 std::cout <<
"s1 = " << s1 <<
" s2 = " << s2 << std::endl;
618 nispSolver.solve(x1, x2, rho_nisp, gamma_nisp);
619 double rho1 = rho_nisp.evaluate(point);
620 double gamma1 = gamma_nisp.evaluate(point);
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;
632 for (
int i=0; i<basis->size(); i++) {
636 siSolver.solve(x1, x2, rho_si, gamma_si);
637 double rho2 = rho_si.evaluate(point);
638 double gamma2 = gamma_si.evaluate(point);
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;
650 for (
int i=0; i<basis->size(); i++) {
656 stSolver.solve(x1, x2, rho_st, gamma_st);
657 double rho3 = rho_st.evaluate(point);
658 double gamma3 = gamma_st.evaluate(point);
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;
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;
678 std::cout <<
"i = " << i <<
"j = " <<
j << std::endl;
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;
685 rho1 = rho_nisp.evaluate(point);
686 gamma1 = gamma_nisp.evaluate(point);
687 nisp << s1 <<
" " << s2 <<
" " << rho1 <<
" " << gamma1 << std::endl;
689 rho2 = rho_si.evaluate(point);
690 gamma2 = gamma_si.evaluate(point);
691 si << s1 <<
" " << s2 <<
" " << rho2 <<
" " << gamma2 << std::endl;
693 rho3 = rho_st.evaluate(point);
694 gamma3 = gamma_st.evaluate(point);
695 st << s1 <<
" " << s2 <<
" " << rho3 <<
" " << gamma3 << std::endl;
704 catch (std::string& s) {
705 std::cout <<
"caught exception: " << s << std::endl;
707 catch (std::exception& e) {
708 std::cout << e.what() << std::endl;
Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > basis
void solve(double xi1, double xi2, double &rho, double &gamma)
Teuchos::RCP< const Stokhos::Quadrature< int, double > > quad
CoupledSolver(double tol_, int max_it_, RhoModel &rhoModel_, GammaModel &gammaModel_)
Teuchos::Array< int > piv
Teuchos::SerialDenseMatrix< int, double > A
Teuchos::SerialDenseVector< int, double > b
Teuchos::LAPACK< int, double > lapack
void computeGammaPCE(const Stokhos::OrthogPolyBasis< int, double > &basis, const Stokhos::Quadrature< int, double > &quad, const Stokhos::OrthogPolyApprox< int, double > &xi1, const Stokhos::OrthogPolyApprox< int, double > &rho, Stokhos::OrthogPolyApprox< int, double > &gamma)
GammaModel(int n, double dx)
void computeGamma(const Teuchos::Array< double > &s1, const Teuchos::Array< double > &rho, Teuchos::Array< double > &gamma)
Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > basis
void solve(const Stokhos::OrthogPolyApprox< int, double > &xi1, const Stokhos::OrthogPolyApprox< int, double > &xi2, Stokhos::OrthogPolyApprox< int, double > &rho_pce, Stokhos::OrthogPolyApprox< int, double > &gamma_pce)
CoupledSolver coupledSolver
NISPCoupledSolver(double tol, int max_it, RhoModel &rhoModel, GammaModel &gammaModel, const Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > &basis_, const Teuchos::RCP< const Stokhos::Quadrature< int, double > > &quad_)
Teuchos::RCP< const Stokhos::Quadrature< int, double > > quad
Teuchos::Array< int > piv
Teuchos::LAPACK< int, double > lapack
void computeRhoPCE(const Stokhos::OrthogPolyBasis< int, double > &basis, const Stokhos::Quadrature< int, double > &quad, const Stokhos::OrthogPolyApprox< int, double > &xi2, const Stokhos::OrthogPolyApprox< int, double > &gamma, Stokhos::OrthogPolyApprox< int, double > &rho)
void computeRho(const Teuchos::Array< double > &s2, const Teuchos::Array< double > &gamma, Teuchos::Array< double > &rho)
Teuchos::SerialDenseMatrix< int, double > A
Teuchos::SerialDenseVector< int, double > b
RhoModel(int n, double dx)
SemiIntrusiveCoupledSolver(double tol_, int max_it_, RhoModel &rhoModel_, GammaModel &gammaModel_, const Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > &basis_, const Teuchos::RCP< const Stokhos::Quadrature< int, double > > &quad_)
void solve(const Stokhos::OrthogPolyApprox< int, double > &xi1, const Stokhos::OrthogPolyApprox< int, double > &xi2, Stokhos::OrthogPolyApprox< int, double > &rho, Stokhos::OrthogPolyApprox< int, double > &gamma)
Teuchos::RCP< const Stokhos::Quadrature< int, double > > quad
Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > basis
Teuchos::RCP< const Stokhos::Quadrature< int, double > > quad
Teuchos::RCP< const Stokhos::ProductBasis< int, double > > basis
StieltjesCoupledSolver(double tol_, int max_it_, RhoModel &rhoModel_, GammaModel &gammaModel_, const Teuchos::RCP< const Stokhos::ProductBasis< int, double > > &basis_, const Teuchos::RCP< const Stokhos::Quadrature< int, double > > &quad_)
void solve(const Stokhos::OrthogPolyApprox< int, double > &xi1, const Stokhos::OrthogPolyApprox< int, double > &xi2, Stokhos::OrthogPolyApprox< int, double > &rho, Stokhos::OrthogPolyApprox< int, double > &gamma)
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Legendre polynomial basis.
Class to store coefficients of a projection onto an orthogonal polynomial basis.
value_type mean() const
Compute mean of expansion.
void init(const value_type &v)
Initialize coefficients to value.
value_type evaluate(const Teuchos::Array< value_type > &point) const
Evaluate polynomial approximation at a point.
Abstract base class for multivariate orthogonal polynomials.
virtual const Teuchos::Array< value_type > & norm_squared() const =0
Return array storing norm-squared of each basis polynomial.
virtual ordinal_type size() const =0
Return total size of basis.
Abstract base class for multivariate orthogonal polynomials generated from tensor products of univari...
Abstract base class for quadrature methods.
virtual const Teuchos::Array< value_type > & getQuadWeights() const =0
Get quadrature weights.
virtual const Teuchos::Array< Teuchos::Array< value_type > > & getBasisAtQuadPoints() const =0
Get values of basis at quadrature points.
virtual const Teuchos::Array< Teuchos::Array< value_type > > & getQuadPoints() const =0
Get quadrature points.
Generates three-term recurrence using the Discretized Stieltjes procedure applied to a polynomial cha...
Defines quadrature for a tensor product basis by tensor products of 1-D quadrature rules.
Specialization for Sacado::UQ::PCE< Storage<...> >
ScalarType g(const Teuchos::Array< ScalarType > &x, const ScalarType &y)
int main(int argc, char **argv)
double rel_error(double a, double b)
Stokhos::LegendreBasis< int, double > basis_type