51 #include "Teuchos_SerialDenseMatrix.hpp" 52 #include "Teuchos_SerialDenseVector.hpp" 53 #include "Teuchos_LAPACK.hpp" 65 const Teuchos::Array<double>& gamma,
66 Teuchos::Array<double>& rho) {
70 for (
int qp=0; qp<m; qp++) {
75 A(0,1) = s2[qp]/(2.*
h) + 1.0/(
h*
h);
76 for (
int i=1; i<n-1; i++) {
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);
81 A(n-1,n-2) = -s2[qp]/(2.*
h) + 1.0/(
h*
h);
82 A(n-1,n-1) = -2.0/(
h*
h);
85 b.putScalar(gamma[qp]);
89 lapack.GESV(n, 1,
A.values(), n, &
piv[0],
b.values(), n, &info);
93 for (
int i=0; i<n; i++)
107 const Teuchos::Array<double>& weights =
109 const Teuchos::Array<Teuchos::Array<double> >& points =
111 const Teuchos::Array<Teuchos::Array<double> >& basis_vals =
113 const Teuchos::Array<double>& norms = basis.
norm_squared();
114 int nqp = weights.size();
115 int sz = basis.
size();
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]);
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];
136 Teuchos::SerialDenseMatrix<int,double>
A;
137 Teuchos::SerialDenseVector<int,double>
b;
147 const Teuchos::Array<double>& rho,
148 Teuchos::Array<double>& gamma) {
152 for (
int qp=0; qp<m; qp++) {
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);
163 A(n-1,n-2) = s1[qp]/(
h*
h);
164 A(n-1,n-1) = -s1[qp]*2.0/(
h*
h) + rho[qp];
168 for (
int i=0; i<n; i++) {
175 lapack.GESV(n, 1,
A.values(), n, &
piv[0],
b.values(), n, &info);
179 for (
int i=0; i<n; i++)
180 gamma[qp] +=
b(i)*
b(i);
181 gamma[qp] *= 100.0*
h;
193 const Teuchos::Array<double>& weights =
195 const Teuchos::Array<Teuchos::Array<double> >& points =
197 const Teuchos::Array<Teuchos::Array<double> >& basis_vals =
199 const Teuchos::Array<double>& norms = basis.
norm_squared();
200 int nqp = weights.size();
201 int sz = basis.
size();
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]);
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];
221 Teuchos::SerialDenseMatrix<int,double>
A;
222 Teuchos::SerialDenseVector<int,double>
b;
231 double tol_,
int max_it_,
235 void solve(
double xi1,
double xi2,
double& rho,
double& gamma) {
238 double rho_error = 1.0;
239 double gamma_error = 1.0;
241 Teuchos::Array<double> s1(1), s2(1), r(1),
g(1);
246 while( (rho_error >
tol || gamma_error >
tol) && it <
max_it) {
266 throw std::string(
"model didn't converge!");
276 Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> >
basis;
277 Teuchos::RCP<const Stokhos::Quadrature<int,double> >
quad;
283 double tol,
int max_it,
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();
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]);
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];
326 Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> >
basis;
327 Teuchos::RCP<const Stokhos::Quadrature<int,double> >
quad;
333 double tol_,
int max_it_,
347 double rho_error = 1.0;
348 double gamma_error = 1.0;
350 int sz =
basis->size();
352 while( (rho_error >
tol || gamma_error >
tol) && it <
max_it) {
363 for (
int i=0; i<sz; i++) {
365 double ge =
rel_error(gamma[i],gamma0[i]);
368 if (ge > gamma_error)
372 std::cout <<
"it = " << it
373 <<
" rho_error = " << rho_error
374 <<
" gamma_error = " << gamma_error << std::endl;
384 throw std::string(
"model didn't converge!");
391 Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> >
basis;
392 Teuchos::RCP<const Stokhos::Quadrature<int,double> >
quad;
398 double tol_,
int max_it_,
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();
421 double rho_error = 1.0;
422 double gamma_error = 1.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;
432 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk;
433 if (project_integrals)
434 Cijk =
basis->computeTripleProductTensor();
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;
442 while( (rho_error >
tol || gamma_error >
tol) && it <
max_it) {
445 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > rho_bases(2);
446 rho_bases[0] = coordinate_bases[1];
449 p, Teuchos::rcp(&gamma,
false), st_quad,
451 normalize, project_integrals, Cijk));
452 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> >
457 Teuchos::RCP<const Stokhos::Quadrature<int,double> > rho_quad =
462 gamma_rho(rho_basis), rho_rho(rho_basis);
463 xi2_rho.term(0, 0) = xi2.
term(1,0);
464 xi2_rho.term(0, 1) = xi2.
term(1,1);
465 gamma_rho.term(1, 0) = gamma.
mean();
467 gamma_rho.term(1, 1) = 1.0;
477 Teuchos::Array<double> rho_point(2);
478 Teuchos::Array<double> rho_basis_vals(rho_basis->size());
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];
490 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > gamma_bases(2);
491 gamma_bases[0] = coordinate_bases[0];
494 p, Teuchos::rcp(&rho,
false), st_quad,
496 normalize, project_integrals, Cijk));
497 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> >
502 Teuchos::RCP<const Stokhos::Quadrature<int,double> > gamma_quad =
507 gamma_gamma(gamma_basis), rho_gamma(gamma_basis);
508 xi1_gamma.term(0, 0) = xi1.
term(0,0);
509 xi1_gamma.term(0, 1) = xi1.
term(0,1);
512 rho_gamma.
term(1, 1) = 1.0;
516 rho_gamma, gamma_gamma);
519 Teuchos::Array<double> gamma_point(2);
520 Teuchos::Array<double> gamma_basis_vals(gamma_basis->size());
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];
534 for (
int i=0; i<sz; i++) {
536 double ge =
rel_error(gamma[i],gamma0[i]);
539 if (ge > gamma_error)
543 std::cout <<
"it = " << it
544 <<
" rho_error = " << rho_error
545 <<
" gamma_error = " << gamma_error << std::endl;
555 throw std::string(
"model didn't converge!");
562 Teuchos::RCP<const Stokhos::ProductBasis<int,double> >
basis;
563 Teuchos::RCP<const Stokhos::Quadrature<int,double> >
quad;
571 const double h = 1.0/(n+1);
574 const double tol = 1.0e-7;
575 const int max_it = 20;
576 CoupledSolver coupledSolver(1e-12, max_it, rhoModel, gammaModel);
581 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
582 for (
int i=0; i<d; i++)
584 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
588 Teuchos::RCP<const Stokhos::Quadrature<int,double> > quad =
595 x1.term(0,0) = (0.5+0.2)/2.0;
596 x1.term(0,1) = (0.5-0.2)/2.0;
599 x2.term(1,0) = (40.0+0)/2.0;
600 x2.term(1,1) = (40.0-0)/2.0;
603 Teuchos::Array<double> point(2);
606 double s1 = x1.evaluate(point);
607 double s2 = x2.evaluate(point);
609 coupledSolver.
solve(s1, s2, rho0, gamma0);
610 std::cout <<
"s1 = " << s1 <<
" s2 = " << s2 << std::endl;
615 nispSolver.
solve(x1, x2, rho_nisp, gamma_nisp);
616 double rho1 = rho_nisp.evaluate(point);
617 double gamma1 = gamma_nisp.
evaluate(point);
619 std::cout <<
"rho_nisp = " << rho1
620 <<
" rho0 = " << rho0 <<
" error = " 622 std::cout <<
"gamma_nisp = " << gamma1
623 <<
" gamma0 = " << gamma0 <<
" error = " 629 for (
int i=0; i<basis->size(); i++) {
633 siSolver.
solve(x1, x2, rho_si, gamma_si);
634 double rho2 = rho_si.evaluate(point);
635 double gamma2 = gamma_si.
evaluate(point);
637 std::cout <<
"rho_si = " << rho2
638 <<
" rho0 = " << rho0 <<
" error = " 640 std::cout <<
"gamma_si = " << gamma2
641 <<
" gamma0 = " << gamma0 <<
" error = " 647 for (
int i=0; i<basis->size(); i++) {
653 stSolver.
solve(x1, x2, rho_st, gamma_st);
654 double rho3 = rho_st.evaluate(point);
655 double gamma3 = gamma_st.
evaluate(point);
657 std::cout <<
"rho_st = " << rho3
658 <<
" rho0 = " << rho0 <<
" error = " 660 std::cout <<
"gamma_st = " << gamma3
661 <<
" gamma0 = " << gamma0 <<
" error = " 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;
675 std::cout <<
"i = " << i <<
"j = " <<
j << std::endl;
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;
682 rho1 = rho_nisp.evaluate(point);
683 gamma1 = gamma_nisp.
evaluate(point);
684 nisp << s1 <<
" " << s2 <<
" " << rho1 <<
" " << gamma1 << std::endl;
686 rho2 = rho_si.evaluate(point);
688 si << s1 <<
" " << s2 <<
" " << rho2 <<
" " << gamma2 << std::endl;
690 rho3 = rho_st.evaluate(point);
692 st << s1 <<
" " << s2 <<
" " << rho3 <<
" " << gamma3 << std::endl;
701 catch (std::string& s) {
702 std::cout <<
"caught exception: " << s << std::endl;
704 catch (std::exception& e) {
705 std::cout << e.what() << std::endl;
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
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_)
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_)
Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > basis
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.