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)
183 lnFunc = rcp(
new LogNormal_Diffusion_Func<KL_Diffusion_Func>(*
klFunc));
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),
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);