42 template <
typename Scalar,
typename MeshScalar,
typename BasisScalar,
48 LocalOrdinal n, LocalOrdinal d,
49 BasisScalar s, BasisScalar mu,
51 bool eliminate_bcs_) :
53 log_normal(log_normal_),
54 eliminate_bcs(eliminate_bcs_)
58 using Teuchos::arrayView;
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;
77 for (GlobalOrdinal
j=0;
j<n;
j++) {
78 MeshScalar y = xyLeft +
j*
h;
79 for (GlobalOrdinal i=0; i<n; i++) {
80 MeshScalar x = xyLeft + i*
h;
81 GlobalOrdinal idx =
j*n+i;
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);
112 x_init->putScalar(0.0);
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);
122 p_init->putScalar(0.0);
126 for (LocalOrdinal i=0;i<d;i++) {
127 std::stringstream ss;
128 ss <<
"KL Random Variable " << i+1;
129 (*p_names)[i] = ss.str();
133 size_t NumMyElements =
x_map->getNodeNumElements();
134 ArrayView<const GlobalOrdinal> MyGlobalElements =
135 x_map->getNodeElementList ();
137 for (
size_t i=0; i<NumMyElements; ++i ) {
140 GlobalOrdinal global_idx = MyGlobalElements[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) {
174 GlobalOrdinal global_idx = MyGlobalElements[i];
175 if (
mesh[global_idx].boundary)
187 template <
typename Scalar,
typename MeshScalar,
typename BasisScalar,
200 template <
typename Scalar,
typename MeshScalar,
typename BasisScalar,
213 template <
typename Scalar,
typename MeshScalar,
typename BasisScalar,
221 get_p_map(LocalOrdinal l)
const
226 "Error! twoD_diffusion_problem::get_p_map(): " <<
227 "Invalid parameter index l = " << l << std::endl);
232 template <
typename Scalar,
typename MeshScalar,
typename BasisScalar,
240 get_g_map(LocalOrdinal l)
const
245 "Error! twoD_diffusion_problem::get_g_map(): " <<
246 "Invalid parameter index l = " << l << std::endl);
251 template <
typename Scalar,
typename MeshScalar,
typename BasisScalar,
256 get_p_names(LocalOrdinal l)
const
261 "Error! twoD_diffusion_problem::get_p_names(): " <<
262 "Invalid parameter index l = " << l << std::endl);
267 template <
typename Scalar,
typename MeshScalar,
typename BasisScalar,
280 template <
typename Scalar,
typename MeshScalar,
typename BasisScalar,
288 get_p_init(LocalOrdinal l)
const
293 "Error! twoD_diffusion_problem::get_p_init(): " <<
294 "Invalid parameter index l = " << l << std::endl);
299 template <
typename Scalar,
typename MeshScalar,
typename BasisScalar,
314 template <
typename Scalar,
typename MeshScalar,
typename BasisScalar,
325 computeA(*lnFunc, p, *
A);
327 computeA(*klFunc, p, *
A);
329 f.update(-1.0, *b, 1.0);
332 template <
typename Scalar,
typename MeshScalar,
typename BasisScalar,
342 computeA(*lnFunc, p, jac);
344 computeA(*klFunc, p, jac);
347 template <
typename Scalar,
typename MeshScalar,
typename BasisScalar,
358 x.meanValue(g_view());
359 g_view[0] *=
Scalar(x.getGlobalLength()) /
Scalar(mesh.size());
362 template <
typename Scalar,
typename MeshScalar,
typename BasisScalar,
364 template <
typename FuncT>
371 using Teuchos::arrayView;
374 jac.setAllToScalar(0.0);
378 size_t NumMyElements = x_map->getNodeNumElements();
379 ArrayView<const GlobalOrdinal> MyGlobalElements =
380 x_map->getNodeElementList ();
384 for(
size_t i=0 ; i<NumMyElements; ++i ) {
387 GlobalOrdinal global_idx = MyGlobalElements[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),
436 template <
typename Scalar,
typename MeshScalar,
typename BasisScalar,
439 Node>::KL_Diffusion_Func::
440 KL_Diffusion_Func(MeshScalar xyLeft, MeshScalar xyRight,
441 BasisScalar mean, BasisScalar std_dev,
442 MeshScalar L, LocalOrdinal num_KL) : point(2)
445 rfParams.
set(
"Number of KL Terms", num_KL);
446 rfParams.
set(
"Mean", mean);
447 rfParams.
set(
"Standard Deviation", std_dev);
448 LocalOrdinal ndim = 2;
450 correlation_length(ndim);
451 for (LocalOrdinal i=0; i<ndim; i++) {
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);
464 template <
typename Scalar,
typename MeshScalar,
typename BasisScalar,
468 Node>::KL_Diffusion_Func::
473 return rf->evaluate(
point, rv);
Teuchos::RCP< Epetra_Map > x_map
Solution vector map.
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Teuchos::RCP< LogNormal_Diffusion_Func< KL_Diffusion_Func > > lnFunc
Teuchos::RCP< Epetra_Map > g_map
Response vector map.
Teuchos::RCP< Epetra_Map > p_map
Parameter vector map.
Teuchos::RCP< Epetra_Vector > p_init
Initial parameters.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::Array< MeshPoint > mesh
Teuchos::RCP< Stokhos::KL::ExponentialRandomField< MeshScalar > > rf
Teuchos::RCP< KL_Diffusion_Func > klFunc
Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > Tpetra_CrsGraph
KokkosClassic::DefaultNode::DefaultNodeType Node
Teuchos::RCP< Epetra_Vector > x_init
Initial guess.
Tpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > Tpetra_Vector
ScalarType f(const Teuchos::Array< ScalarType > &x, double a, double b)
Teuchos::Array< double > point
Array to store a point for basis evaluation.
Teuchos::RCP< Epetra_Vector > b
Deterministic RHS.
Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > Tpetra_CrsMatrix
Teuchos::RCP< Teuchos::Array< std::string > > p_names
Parameter names.
Teuchos::RCP< Epetra_CrsGraph > graph
Jacobian graph.
ScalarType g(const Teuchos::Array< ScalarType > &x, const ScalarType &y)
A linear 2-D diffusion problem.