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 double operator() (
double x,
double y,
int k)
const {
86 return rf->evaluate_eigenfunction(point, k-1);
92 template <
typename DiffusionFunc>
93 struct Normalized_KL_Diffusion_Func {
94 const DiffusionFunc& func;
96 Teuchos::Array<double> psi_0, psi_1;
98 Normalized_KL_Diffusion_Func(
99 const DiffusionFunc& func_,
102 d(basis.dimension()),
106 Teuchos::Array<double> zero(d), one(d);
107 for(
int k=0; k<d; k++) {
115 double operator() (
double x,
double y,
int k)
const {
117 double val = func(
x,
y, 0);
118 for (
int i=1; i<=d; i++)
119 val -= psi_0[i]/(psi_1[i]-psi_0[i])*func(
x,
y, i);
124 return 1.0/(psi_1[k]-psi_0[k])*func(
x,
y, k);
129 template <
typename DiffusionFunc>
130 struct LogNormal_Diffusion_Func {
132 const DiffusionFunc& func;
133 Teuchos::RCP<const Stokhos::ProductBasis<int, double> > prodbasis;
135 LogNormal_Diffusion_Func(
137 const DiffusionFunc& func_,
139 : mean(mean_), func(func_), prodbasis(prodbasis_) {}
141 double operator() (
double x,
double y,
int k)
const {
142 int d = prodbasis->dimension();
143 const Teuchos::Array<double>& norms = prodbasis->norm_squared();
144 Teuchos::Array<int> multiIndex = prodbasis->getTerm(k);
145 double sum_g = 0.0, efval;
146 for (
int l=0; l<d; l++) {
149 efval =
std::exp(mean + 0.5*sum_g)/norms[k];
150 for (
int l=0; l<d; l++) {
151 efval *=
std::pow(func(
x,
y,l+1),multiIndex[l]);
161 const Teuchos::RCP<const Epetra_Comm>& comm,
int n,
int d,
165 bool eliminate_bcs_) :
168 log_normal(log_normal_),
169 eliminate_bcs(eliminate_bcs_)
182 h = (xyRight - xyLeft)/((
double)(n-1));
183 Teuchos::Array<int> global_dof_indices;
184 for (
int j=0;
j<n;
j++) {
185 double y = xyLeft +
j*
h;
186 for (
int i=0; i<n; i++) {
187 double x = xyLeft + i*
h;
191 if (i == 0 || i == n-1 ||
j == 0 ||
j == n-1)
192 mesh[idx].boundary =
true;
194 mesh[idx].left = idx-1;
196 mesh[idx].right = idx+1;
198 mesh[idx].down = idx-n;
200 mesh[idx].up = idx+n;
202 global_dof_indices.push_back(idx);
207 int n_global_dof = global_dof_indices.size();
208 int n_proc = comm->NumProc();
209 int proc_id = comm->MyPID();
210 int n_my_dof = n_global_dof / n_proc;
211 if (proc_id == n_proc-1)
212 n_my_dof += n_global_dof % n_proc;
213 int *my_dof = global_dof_indices.getRawPtr() + proc_id*(n_global_dof / n_proc);
215 Teuchos::rcp(
new Epetra_Map(n_global_dof, n_my_dof, my_dof, 0, *comm));
222 p_map = Teuchos::rcp(
new Epetra_LocalMap(d, 0, *comm));
225 g_map = Teuchos::rcp(
new Epetra_LocalMap(1, 0, *comm));
232 p_names = Teuchos::rcp(
new Teuchos::Array<std::string>(d));
233 for (
int i=0;i<d;i++) {
234 std::stringstream ss;
235 ss <<
"KL Random Variable " << i+1;
236 (*p_names)[i] = ss.str();
240 int NumMyElements =
x_map->NumMyElements();
241 int *MyGlobalElements =
x_map->MyGlobalElements();
242 graph = Teuchos::rcp(
new Epetra_CrsGraph(Copy, *
x_map, 5));
243 for (
int i=0; i<NumMyElements; ++i ) {
246 int global_idx = MyGlobalElements[i];
247 graph->InsertGlobalIndices(global_idx, 1, &global_idx);
249 if (!
mesh[global_idx].boundary) {
252 graph->InsertGlobalIndices(global_idx, 1, &
mesh[global_idx].down);
256 graph->InsertGlobalIndices(global_idx, 1, &
mesh[global_idx].left);
260 graph->InsertGlobalIndices(global_idx, 1, &
mesh[global_idx].right);
264 graph->InsertGlobalIndices(global_idx, 1, &
mesh[global_idx].up);
267 graph->FillComplete();
268 graph->OptimizeStorage();
273 if (
basis == Teuchos::null) {
277 Normalized_KL_Diffusion_Func<KL_Diffusion_Func> nklFunc(
klFunc, *
basis);
283 int sz =
basis->size();
284 Teuchos::RCP<const Stokhos::ProductBasis<int, double> > prodbasis =
292 A = Teuchos::rcp(
new Epetra_CrsMatrix(Copy, *
graph));
295 b = Teuchos::rcp(
new Epetra_Vector(*
x_map));
296 for(
int i=0 ; i<NumMyElements; ++i ) {
297 int global_idx = MyGlobalElements[i];
298 if (
mesh[global_idx].boundary)
304 if (
basis != Teuchos::null) {
312 Teuchos::RCP<const Epetra_Map>
319 Teuchos::RCP<const Epetra_Map>
326 Teuchos::RCP<const Epetra_Map>
330 TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
333 "Error! twoD_diffusion_problem::get_p_map(): " <<
334 "Invalid parameter index l = " << l << std::endl);
339 Teuchos::RCP<const Epetra_Map>
343 TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
346 "Error! twoD_diffusion_problem::get_g_map(): " <<
347 "Invalid parameter index l = " << l << std::endl);
352 Teuchos::RCP<const Teuchos::Array<std::string> >
356 TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
359 "Error! twoD_diffusion_problem::get_p_names(): " <<
360 "Invalid parameter index l = " << l << std::endl);
365 Teuchos::RCP<const Epetra_Vector>
372 Teuchos::RCP<const Epetra_Vector>
376 TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
379 "Error! twoD_diffusion_problem::get_p_init(): " <<
380 "Invalid parameter index l = " << l << std::endl);
385 Teuchos::RCP<Epetra_Operator>
389 Teuchos::RCP<Epetra_CrsMatrix> AA =
390 Teuchos::rcp(
new Epetra_CrsMatrix(Copy, *
graph));
392 AA->OptimizeStorage();
399 const Epetra_Vector& p,
405 f.Update(-1.0, *
b, 1.0);
411 const Epetra_Vector& p,
416 Epetra_CrsMatrix& jac =
dynamic_cast<Epetra_CrsMatrix&
>(J);
419 jac.OptimizeStorage();
425 const Epetra_Vector& p,
443 const Teuchos::Array<double>& norms =
basis->norm_squared();
447 for (
int i=0;i<
basis->size();i++) {
453 Cijk_type::k_iterator k_begin = Cijk->k_begin();
454 Cijk_type::k_iterator k_end = Cijk->k_end();
455 for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
456 int k = Stokhos::index(k_it);
457 for (Cijk_type::kj_iterator j_it = Cijk->j_begin(k_it);
458 j_it != Cijk->j_end(k_it); ++j_it) {
459 int j = Stokhos::index(j_it);
462 for (Cijk_type::kj_iterator j_it = Cijk->j_begin(k_it);
463 j_it != Cijk->j_end(k_it); ++j_it) {
464 int j = Stokhos::index(j_it);
465 for (Cijk_type::kji_iterator i_it = Cijk->i_begin(j_it);
466 i_it != Cijk->i_end(j_it); ++i_it) {
467 int i = Stokhos::index(i_it);
468 double c = Stokhos::value(i_it);
473 f_sg[0].Update(-1.0,*
b,1.0);
482 for (
int i=0; i<J_sg.
size(); i++) {
483 Teuchos::RCP<Epetra_CrsMatrix> jac =
484 Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(J_sg.
getCoeffPtr(i),
true);
487 jac->OptimizeStorage();
497 int sz = x_sg.
size();
498 for (
int i=0; i<sz; i++) {
499 x_sg[i].MeanValue(&g_sg[i][0]);
504 template <
typename FuncT>
509 int NumMyElements =
x_map->NumMyElements();
510 int *MyGlobalElements =
x_map->MyGlobalElements();
515 for (
int k=0; k<sz; k++) {
516 A_k[k] = Teuchos::rcp(
new Epetra_CrsMatrix(Copy, *
graph));
517 for(
int i=0 ; i<NumMyElements; ++i ) {
520 int global_idx = MyGlobalElements[i];
521 if (
mesh[global_idx].boundary) {
526 A_k[k]->ReplaceGlobalValues(global_idx, 1, &
val, &global_idx);
530 -func(
mesh[global_idx].
x,
mesh[global_idx].
y-
h/2.0, k)/h2;
532 -func(
mesh[global_idx].
x-
h/2.0,
mesh[global_idx].
y, k)/h2;
534 -func(
mesh[global_idx].
x+
h/2.0,
mesh[global_idx].
y, k)/h2;
536 -func(
mesh[global_idx].
x,
mesh[global_idx].
y+
h/2.0, k)/h2;
539 val = -(a_down + a_left + a_right + a_up);
540 A_k[k]->ReplaceGlobalValues(global_idx, 1, &
val, &global_idx);
544 A_k[k]->ReplaceGlobalValues(global_idx, 1, &a_down,
545 &
mesh[global_idx].down);
549 A_k[k]->ReplaceGlobalValues(global_idx, 1, &a_left,
550 &
mesh[global_idx].left);
554 A_k[k]->ReplaceGlobalValues(global_idx, 1, &a_right,
555 &
mesh[global_idx].right);
559 A_k[k]->ReplaceGlobalValues(global_idx, 1, &a_up,
560 &
mesh[global_idx].up);
563 A_k[k]->FillComplete();
564 A_k[k]->OptimizeStorage();
572 if (
basis != Teuchos::null) {
573 for (
int i=0; i<
point.size(); i++)
577 for (
int k=0;k<
A_k.size();k++)
578 EpetraExt::MatrixMatrix::Add((*
A_k[k]),
false,
basis_vals[k], *
A, 1.0);
582 for (
int k=1;k<
A_k.size();k++)
583 EpetraExt::MatrixMatrix::Add((*
A_k[k]),
false, p[k-1], *
A, 1.0);
586 A->OptimizeStorage();
ordinal_type size() const
Return size.
A container class storing an orthogonal polynomial whose coefficients are vectors, operators, or in general any type that would have an expensive copy constructor.
void computeSGJacobian(const Stokhos::EpetraVectorOrthogPoly &x_sg, const Stokhos::EpetraVectorOrthogPoly &p_sg, Stokhos::EpetraOperatorOrthogPoly &J_sg)
Compute Jacobian.
Teuchos::Array< Teuchos::RCP< Epetra_CrsMatrix > > A_k
KL coefficients of operator.
Teuchos::RCP< Epetra_Map > x_map
Solution vector map.
KOKKOS_INLINE_FUNCTION PCE< Storage > pow(const PCE< Storage > &a, const PCE< Storage > &b)
void init(const value_type &val)
Initialize coefficients.
Teuchos::RCP< const Epetra_Vector > get_p_init(int l) const
Return initial parameters.
Teuchos::RCP< Epetra_CrsMatrix > A
Matrix to store deterministic operator.
void computeResidual(const Epetra_Vector &x, const Epetra_Vector &p, Epetra_Vector &f)
Compute residual.
void computeResponse(const Epetra_Vector &x, const Epetra_Vector &p, Epetra_Vector &g)
Compute response.
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< LogNormal_Diffusion_Func< KL_Diffusion_Func > > lnFunc
void computeSGResponse(const Stokhos::EpetraVectorOrthogPoly &x_sg, const Stokhos::EpetraVectorOrthogPoly &p_sg, Stokhos::EpetraVectorOrthogPoly &g_sg)
Compute SG response.
Teuchos::RCP< Epetra_Map > g_map
Response vector map.
Teuchos::RCP< Epetra_Map > p_map
Parameter vector map.
void computeJacobian(const Epetra_Vector &x, const Epetra_Vector &p, Epetra_Operator &J)
Compute Jacobian.
Teuchos::RCP< const Epetra_Map > get_p_map(int l) const
Return parameter vector map.
Teuchos::RCP< Epetra_Vector > p_init
Initial parameters.
Abstract base class for orthogonal polynomial-based expansions.
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int l) const
Return array of parameter names.
Teuchos::RCP< const Epetra_Vector > get_x_init() const
Return initial solution.
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
Teuchos::Array< MeshPoint > mesh
A container class storing an orthogonal polynomial whose coefficients are vectors, operators, or in general any type that would have an expensive copy constructor.
Teuchos::RCP< KL_Diffusion_Func > klFunc
Teuchos::Array< double > basis_vals
Array to store values of basis at a point.
Teuchos::RCP< Epetra_Operator > create_W() const
Create W = alpha*M + beta*J matrix.
Teuchos::RCP< Epetra_Vector > x_init
Initial guess.
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.
Abstract base class for multivariate orthogonal polynomials generated from tensor products of univari...
twoD_diffusion_problem(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)
Constructor.
KOKKOS_INLINE_FUNCTION PCE< Storage > exp(const PCE< Storage > &a)
Stokhos::Sparse3Tensor< int, double > Cijk_type
Teuchos::RCP< coeff_type > getCoeffPtr(ordinal_type i)
Return ref-count pointer to coefficient i.
Teuchos::RCP< const Epetra_Map > get_g_map(int j) const
Return response function map.
ScalarType f(const Teuchos::Array< ScalarType > &x, double a, double b)
Teuchos::Array< double > point
Array to store a point for basis evaluation.
Class representing a KL expansion of an exponential random field.
Teuchos::RCP< Epetra_Vector > b
Deterministic RHS.
Teuchos::RCP< const Epetra_Map > get_f_map() const
Return residual vector map.
void compute_A(const Epetra_Vector &p)
Compute A matrix.
void fillMatrices(const FuncT &func, int sz)
Fill coefficient matrix given function to evaluate diffusion coefficient.
Teuchos::RCP< Teuchos::Array< std::string > > p_names
Parameter names.
Teuchos::RCP< Epetra_CrsGraph > graph
Jacobian graph.
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType ValueType * y
Teuchos::RCP< const Epetra_Map > get_x_map() const
Return solution vector map.
ScalarType g(const Teuchos::Array< ScalarType > &x, const ScalarType &y)
void computeSGResidual(const Stokhos::EpetraVectorOrthogPoly &x_sg, const Stokhos::EpetraVectorOrthogPoly &p_sg, Stokhos::OrthogPolyExpansion< int, double > &expn, Stokhos::EpetraVectorOrthogPoly &f_sg)
Compute SG residual.
virtual Teuchos::RCP< const Sparse3Tensor< ordinal_type, value_type > > getTripleProduct() const =0
Get triple product.