42 template <
typename Scalar,
typename MeshScalar,
typename BasisScalar,
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;
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);
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 ) {
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)
187 template <
typename Scalar,
typename MeshScalar,
typename BasisScalar,
200 template <
typename Scalar,
typename MeshScalar,
typename BasisScalar,
213 template <
typename Scalar,
typename MeshScalar,
typename BasisScalar,
226 "Error! twoD_diffusion_problem::get_p_map(): " <<
227 "Invalid parameter index l = " << l << std::endl);
232 template <
typename Scalar,
typename MeshScalar,
typename BasisScalar,
245 "Error! twoD_diffusion_problem::get_g_map(): " <<
246 "Invalid parameter index l = " << l << std::endl);
251 template <
typename Scalar,
typename MeshScalar,
typename BasisScalar,
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,
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 ) {
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,
445 rfParams.
set(
"Number of KL Terms", num_KL);
446 rfParams.
set(
"Mean", mean);
447 rfParams.
set(
"Standard Deviation", std_dev);
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);
464 template <
typename Scalar,
typename MeshScalar,
typename BasisScalar,
468 Node>::KL_Diffusion_Func::
473 return rf->evaluate(
point, rv);