Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
twoD_diffusion_problem_tpetra_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Stokhos Package
4 //
5 // Copyright 2009 NTESS and the Stokhos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 template <typename Scalar, typename MeshScalar, typename BasisScalar,
11  typename LocalOrdinal, typename GlobalOrdinal, typename Node>
12 twoD_diffusion_problem<Scalar,MeshScalar,BasisScalar,LocalOrdinal,GlobalOrdinal,
13  Node>::
14 twoD_diffusion_problem(
15  const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
16  LocalOrdinal n, LocalOrdinal d,
17  BasisScalar s, BasisScalar mu,
18  bool log_normal_,
19  bool eliminate_bcs_) :
20  mesh(n*n),
21  log_normal(log_normal_),
22  eliminate_bcs(eliminate_bcs_)
23 {
24  using Teuchos::Array;
25  using Teuchos::ArrayView;
26  using Teuchos::arrayView;
27  using Teuchos::ArrayRCP;
28  using Teuchos::RCP;
29  using Teuchos::rcp;
30  using Tpetra::global_size_t;
31 
33  // Construct the mesh.
34  // The mesh is uniform and the nodes are numbered
35  // LEFT to RIGHT, DOWN to UP.
36  //
37  // 5-6-7-8-9
38  // | | | | |
39  // 0-1-2-3-4
41  MeshScalar xyLeft = -.5;
42  MeshScalar xyRight = .5;
43  h = (xyRight - xyLeft)/((MeshScalar)(n-1));
44  Array<GlobalOrdinal> global_dof_indices;
45  for (GlobalOrdinal j=0; j<n; j++) {
46  MeshScalar y = xyLeft + j*h;
47  for (GlobalOrdinal i=0; i<n; i++) {
48  MeshScalar x = xyLeft + i*h;
49  GlobalOrdinal idx = j*n+i;
50  mesh[idx].x = x;
51  mesh[idx].y = y;
52  if (i == 0 || i == n-1 || j == 0 || j == n-1)
53  mesh[idx].boundary = true;
54  if (i != 0)
55  mesh[idx].left = idx-1;
56  if (i != n-1)
57  mesh[idx].right = idx+1;
58  if (j != 0)
59  mesh[idx].down = idx-n;
60  if (j != n-1)
61  mesh[idx].up = idx+n;
62  if (!(eliminate_bcs && mesh[idx].boundary))
63  global_dof_indices.push_back(idx);
64  }
65  }
66 
67  // Solution vector map
68  global_size_t n_global_dof = global_dof_indices.size();
69  int n_proc = comm->getSize();
70  int proc_id = comm->getRank();
71  size_t n_my_dof = n_global_dof / n_proc;
72  if (proc_id == n_proc-1)
73  n_my_dof += n_global_dof % n_proc;
74  ArrayView<GlobalOrdinal> my_dof =
75  global_dof_indices.view(proc_id*(n_global_dof / n_proc), n_my_dof);
76  x_map = Tpetra::createNonContigMap<LocalOrdinal,GlobalOrdinal>(my_dof, comm);
77 
78  // Initial guess, initialized to 0.0
79  x_init = Tpetra::createVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(x_map);
80  x_init->putScalar(0.0);
81 
82  // Parameter vector map
83  p_map = Tpetra::createLocalMap<LocalOrdinal,GlobalOrdinal>(d, comm);
84 
85  // Response vector map
86  g_map = Tpetra::createLocalMap<LocalOrdinal,GlobalOrdinal>(1, comm);
87 
88  // Initial parameters
89  p_init = Tpetra::createVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(p_map);
90  p_init->putScalar(0.0);
91 
92  // Parameter names
93  p_names = Teuchos::rcp(new Array<std::string>(d));
94  for (LocalOrdinal i=0;i<d;i++) {
95  std::stringstream ss;
96  ss << "KL Random Variable " << i+1;
97  (*p_names)[i] = ss.str();
98  }
99 
100  // Build Jacobian graph
101  size_t NumMyElements = x_map->getLocalNumElements();
102  ArrayView<const GlobalOrdinal> MyGlobalElements =
103  x_map->getLocalElementList ();
104  graph = rcp(new Tpetra_CrsGraph(x_map, 5));
105  for (size_t i=0; i<NumMyElements; ++i ) {
106 
107  // Center
108  GlobalOrdinal global_idx = MyGlobalElements[i];
109  graph->insertGlobalIndices(global_idx, arrayView(&global_idx, 1));
110 
111  if (!mesh[global_idx].boundary) {
112  // Down
113  if (!(eliminate_bcs && mesh[mesh[global_idx].down].boundary))
114  graph->insertGlobalIndices(global_idx,
115  arrayView(&mesh[global_idx].down,1));
116 
117  // Left
118  if (!(eliminate_bcs && mesh[mesh[global_idx].left].boundary))
119  graph->insertGlobalIndices(global_idx,
120  arrayView(&mesh[global_idx].left,1));
121 
122  // Right
123  if (!(eliminate_bcs && mesh[mesh[global_idx].right].boundary))
124  graph->insertGlobalIndices(global_idx,
125  arrayView(&mesh[global_idx].right,1));
126 
127  // Up
128  if (!(eliminate_bcs && mesh[mesh[global_idx].up].boundary))
129  graph->insertGlobalIndices(global_idx,
130  arrayView(&mesh[global_idx].up,1));
131  }
132  }
133  graph->fillComplete();
134 
135  // Construct deterministic operator
136  A = rcp(new Tpetra_CrsMatrix(graph));
137 
138  // Construct the RHS vector.
139  b = Tpetra::createVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(x_map);
140  ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
141  for(size_t i=0; i<NumMyElements; ++i) {
142  GlobalOrdinal global_idx = MyGlobalElements[i];
143  if (mesh[global_idx].boundary)
144  b_view[i] = 0;
145  else
146  b_view[i] = 1;
147  }
148 
149  // Diffusion functions
150  klFunc = rcp(new KL_Diffusion_Func(xyLeft, xyRight, mu, s, 1.0, d));
151  lnFunc = rcp(new LogNormal_Diffusion_Func<KL_Diffusion_Func>(*klFunc));
152 }
153 
154 // Overridden from EpetraExt::ModelEvaluator
155 template <typename Scalar, typename MeshScalar, typename BasisScalar,
156  typename LocalOrdinal, typename GlobalOrdinal, typename Node>
157 Teuchos::RCP<const
158  typename twoD_diffusion_problem<Scalar,MeshScalar,BasisScalar,
159  LocalOrdinal,GlobalOrdinal,Node
160  >::Tpetra_Map>
161 twoD_diffusion_problem<Scalar,MeshScalar,BasisScalar,LocalOrdinal,GlobalOrdinal,
162  Node>::
163 get_x_map() const
164 {
165  return x_map;
166 }
167 
168 template <typename Scalar, typename MeshScalar, typename BasisScalar,
169  typename LocalOrdinal, typename GlobalOrdinal, typename Node>
170 Teuchos::RCP<const
171  typename twoD_diffusion_problem<Scalar,MeshScalar,BasisScalar,
172  LocalOrdinal,GlobalOrdinal,Node
173  >::Tpetra_Map>
174 twoD_diffusion_problem<Scalar,MeshScalar,BasisScalar,LocalOrdinal,GlobalOrdinal,
175  Node>::
176 get_f_map() const
177 {
178  return x_map;
179 }
180 
181 template <typename Scalar, typename MeshScalar, typename BasisScalar,
182  typename LocalOrdinal, typename GlobalOrdinal, typename Node>
183 Teuchos::RCP<const
184  typename twoD_diffusion_problem<Scalar,MeshScalar,BasisScalar,
185  LocalOrdinal,GlobalOrdinal,Node
186  >::Tpetra_Map>
187 twoD_diffusion_problem<Scalar,MeshScalar,BasisScalar,LocalOrdinal,GlobalOrdinal,
188  Node>::
189 get_p_map(LocalOrdinal l) const
190 {
191  TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
192  std::logic_error,
193  std::endl <<
194  "Error! twoD_diffusion_problem::get_p_map(): " <<
195  "Invalid parameter index l = " << l << std::endl);
196 
197  return p_map;
198 }
199 
200 template <typename Scalar, typename MeshScalar, typename BasisScalar,
201  typename LocalOrdinal, typename GlobalOrdinal, typename Node>
202 Teuchos::RCP<const
203  typename twoD_diffusion_problem<Scalar,MeshScalar,BasisScalar,
204  LocalOrdinal,GlobalOrdinal,Node
205  >::Tpetra_Map>
206 twoD_diffusion_problem<Scalar,MeshScalar,BasisScalar,LocalOrdinal,GlobalOrdinal,
207  Node>::
208 get_g_map(LocalOrdinal l) const
209 {
210  TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
211  std::logic_error,
212  std::endl <<
213  "Error! twoD_diffusion_problem::get_g_map(): " <<
214  "Invalid parameter index l = " << l << std::endl);
215 
216  return g_map;
217 }
218 
219 template <typename Scalar, typename MeshScalar, typename BasisScalar,
220  typename LocalOrdinal, typename GlobalOrdinal, typename Node>
221 Teuchos::RCP<const Teuchos::Array<std::string> >
222 twoD_diffusion_problem<Scalar,MeshScalar,BasisScalar,LocalOrdinal,GlobalOrdinal,
223  Node>::
224 get_p_names(LocalOrdinal l) const
225 {
226  TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
227  std::logic_error,
228  std::endl <<
229  "Error! twoD_diffusion_problem::get_p_names(): " <<
230  "Invalid parameter index l = " << l << std::endl);
231 
232  return p_names;
233 }
234 
235 template <typename Scalar, typename MeshScalar, typename BasisScalar,
236  typename LocalOrdinal, typename GlobalOrdinal, typename Node>
237 Teuchos::RCP<const
238  typename twoD_diffusion_problem<Scalar,MeshScalar,BasisScalar,
239  LocalOrdinal,GlobalOrdinal,Node
240  >::Tpetra_Vector>
241 twoD_diffusion_problem<Scalar,MeshScalar,BasisScalar,LocalOrdinal,GlobalOrdinal,
242  Node>::
243 get_x_init() const
244 {
245  return x_init;
246 }
247 
248 template <typename Scalar, typename MeshScalar, typename BasisScalar,
249  typename LocalOrdinal, typename GlobalOrdinal, typename Node>
250 Teuchos::RCP<const
251  typename twoD_diffusion_problem<Scalar,MeshScalar,BasisScalar,
252  LocalOrdinal,GlobalOrdinal,Node
253  >::Tpetra_Vector>
254 twoD_diffusion_problem<Scalar,MeshScalar,BasisScalar,LocalOrdinal,GlobalOrdinal,
255  Node>::
256 get_p_init(LocalOrdinal l) const
257 {
258  TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
259  std::logic_error,
260  std::endl <<
261  "Error! twoD_diffusion_problem::get_p_init(): " <<
262  "Invalid parameter index l = " << l << std::endl);
263 
264  return p_init;
265 }
266 
267 template <typename Scalar, typename MeshScalar, typename BasisScalar,
268  typename LocalOrdinal, typename GlobalOrdinal, typename Node>
269 Teuchos::RCP<typename twoD_diffusion_problem<Scalar,MeshScalar,BasisScalar,
270  LocalOrdinal,GlobalOrdinal,Node
271  >::Tpetra_CrsMatrix>
272 twoD_diffusion_problem<Scalar,MeshScalar,BasisScalar,LocalOrdinal,GlobalOrdinal,
273  Node>::
274 create_W() const
275 {
276  Teuchos::RCP<Tpetra_CrsMatrix> AA =
277  Teuchos::rcp(new Tpetra_CrsMatrix(graph));
278  AA->fillComplete();
279  return AA;
280 }
281 
282 template <typename Scalar, typename MeshScalar, typename BasisScalar,
283  typename LocalOrdinal, typename GlobalOrdinal, typename Node>
284 void
285 twoD_diffusion_problem<Scalar,MeshScalar,BasisScalar,LocalOrdinal,GlobalOrdinal,
286  Node>::
287 computeResidual(const Tpetra_Vector& x,
288  const Tpetra_Vector& p,
289  Tpetra_Vector& f)
290 {
291  // f = A*x - b
292  if (log_normal)
293  computeA(*lnFunc, p, *A);
294  else
295  computeA(*klFunc, p, *A);
296  A->apply(x,f);
297  f.update(-1.0, *b, 1.0);
298 }
299 
300 template <typename Scalar, typename MeshScalar, typename BasisScalar,
301  typename LocalOrdinal, typename GlobalOrdinal, typename Node>
302 void
303 twoD_diffusion_problem<Scalar,MeshScalar,BasisScalar,LocalOrdinal,GlobalOrdinal,
304  Node>::
305 computeJacobian(const Tpetra_Vector& x,
306  const Tpetra_Vector& p,
307  Tpetra_CrsMatrix& jac)
308 {
309  if (log_normal)
310  computeA(*lnFunc, p, jac);
311  else
312  computeA(*klFunc, p, jac);
313 }
314 
315 template <typename Scalar, typename MeshScalar, typename BasisScalar,
316  typename LocalOrdinal, typename GlobalOrdinal, typename Node>
317 void
318 twoD_diffusion_problem<Scalar,MeshScalar,BasisScalar,LocalOrdinal,GlobalOrdinal,
319  Node>::
320 computeResponse(const Tpetra_Vector& x,
321  const Tpetra_Vector& p,
322  Tpetra_Vector& g)
323 {
324  // g = average of x
325  Teuchos::ArrayRCP<Scalar> g_view = g.get1dViewNonConst();
326  x.meanValue(g_view());
327  g_view[0] *= Scalar(x.getGlobalLength()) / Scalar(mesh.size());
328 }
329 
330 template <typename Scalar, typename MeshScalar, typename BasisScalar,
331  typename LocalOrdinal, typename GlobalOrdinal, typename Node>
332 template <typename FuncT>
333 void
334 twoD_diffusion_problem<Scalar,MeshScalar,BasisScalar,LocalOrdinal,GlobalOrdinal,
335  Node>::
336 computeA(const FuncT& func, const Tpetra_Vector& p, Tpetra_CrsMatrix& jac)
337 {
338  using Teuchos::ArrayView;
339  using Teuchos::arrayView;
340 
341  jac.resumeFill();
342  jac.setAllToScalar(0.0);
343 
344  Teuchos::ArrayRCP<const Scalar> p_view = p.get1dView();
345  Teuchos::Array<Scalar> rv(p_view());
346  size_t NumMyElements = x_map->getLocalNumElements();
347  ArrayView<const GlobalOrdinal> MyGlobalElements =
348  x_map->getLocalElementList ();
349  MeshScalar h2 = h*h;
350  Scalar val;
351 
352  for(size_t i=0 ; i<NumMyElements; ++i ) {
353 
354  // Center
355  GlobalOrdinal global_idx = MyGlobalElements[i];
356  if (mesh[global_idx].boundary) {
357  val = 1.0;
358  jac.replaceGlobalValues(global_idx, arrayView(&global_idx,1),
359  arrayView(&val,1));
360  }
361  else {
362  Scalar a_down =
363  -func(mesh[global_idx].x, mesh[global_idx].y-h/2.0, rv)/h2;
364  Scalar a_left =
365  -func(mesh[global_idx].x-h/2.0, mesh[global_idx].y, rv)/h2;
366  Scalar a_right =
367  -func(mesh[global_idx].x+h/2.0, mesh[global_idx].y, rv)/h2;
368  Scalar a_up =
369  -func(mesh[global_idx].x, mesh[global_idx].y+h/2.0, rv)/h2;
370 
371  // Center
372  val = -(a_down + a_left + a_right + a_up);
373  jac.replaceGlobalValues(global_idx, arrayView(&global_idx,1),
374  arrayView(&val,1));
375 
376  // Down
377  if (!(eliminate_bcs && mesh[mesh[global_idx].down].boundary))
378  jac.replaceGlobalValues(global_idx,
379  arrayView(&mesh[global_idx].down,1),
380  arrayView(&a_down,1));
381 
382  // Left
383  if (!(eliminate_bcs && mesh[mesh[global_idx].left].boundary))
384  jac.replaceGlobalValues(global_idx,
385  arrayView(&mesh[global_idx].left,1),
386  arrayView(&a_left,1));
387 
388  // Right
389  if (!(eliminate_bcs && mesh[mesh[global_idx].right].boundary))
390  jac.replaceGlobalValues(global_idx,
391  arrayView(&mesh[global_idx].right,1),
392  arrayView(&a_right,1));
393 
394  // Up
395  if (!(eliminate_bcs && mesh[mesh[global_idx].up].boundary))
396  jac.replaceGlobalValues(global_idx,
397  arrayView(&mesh[global_idx].up,1),
398  arrayView(&a_up,1));
399  }
400  }
401  jac.fillComplete();
402 }
403 
404 template <typename Scalar, typename MeshScalar, typename BasisScalar,
405  typename LocalOrdinal, typename GlobalOrdinal, typename Node>
406 twoD_diffusion_problem<Scalar,MeshScalar,BasisScalar,LocalOrdinal,GlobalOrdinal,
407  Node>::KL_Diffusion_Func::
408 KL_Diffusion_Func(MeshScalar xyLeft, MeshScalar xyRight,
409  BasisScalar mean, BasisScalar std_dev,
410  MeshScalar L, LocalOrdinal num_KL) : point(2)
411 {
412  Teuchos::ParameterList rfParams;
413  rfParams.set("Number of KL Terms", num_KL);
414  rfParams.set("Mean", mean);
415  rfParams.set("Standard Deviation", std_dev);
416  LocalOrdinal ndim = 2;
417  Teuchos::Array<MeshScalar> domain_upper(ndim), domain_lower(ndim),
418  correlation_length(ndim);
419  for (LocalOrdinal i=0; i<ndim; i++) {
420  domain_upper[i] = xyRight;
421  domain_lower[i] = xyLeft;
422  correlation_length[i] = L;
423  }
424  rfParams.set("Domain Upper Bounds", domain_upper);
425  rfParams.set("Domain Lower Bounds", domain_lower);
426  rfParams.set("Correlation Lengths", correlation_length);
427 
428  rf =
429  Teuchos::rcp(new Stokhos::KL::ExponentialRandomField<MeshScalar>(rfParams));
430 }
431 
432 template <typename Scalar, typename MeshScalar, typename BasisScalar,
433  typename LocalOrdinal, typename GlobalOrdinal, typename Node>
434 Scalar
435 twoD_diffusion_problem<Scalar,MeshScalar,BasisScalar,LocalOrdinal,GlobalOrdinal,
436  Node>::KL_Diffusion_Func::
437 operator() (MeshScalar x, MeshScalar y, const Teuchos::Array<Scalar>& rv) const
438 {
439  point[0] = x;
440  point[1] = y;
441  return rf->evaluate(point, rv);
442 }
twoD_diffusion_problem::x_map
Teuchos::RCP< Epetra_Map > x_map
Solution vector map.
Definition: twoD_diffusion_problem.hpp:128
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
twoD_diffusion_problem::LogNormal_Diffusion_Func
Definition: twoD_diffusion_problem_tpetra.hpp:163
Teuchos::ParameterList::set
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
twoD_diffusion_problem::lnFunc
Teuchos::RCP< LogNormal_Diffusion_Func< KL_Diffusion_Func > > lnFunc
Definition: twoD_diffusion_problem_tpetra.hpp:173
twoD_diffusion_problem::g_map
Teuchos::RCP< Epetra_Map > g_map
Response vector map.
Definition: twoD_diffusion_problem.hpp:140
twoD_diffusion_problem::p_map
Teuchos::RCP< Epetra_Map > p_map
Parameter vector map.
Definition: twoD_diffusion_problem.hpp:137
j
j
Definition: Sacado_Fad_Exp_MP_Vector.hpp:495
Scalar
pce_type Scalar
Definition: linear2d_diffusion_scalar_types.hpp:17
twoD_diffusion_problem::p_init
Teuchos::RCP< Epetra_Vector > p_init
Initial parameters.
Definition: twoD_diffusion_problem.hpp:143
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
twoD_diffusion_problem::mesh
Teuchos::Array< MeshPoint > mesh
Definition: twoD_diffusion_problem.hpp:120
twoD_diffusion_problem::KL_Diffusion_Func::rf
Teuchos::RCP< Stokhos::KL::ExponentialRandomField< MeshScalar > > rf
Definition: twoD_diffusion_problem_tpetra.hpp:152
twoD_diffusion_problem::klFunc
Teuchos::RCP< KL_Diffusion_Func > klFunc
Definition: twoD_diffusion_problem_tpetra.hpp:172
twoD_diffusion_problem::Tpetra_CrsGraph
Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > Tpetra_CrsGraph
Definition: twoD_diffusion_problem_tpetra.hpp:41
twoD_diffusion_problem::h
double h
Definition: twoD_diffusion_problem.hpp:113
LocalOrdinal
int LocalOrdinal
Definition: linear2d_diffusion_scalar_types.hpp:18
twoD_diffusion_problem::KL_Diffusion_Func
Definition: twoD_diffusion_problem_tpetra.hpp:150
twoD_diffusion_problem::x_init
Teuchos::RCP< Epetra_Vector > x_init
Initial guess.
Definition: twoD_diffusion_problem.hpp:134
GlobalOrdinal
int GlobalOrdinal
Definition: linear2d_diffusion_scalar_types.hpp:19
twoD_diffusion_problem::Tpetra_Vector
Tpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > Tpetra_Vector
Definition: twoD_diffusion_problem_tpetra.hpp:38
val
expr val()
f
ScalarType f(const Teuchos::Array< ScalarType > &x, double a, double b)
Definition: gram_schmidt_example3.cpp:64
twoD_diffusion_problem::point
Teuchos::Array< double > point
Array to store a point for basis evaluation.
Definition: twoD_diffusion_problem.hpp:164
twoD_diffusion_problem::b
Teuchos::RCP< Epetra_Vector > b
Deterministic RHS.
Definition: twoD_diffusion_problem.hpp:155
Node
Tpetra::KokkosClassic::DefaultNode::DefaultNodeType Node
Definition: experimental/linear2d_diffusion_pce.cpp:47
twoD_diffusion_problem::Tpetra_CrsMatrix
Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > Tpetra_CrsMatrix
Definition: twoD_diffusion_problem_tpetra.hpp:40
twoD_diffusion_problem::p_names
Teuchos::RCP< Teuchos::Array< std::string > > p_names
Parameter names.
Definition: twoD_diffusion_problem.hpp:146
twoD_diffusion_problem::graph
Teuchos::RCP< Epetra_CrsGraph > graph
Jacobian graph.
Definition: twoD_diffusion_problem.hpp:149
g
ScalarType g(const Teuchos::Array< ScalarType > &x, const ScalarType &y)
Definition: gram_schmidt_example3.cpp:75
twoD_diffusion_problem::eliminate_bcs
bool eliminate_bcs
Definition: twoD_diffusion_problem.hpp:125
twoD_diffusion_problem
A linear 2-D diffusion problem.
Definition: twoD_diffusion_problem.hpp:26
Generated on Fri Nov 22 2024 09:18:16 for Stokhos Package Browser (Single Doxygen Collection) by   doxygen 1.8.5