42template <
typename Scalar,
typename MeshScalar,
typename BasisScalar,
46twoD_diffusion_problem(
47 const Teuchos::RCP<
const Teuchos::Comm<int> >& comm,
49 BasisScalar s, BasisScalar mu,
51 bool eliminate_bcs_) :
53 log_normal(log_normal_),
54 eliminate_bcs(eliminate_bcs_)
57 using Teuchos::ArrayView;
58 using Teuchos::arrayView;
59 using Teuchos::ArrayRCP;
62 using Tpetra::global_size_t;
73 MeshScalar xyLeft = -.5;
74 MeshScalar xyRight = .5;
75 h = (xyRight - xyLeft)/((MeshScalar)(n-1));
76 Array<GlobalOrdinal> global_dof_indices;
78 MeshScalar y = xyLeft +
j*
h;
80 MeshScalar x = xyLeft + i*
h;
84 if (i == 0 || i == n-1 ||
j == 0 ||
j == n-1)
85 mesh[idx].boundary =
true;
87 mesh[idx].left = idx-1;
89 mesh[idx].right = idx+1;
91 mesh[idx].down = idx-n;
95 global_dof_indices.push_back(idx);
100 global_size_t n_global_dof = global_dof_indices.size();
101 int n_proc = comm->getSize();
102 int proc_id = comm->getRank();
103 size_t n_my_dof = n_global_dof / n_proc;
104 if (proc_id == n_proc-1)
105 n_my_dof += n_global_dof % n_proc;
106 ArrayView<GlobalOrdinal> my_dof =
107 global_dof_indices.view(proc_id*(n_global_dof / n_proc), n_my_dof);
108 x_map = Tpetra::createNonContigMap<LocalOrdinal,GlobalOrdinal>(my_dof, comm);
111 x_init = Tpetra::createVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(
x_map);
115 p_map = Tpetra::createLocalMap<LocalOrdinal,GlobalOrdinal>(d, comm);
118 g_map = Tpetra::createLocalMap<LocalOrdinal,GlobalOrdinal>(1, comm);
121 p_init = Tpetra::createVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(
p_map);
125 p_names = Teuchos::rcp(
new Array<std::string>(d));
127 std::stringstream ss;
128 ss <<
"KL Random Variable " << i+1;
129 (*p_names)[i] = ss.str();
133 size_t NumMyElements =
x_map->getLocalNumElements();
134 ArrayView<const GlobalOrdinal> MyGlobalElements =
135 x_map->getLocalElementList ();
137 for (
size_t i=0; i<NumMyElements; ++i ) {
141 graph->insertGlobalIndices(global_idx, arrayView(&global_idx, 1));
143 if (!
mesh[global_idx].boundary) {
146 graph->insertGlobalIndices(global_idx,
147 arrayView(&
mesh[global_idx].down,1));
151 graph->insertGlobalIndices(global_idx,
152 arrayView(&
mesh[global_idx].left,1));
156 graph->insertGlobalIndices(global_idx,
157 arrayView(&
mesh[global_idx].right,1));
161 graph->insertGlobalIndices(global_idx,
162 arrayView(&
mesh[global_idx].up,1));
165 graph->fillComplete();
171 b = Tpetra::createVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(
x_map);
172 ArrayRCP<Scalar> b_view =
b->get1dViewNonConst();
173 for(
size_t i=0; i<NumMyElements; ++i) {
175 if (
mesh[global_idx].boundary)
187template <
typename Scalar,
typename MeshScalar,
typename BasisScalar,
200template <
typename Scalar,
typename MeshScalar,
typename BasisScalar,
213template <
typename Scalar,
typename MeshScalar,
typename BasisScalar,
223 TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
226 "Error! twoD_diffusion_problem::get_p_map(): " <<
227 "Invalid parameter index l = " << l << std::endl);
232template <
typename Scalar,
typename MeshScalar,
typename BasisScalar,
242 TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
245 "Error! twoD_diffusion_problem::get_g_map(): " <<
246 "Invalid parameter index l = " << l << std::endl);
251template <
typename Scalar,
typename MeshScalar,
typename BasisScalar,
253Teuchos::RCP<const Teuchos::Array<std::string> >
258 TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
261 "Error! twoD_diffusion_problem::get_p_names(): " <<
262 "Invalid parameter index l = " << l << std::endl);
267template <
typename Scalar,
typename MeshScalar,
typename BasisScalar,
280template <
typename Scalar,
typename MeshScalar,
typename BasisScalar,
290 TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
293 "Error! twoD_diffusion_problem::get_p_init(): " <<
294 "Invalid parameter index l = " << l << std::endl);
299template <
typename Scalar,
typename MeshScalar,
typename BasisScalar,
308 Teuchos::RCP<Tpetra_CrsMatrix> AA =
309 Teuchos::rcp(
new Tpetra_CrsMatrix(graph));
314template <
typename Scalar,
typename MeshScalar,
typename BasisScalar,
325 computeA(*lnFunc, p, *A);
327 computeA(*klFunc, p, *A);
329 f.update(-1.0, *b, 1.0);
332template <
typename Scalar,
typename MeshScalar,
typename BasisScalar,
342 computeA(*lnFunc, p, jac);
344 computeA(*klFunc, p, jac);
347template <
typename Scalar,
typename MeshScalar,
typename BasisScalar,
357 Teuchos::ArrayRCP<Scalar> g_view =
g.get1dViewNonConst();
358 x.meanValue(g_view());
359 g_view[0] *=
Scalar(x.getGlobalLength()) /
Scalar(mesh.size());
362template <
typename Scalar,
typename MeshScalar,
typename BasisScalar,
364template <
typename FuncT>
370 using Teuchos::ArrayView;
371 using Teuchos::arrayView;
374 jac.setAllToScalar(0.0);
376 Teuchos::ArrayRCP<const Scalar> p_view = p.get1dView();
377 Teuchos::Array<Scalar> rv(p_view());
378 size_t NumMyElements = x_map->getLocalNumElements();
379 ArrayView<const GlobalOrdinal> MyGlobalElements =
380 x_map->getLocalElementList ();
384 for(
size_t i=0 ; i<NumMyElements; ++i ) {
388 if (mesh[global_idx].boundary) {
390 jac.replaceGlobalValues(global_idx, arrayView(&global_idx,1),
395 -func(mesh[global_idx].x, mesh[global_idx].y-h/2.0, rv)/h2;
397 -func(mesh[global_idx].x-h/2.0, mesh[global_idx].y, rv)/h2;
399 -func(mesh[global_idx].x+h/2.0, mesh[global_idx].y, rv)/h2;
401 -func(mesh[global_idx].x, mesh[global_idx].y+h/2.0, rv)/h2;
404 val = -(a_down + a_left + a_right + a_up);
405 jac.replaceGlobalValues(global_idx, arrayView(&global_idx,1),
409 if (!(eliminate_bcs && mesh[mesh[global_idx].down].boundary))
410 jac.replaceGlobalValues(global_idx,
411 arrayView(&mesh[global_idx].down,1),
412 arrayView(&a_down,1));
415 if (!(eliminate_bcs && mesh[mesh[global_idx].left].boundary))
416 jac.replaceGlobalValues(global_idx,
417 arrayView(&mesh[global_idx].left,1),
418 arrayView(&a_left,1));
421 if (!(eliminate_bcs && mesh[mesh[global_idx].right].boundary))
422 jac.replaceGlobalValues(global_idx,
423 arrayView(&mesh[global_idx].right,1),
424 arrayView(&a_right,1));
427 if (!(eliminate_bcs && mesh[mesh[global_idx].up].boundary))
428 jac.replaceGlobalValues(global_idx,
429 arrayView(&mesh[global_idx].up,1),
436template <
typename Scalar,
typename MeshScalar,
typename BasisScalar,
440KL_Diffusion_Func(MeshScalar xyLeft, MeshScalar xyRight,
441 BasisScalar mean, BasisScalar std_dev,
444 Teuchos::ParameterList rfParams;
445 rfParams.set(
"Number of KL Terms", num_KL);
446 rfParams.set(
"Mean", mean);
447 rfParams.set(
"Standard Deviation", std_dev);
449 Teuchos::Array<MeshScalar> domain_upper(ndim), domain_lower(ndim),
450 correlation_length(ndim);
452 domain_upper[i] = xyRight;
453 domain_lower[i] = xyLeft;
454 correlation_length[i] = L;
456 rfParams.set(
"Domain Upper Bounds", domain_upper);
457 rfParams.set(
"Domain Lower Bounds", domain_lower);
458 rfParams.set(
"Correlation Lengths", correlation_length);
464template <
typename Scalar,
typename MeshScalar,
typename BasisScalar,
469operator() (MeshScalar x, MeshScalar y,
const Teuchos::Array<Scalar>& rv)
const
473 return rf->evaluate(
point, rv);
Class representing a KL expansion of an exponential random field.
A linear 2-D diffusion problem.
Teuchos::RCP< Epetra_Vector > p_init
Initial parameters.
Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > Tpetra_CrsMatrix
Teuchos::RCP< Epetra_Map > x_map
Solution vector map.
Teuchos::RCP< Epetra_Vector > x_init
Initial guess.
Teuchos::RCP< Epetra_CrsMatrix > A
Matrix to store deterministic operator.
Teuchos::RCP< Teuchos::Array< std::string > > p_names
Parameter names.
Teuchos::RCP< Epetra_Map > p_map
Parameter vector map.
Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > Tpetra_CrsGraph
Teuchos::RCP< LogNormal_Diffusion_Func< KL_Diffusion_Func > > lnFunc
Teuchos::Array< MeshPoint > mesh
Teuchos::Array< double > point
Array to store a point for basis evaluation.
Tpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > Tpetra_Vector
Teuchos::RCP< KL_Diffusion_Func > klFunc
Teuchos::RCP< Epetra_Vector > b
Deterministic RHS.
Teuchos::RCP< Epetra_CrsGraph > graph
Jacobian graph.
Teuchos::RCP< Epetra_Map > g_map
Response vector map.
KokkosClassic::DefaultNode::DefaultNodeType Node
ScalarType g(const Teuchos::Array< ScalarType > &x, const ScalarType &y)
ScalarType f(const Teuchos::Array< ScalarType > &x, double a, double b)
Teuchos::RCP< Stokhos::KL::ExponentialRandomField< MeshScalar > > rf