45 #include "EpetraExt_MatrixMatrix.h" 46 #include "Teuchos_Assert.hpp" 52 struct KL_Diffusion_Func {
54 mutable Teuchos::Array<double> point;
55 Teuchos::RCP< Stokhos::KL::ExponentialRandomField<double> > rf;
57 KL_Diffusion_Func(
double xyLeft,
double xyRight,
58 double mean_,
double std_dev,
59 double L,
int num_KL) : mean(mean_), point(2)
61 Teuchos::ParameterList rfParams;
62 rfParams.set(
"Number of KL Terms", num_KL);
63 rfParams.set(
"Mean", mean);
64 rfParams.set(
"Standard Deviation", std_dev);
66 Teuchos::Array<double> domain_upper(ndim), domain_lower(ndim),
67 correlation_length(ndim);
68 for (
int i=0; i<ndim; i++) {
69 domain_upper[i] = xyRight;
70 domain_lower[i] = xyLeft;
71 correlation_length[i] = L;
73 rfParams.set(
"Domain Upper Bounds", domain_upper);
74 rfParams.set(
"Domain Lower Bounds", domain_lower);
75 rfParams.set(
"Correlation Lengths", correlation_length);
81 ~KL_Diffusion_Func() {
84 double operator() (
double x,
double y,
int k)
const {
89 return rf->evaluate_eigenfunction(point, k-1);
95 template <
typename DiffusionFunc>
96 struct Normalized_KL_Diffusion_Func {
97 const DiffusionFunc& func;
99 Teuchos::Array<double> psi_0, psi_1;
101 Normalized_KL_Diffusion_Func(
102 const DiffusionFunc& func_,
105 d(basis.dimension()),
109 Teuchos::Array<double> zero(d), one(d);
110 for(
int k=0; k<d; k++) {
118 double operator() (
double x,
double y,
int k)
const {
120 double val = func(
x,
y, 0);
121 for (
int i=1; i<=d; i++)
122 val -= psi_0[i]/(psi_1[i]-psi_0[i])*func(
x,
y, i);
127 return 1.0/(psi_1[k]-psi_0[k])*func(
x,
y, k);
132 template <
typename DiffusionFunc>
133 struct LogNormal_Diffusion_Func {
135 const DiffusionFunc& func;
136 Teuchos::RCP<const Stokhos::ProductBasis<int, double> > prodbasis;
138 LogNormal_Diffusion_Func(
140 const DiffusionFunc& func_,
142 : mean(mean_), func(func_), prodbasis(prodbasis_) {}
144 double operator() (
double x,
double y,
int k)
const {
145 int d = prodbasis->dimension();
146 const Teuchos::Array<double>& norms = prodbasis->norm_squared();
148 double sum_g = 0.0, efval;
149 for (
int l=0; l<d; l++) {
152 efval =
std::exp(mean + 0.5*sum_g)/norms[k];
153 for (
int l=0; l<d; l++) {
154 efval *=
std::pow(func(
x,
y,l+1),multiIndex[l]);
164 const Teuchos::RCP<const Epetra_Comm>& comm,
int n,
int d,
169 const Teuchos::RCP<Teuchos::ParameterList>& precParams_) :
172 log_normal(log_normal_),
173 eliminate_bcs(eliminate_bcs_),
174 precParams(precParams_)
176 Kokkos::initialize();
189 h = (xyRight - xyLeft)/((
double)(n-1));
190 Teuchos::Array<int> global_dof_indices;
191 for (
int j=0;
j<n;
j++) {
192 double y = xyLeft +
j*
h;
193 for (
int i=0; i<n; i++) {
194 double x = xyLeft + i*
h;
198 if (i == 0 || i == n-1 ||
j == 0 ||
j == n-1)
199 mesh[idx].boundary =
true;
201 mesh[idx].left = idx-1;
203 mesh[idx].right = idx+1;
205 mesh[idx].down = idx-n;
207 mesh[idx].up = idx+n;
209 global_dof_indices.push_back(idx);
214 int n_global_dof = global_dof_indices.size();
215 int n_proc = comm->NumProc();
216 int proc_id = comm->MyPID();
217 int n_my_dof = n_global_dof / n_proc;
218 if (proc_id == n_proc-1)
219 n_my_dof += n_global_dof % n_proc;
220 int *my_dof = global_dof_indices.getRawPtr() + proc_id*(n_global_dof / n_proc);
222 Teuchos::rcp(
new Epetra_Map(n_global_dof, n_my_dof, my_dof, 0, *comm));
229 p_map = Teuchos::rcp(
new Epetra_LocalMap(d, 0, *comm));
232 g_map = Teuchos::rcp(
new Epetra_LocalMap(1, 0, *comm));
239 p_names = Teuchos::rcp(
new Teuchos::Array<std::string>(d));
240 for (
int i=0;i<d;i++) {
241 std::stringstream ss;
242 ss <<
"KL Random Variable " << i+1;
243 (*p_names)[i] = ss.str();
247 int NumMyElements =
x_map->NumMyElements();
248 int *MyGlobalElements =
x_map->MyGlobalElements();
249 graph = Teuchos::rcp(
new Epetra_CrsGraph(Copy, *
x_map, 5));
250 for (
int i=0; i<NumMyElements; ++i ) {
253 int global_idx = MyGlobalElements[i];
254 graph->InsertGlobalIndices(global_idx, 1, &global_idx);
256 if (!
mesh[global_idx].boundary) {
259 graph->InsertGlobalIndices(global_idx, 1, &
mesh[global_idx].down);
263 graph->InsertGlobalIndices(global_idx, 1, &
mesh[global_idx].left);
267 graph->InsertGlobalIndices(global_idx, 1, &
mesh[global_idx].right);
271 graph->InsertGlobalIndices(global_idx, 1, &
mesh[global_idx].up);
274 graph->FillComplete();
275 graph->OptimizeStorage();
277 KL_Diffusion_Func klFunc(xyLeft, xyRight, mu, s, 1.0, d);
280 if (
basis == Teuchos::null) {
284 Normalized_KL_Diffusion_Func<KL_Diffusion_Func> nklFunc(klFunc, *
basis);
290 int sz =
basis->size();
291 Teuchos::RCP<const Stokhos::ProductBasis<int, double> > prodbasis =
294 LogNormal_Diffusion_Func<KL_Diffusion_Func> lnFunc(mu, klFunc, prodbasis);
299 A = Teuchos::rcp(
new Epetra_CrsMatrix(Copy, *
graph));
302 b = Teuchos::rcp(
new Epetra_Vector(*
x_map));
303 for(
int i=0 ; i<NumMyElements; ++i ) {
304 int global_idx = MyGlobalElements[i];
305 if (
mesh[global_idx].boundary)
311 if (
basis != Teuchos::null) {
317 std::string name =
precParams->get(
"Preconditioner Type",
"Ifpack");
318 Teuchos::RCP<Teuchos::ParameterList> p =
319 Teuchos::rcp(&(
precParams->sublist(
"Preconditioner Parameters")),
false);
333 Teuchos::RCP<const Epetra_Map>
340 Teuchos::RCP<const Epetra_Map>
347 Teuchos::RCP<const Epetra_Map>
351 TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
354 "Error! twoD_diffusion_ME::get_p_map(): " <<
355 "Invalid parameter index l = " << l << std::endl);
360 Teuchos::RCP<const Epetra_Map>
364 TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
367 "Error! twoD_diffusion_ME::get_g_map(): " <<
368 "Invalid parameter index l = " << l << std::endl);
373 Teuchos::RCP<const Teuchos::Array<std::string> >
377 TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
380 "Error! twoD_diffusion_ME::get_p_names(): " <<
381 "Invalid parameter index l = " << l << std::endl);
386 Teuchos::RCP<const Epetra_Vector>
393 Teuchos::RCP<const Epetra_Vector>
397 TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
400 "Error! twoD_diffusion_ME::get_p_init(): " <<
401 "Invalid parameter index l = " << l << std::endl);
406 Teuchos::RCP<Epetra_Operator>
410 Teuchos::RCP<Epetra_CrsMatrix> AA =
411 Teuchos::rcp(
new Epetra_CrsMatrix(Copy, *
graph));
413 AA->OptimizeStorage();
417 Teuchos::RCP<EpetraExt::ModelEvaluator::Preconditioner>
422 Teuchos::RCP<Epetra_Operator> precOp =
424 return Teuchos::rcp(
new EpetraExt::ModelEvaluator::Preconditioner(precOp,
427 return Teuchos::null;
430 EpetraExt::ModelEvaluator::InArgs
435 inArgs.setModelEvalDescription(
"TwoD Diffusion Model Evaluator");
438 inArgs.setSupports(IN_ARG_x,
true);
442 inArgs.setSupports(IN_ARG_x_sg,
true);
443 inArgs.setSupports(IN_ARG_p_sg, 0,
true);
444 inArgs.setSupports(IN_ARG_sg_basis,
true);
445 inArgs.setSupports(IN_ARG_sg_quadrature,
true);
446 inArgs.setSupports(IN_ARG_sg_expansion,
true);
449 inArgs.setSupports(IN_ARG_x_mp,
true);
450 inArgs.setSupports(IN_ARG_p_mp, 0,
true);
455 EpetraExt::ModelEvaluator::OutArgs
459 OutArgsSetup outArgs;
460 outArgs.setModelEvalDescription(
"TwoD Diffusion Model Evaluator");
463 outArgs.set_Np_Ng(1, 1);
464 outArgs.setSupports(OUT_ARG_f,
true);
465 outArgs.setSupports(OUT_ARG_W,
true);
467 outArgs.setSupports(OUT_ARG_WPrec,
true);
470 outArgs.setSupports(OUT_ARG_f_sg,
true);
471 outArgs.setSupports(OUT_ARG_W_sg,
true);
472 outArgs.setSupports(OUT_ARG_g_sg, 0,
true);
475 outArgs.setSupports(OUT_ARG_f_mp,
true);
476 outArgs.setSupports(OUT_ARG_W_mp,
true);
477 outArgs.setSupports(OUT_ARG_g_mp, 0,
true);
484 evalModel(
const InArgs& inArgs,
const OutArgs& outArgs)
const 492 Teuchos::RCP<const Epetra_Vector> det_x = inArgs.get_x();
495 Teuchos::RCP<const Epetra_Vector> p = inArgs.get_p(0);
496 if (p == Teuchos::null)
499 Teuchos::RCP<Epetra_Vector>
f = outArgs.get_f();
500 Teuchos::RCP<Epetra_Operator> W = outArgs.get_W();
501 Teuchos::RCP<Epetra_Operator> WPrec = outArgs.get_WPrec();
502 if (
f != Teuchos::null || W != Teuchos::null || WPrec != Teuchos::null) {
503 if (
basis != Teuchos::null) {
504 for (
int i=0; i<
point.size(); i++)
508 for (
int k=0;k<
A_k.size();k++)
509 EpetraExt::MatrixMatrix::Add((*
A_k[k]),
false,
basis_vals[k], *
A, 1.0);
513 for (
int k=1;k<
A_k.size();k++)
514 EpetraExt::MatrixMatrix::Add((*
A_k[k]),
false, (*p)[k-1], *
A, 1.0);
517 A->OptimizeStorage();
521 if (
f != Teuchos::null) {
522 Teuchos::RCP<Epetra_Vector> kx = Teuchos::rcp(
new Epetra_Vector(*
x_map));
523 A->Apply(*det_x,*kx);
524 f->Update(1.0,*kx,-1.0, *
b, 0.0);
528 if (W != Teuchos::null) {
529 Teuchos::RCP<Epetra_CrsMatrix> jac =
530 Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(W,
true);
533 jac->OptimizeStorage();
537 if (WPrec != Teuchos::null)
541 Teuchos::RCP<Epetra_Vector>
g = outArgs.get_g(0);
542 if (
g != Teuchos::null) {
543 (det_x->MeanValue(&(*
g)[0]));
544 (*g)[0] *=
double(det_x->GlobalLength()) /
double(
mesh.size());
552 InArgs::sg_const_vector_t x_sg = inArgs.get_x_sg();
555 InArgs::sg_const_vector_t p_sg = inArgs.get_p_sg(0);
558 OutArgs::sg_vector_t f_sg = outArgs.get_f_sg();
559 if (f_sg != Teuchos::null) {
562 Teuchos::RCP<Stokhos::OrthogPolyExpansion<int,double> > expn =
563 inArgs.get_sg_expansion();
565 Teuchos::RCP<const Cijk_type> Cijk = expn->getTripleProduct();
566 const Teuchos::Array<double>& norms =
basis->norm_squared();
570 for (
int i=0;i<
basis->size();i++) {
576 Cijk_type::k_iterator k_begin = Cijk->k_begin();
577 Cijk_type::k_iterator k_end = Cijk->k_end();
578 for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
579 int k = Stokhos::index(k_it);
580 for (Cijk_type::kj_iterator j_it = Cijk->j_begin(k_it);
581 j_it != Cijk->j_end(k_it); ++j_it) {
582 int j = Stokhos::index(j_it);
585 for (Cijk_type::kj_iterator j_it = Cijk->j_begin(k_it);
586 j_it != Cijk->j_end(k_it); ++j_it) {
587 int j = Stokhos::index(j_it);
588 for (Cijk_type::kji_iterator i_it = Cijk->i_begin(j_it);
589 i_it != Cijk->i_end(j_it); ++i_it) {
590 int i = Stokhos::index(i_it);
591 double c = Stokhos::value(i_it);
596 (*f_sg)[0].Update(-1.0,*
b,1.0);
600 OutArgs::sg_operator_t W_sg = outArgs.get_W_sg();
601 if (W_sg != Teuchos::null) {
602 for (
int i=0; i<W_sg->size(); i++) {
603 Teuchos::RCP<Epetra_CrsMatrix> jac =
604 Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(W_sg->getCoeffPtr(i),
true);
607 jac->OptimizeStorage();
612 Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > g_sg =
614 if (g_sg != Teuchos::null) {
615 int sz = x_sg->size();
616 for (
int i=0; i<sz; i++) {
617 (*x_sg)[i].MeanValue(&(*g_sg)[i][0]);
627 mp_const_vector_t x_mp = inArgs.get_x_mp();
630 mp_const_vector_t p_mp = inArgs.get_p_mp(0);
633 mp_vector_t f_mp = outArgs.get_f_mp();
634 mp_operator_t W_mp = outArgs.get_W_mp();
635 if (f_mp != Teuchos::null || W_mp != Teuchos::null) {
636 int num_mp = x_mp->size();
637 for (
int i=0; i<num_mp; i++) {
639 if (
basis != Teuchos::null) {
640 for (
int k=0; k<
point.size(); k++)
641 point[k] = (*p_mp)[i][k];
644 for (
int k=0;k<
A_k.size();k++)
650 for (
int k=1;k<
A_k.size();k++)
651 EpetraExt::MatrixMatrix::Add((*
A_k[k]),
false, (*p_mp)[i][k-1], *
A,
656 A->OptimizeStorage();
659 if (f_mp != Teuchos::null) {
660 A->Apply((*x_mp)[i], (*f_mp)[i]);
661 (*f_mp)[i].Update(-1.0, *
b, 1.0);
665 if (W_mp != Teuchos::null) {
666 Teuchos::RCP<Epetra_CrsMatrix> jac =
667 Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(W_mp->getCoeffPtr(i),
671 jac->OptimizeStorage();
677 mp_vector_t g_mp = outArgs.get_g_mp(0);
678 if (g_mp != Teuchos::null) {
679 int sz = x_mp->size();
680 for (
int i=0; i<sz; i++) {
681 (*x_mp)[i].MeanValue(&(*g_mp)[i][0]);
KOKKOS_INLINE_FUNCTION PCE< Storage > pow(const PCE< Storage > &a, const PCE< Storage > &b)
Teuchos::Array< Teuchos::RCP< Epetra_CrsMatrix > > A_k
KL coefficients of operator.
Teuchos::RCP< Epetra_Vector > x_init
Initial guess.
Teuchos::RCP< Epetra_Map > p_map
Parameter vector map.
Teuchos::RCP< const Epetra_Map > get_f_map() const
Return residual vector map.
virtual void evaluateBases(const Teuchos::ArrayView< const value_type > &point, Teuchos::Array< value_type > &basis_vals) const =0
Evaluate basis polynomials at given point point.
Teuchos::RCP< const Epetra_Vector > get_p_init(int l) const
Return initial parameters.
A multidimensional index.
Teuchos::RCP< Stokhos::AbstractPreconditionerFactory > precFactory
Teuchos::RCP< Epetra_CrsGraph > graph
Jacobian graph.
Teuchos::RCP< const Epetra_Map > get_g_map(int j) const
Return response function map.
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
Teuchos::Array< MeshPoint > mesh
Teuchos::RCP< Teuchos::ParameterList > precParams
Teuchos::Array< double > basis_vals
Array to store values of basis at a point.
Teuchos::RCP< const Epetra_Map > get_x_map() const
Return solution vector map.
Abstract base class for multivariate orthogonal polynomials generated from tensor products of univari...
~twoD_diffusion_ME()
Destructor.
KOKKOS_INLINE_FUNCTION PCE< Storage > exp(const PCE< Storage > &a)
Teuchos::RCP< Epetra_Map > x_map
Solution vector map.
Stokhos::Sparse3Tensor< int, double > Cijk_type
Teuchos::RCP< Epetra_CrsMatrix > A
Matrix to store deterministic operator.
twoD_diffusion_ME(const Teuchos::RCP< const Epetra_Comm > &comm, int n, int d, double s=0.1, double mu=0.2, const Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > &basis=Teuchos::null, bool log_normal=false, bool eliminate_bcs=false, const Teuchos::RCP< Teuchos::ParameterList > &precParams=Teuchos::null)
Constructor.
An class for building preconditioners.
void evalModel(const InArgs &inArgs, const OutArgs &outArgs) const
Evaluate model on InArgs.
Teuchos::RCP< Epetra_Vector > p_init
Initial parameters.
ScalarType f(const Teuchos::Array< ScalarType > &x, double a, double b)
Teuchos::RCP< Epetra_Operator > create_W() const
Create W = alpha*M + beta*J matrix.
InArgs createInArgs() const
Create InArgs.
Teuchos::RCP< const Epetra_Vector > get_x_init() const
Return initial solution.
Class representing a KL expansion of an exponential random field.
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int l) const
Return array of parameter names.
Teuchos::RCP< const Epetra_Map > get_p_map(int l) const
Return parameter vector map.
Teuchos::RCP< Epetra_Vector > b
Deterministic RHS.
Teuchos::RCP< Teuchos::Array< std::string > > p_names
Parameter names.
OutArgs createOutArgs() const
Create OutArgs.
Teuchos::Array< double > point
Array to store a point for basis evaluation.
void fillMatrices(const FuncT &func, int sz)
Fill coefficient matrix given function to evaluate diffusion coefficient.
Teuchos::RCP< Epetra_Map > g_map
Response vector map.
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType ValueType * y
Teuchos::RCP< EpetraExt::ModelEvaluator::Preconditioner > create_WPrec() const
Create preconditioner for W.
Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > basis
Teuchos::Array< Teuchos::RCP< Epetra_Vector > > sg_kx_vec_all
Vectors to store matrix-vector products in SG residual calculation.
ScalarType g(const Teuchos::Array< ScalarType > &x, const ScalarType &y)