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 //
4 // Stokhos Package
5 // Copyright (2009) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 template <typename Scalar, typename MeshScalar, typename BasisScalar,
43  typename LocalOrdinal, typename GlobalOrdinal, typename Node>
45  Node>::
47  const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
48  LocalOrdinal n, LocalOrdinal d,
49  BasisScalar s, BasisScalar mu,
50  bool log_normal_,
51  bool eliminate_bcs_) :
52  mesh(n*n),
53  log_normal(log_normal_),
54  eliminate_bcs(eliminate_bcs_)
55 {
56  using Teuchos::Array;
57  using Teuchos::ArrayView;
58  using Teuchos::arrayView;
59  using Teuchos::ArrayRCP;
60  using Teuchos::RCP;
61  using Teuchos::rcp;
62  using Tpetra::global_size_t;
63 
65  // Construct the mesh.
66  // The mesh is uniform and the nodes are numbered
67  // LEFT to RIGHT, DOWN to UP.
68  //
69  // 5-6-7-8-9
70  // | | | | |
71  // 0-1-2-3-4
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;
82  mesh[idx].x = x;
83  mesh[idx].y = y;
84  if (i == 0 || i == n-1 || j == 0 || j == n-1)
85  mesh[idx].boundary = true;
86  if (i != 0)
87  mesh[idx].left = idx-1;
88  if (i != n-1)
89  mesh[idx].right = idx+1;
90  if (j != 0)
91  mesh[idx].down = idx-n;
92  if (j != n-1)
93  mesh[idx].up = idx+n;
94  if (!(eliminate_bcs && mesh[idx].boundary))
95  global_dof_indices.push_back(idx);
96  }
97  }
98 
99  // Solution vector map
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);
109 
110  // Initial guess, initialized to 0.0
111  x_init = Tpetra::createVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(x_map);
112  x_init->putScalar(0.0);
113 
114  // Parameter vector map
115  p_map = Tpetra::createLocalMap<LocalOrdinal,GlobalOrdinal>(d, comm);
116 
117  // Response vector map
118  g_map = Tpetra::createLocalMap<LocalOrdinal,GlobalOrdinal>(1, comm);
119 
120  // Initial parameters
121  p_init = Tpetra::createVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(p_map);
122  p_init->putScalar(0.0);
123 
124  // Parameter names
125  p_names = Teuchos::rcp(new Array<std::string>(d));
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();
130  }
131 
132  // Build Jacobian graph
133  size_t NumMyElements = x_map->getNodeNumElements();
134  ArrayView<const GlobalOrdinal> MyGlobalElements =
135  x_map->getNodeElementList ();
136  graph = rcp(new Tpetra_CrsGraph(x_map, 5));
137  for (size_t i=0; i<NumMyElements; ++i ) {
138 
139  // Center
140  GlobalOrdinal global_idx = MyGlobalElements[i];
141  graph->insertGlobalIndices(global_idx, arrayView(&global_idx, 1));
142 
143  if (!mesh[global_idx].boundary) {
144  // Down
145  if (!(eliminate_bcs && mesh[mesh[global_idx].down].boundary))
146  graph->insertGlobalIndices(global_idx,
147  arrayView(&mesh[global_idx].down,1));
148 
149  // Left
150  if (!(eliminate_bcs && mesh[mesh[global_idx].left].boundary))
151  graph->insertGlobalIndices(global_idx,
152  arrayView(&mesh[global_idx].left,1));
153 
154  // Right
155  if (!(eliminate_bcs && mesh[mesh[global_idx].right].boundary))
156  graph->insertGlobalIndices(global_idx,
157  arrayView(&mesh[global_idx].right,1));
158 
159  // Up
160  if (!(eliminate_bcs && mesh[mesh[global_idx].up].boundary))
161  graph->insertGlobalIndices(global_idx,
162  arrayView(&mesh[global_idx].up,1));
163  }
164  }
165  graph->fillComplete();
166 
167  // Construct deterministic operator
168  A = rcp(new Tpetra_CrsMatrix(graph));
169 
170  // Construct the RHS vector.
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)
176  b_view[i] = 0;
177  else
178  b_view[i] = 1;
179  }
180 
181  // Diffusion functions
182  klFunc = rcp(new KL_Diffusion_Func(xyLeft, xyRight, mu, s, 1.0, d));
184 }
185 
186 // Overridden from EpetraExt::ModelEvaluator
187 template <typename Scalar, typename MeshScalar, typename BasisScalar,
188  typename LocalOrdinal, typename GlobalOrdinal, typename Node>
189 Teuchos::RCP<const
190  typename twoD_diffusion_problem<Scalar,MeshScalar,BasisScalar,
192  >::Tpetra_Map>
194  Node>::
195 get_x_map() const
196 {
197  return x_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,
205  >::Tpetra_Map>
207  Node>::
208 get_f_map() const
209 {
210  return x_map;
211 }
212 
213 template <typename Scalar, typename MeshScalar, typename BasisScalar,
214  typename LocalOrdinal, typename GlobalOrdinal, typename Node>
215 Teuchos::RCP<const
216  typename twoD_diffusion_problem<Scalar,MeshScalar,BasisScalar,
218  >::Tpetra_Map>
220  Node>::
221 get_p_map(LocalOrdinal l) const
222 {
224  std::logic_error,
225  std::endl <<
226  "Error! twoD_diffusion_problem::get_p_map(): " <<
227  "Invalid parameter index l = " << l << std::endl);
228 
229  return p_map;
230 }
231 
232 template <typename Scalar, typename MeshScalar, typename BasisScalar,
233  typename LocalOrdinal, typename GlobalOrdinal, typename Node>
234 Teuchos::RCP<const
235  typename twoD_diffusion_problem<Scalar,MeshScalar,BasisScalar,
237  >::Tpetra_Map>
239  Node>::
240 get_g_map(LocalOrdinal l) const
241 {
243  std::logic_error,
244  std::endl <<
245  "Error! twoD_diffusion_problem::get_g_map(): " <<
246  "Invalid parameter index l = " << l << std::endl);
247 
248  return g_map;
249 }
250 
251 template <typename Scalar, typename MeshScalar, typename BasisScalar,
252  typename LocalOrdinal, typename GlobalOrdinal, typename Node>
255  Node>::
256 get_p_names(LocalOrdinal l) const
257 {
259  std::logic_error,
260  std::endl <<
261  "Error! twoD_diffusion_problem::get_p_names(): " <<
262  "Invalid parameter index l = " << l << std::endl);
263 
264  return p_names;
265 }
266 
267 template <typename Scalar, typename MeshScalar, typename BasisScalar,
268  typename LocalOrdinal, typename GlobalOrdinal, typename Node>
269 Teuchos::RCP<const
270  typename twoD_diffusion_problem<Scalar,MeshScalar,BasisScalar,
272  >::Tpetra_Vector>
274  Node>::
275 get_x_init() const
276 {
277  return x_init;
278 }
279 
280 template <typename Scalar, typename MeshScalar, typename BasisScalar,
281  typename LocalOrdinal, typename GlobalOrdinal, typename Node>
282 Teuchos::RCP<const
283  typename twoD_diffusion_problem<Scalar,MeshScalar,BasisScalar,
285  >::Tpetra_Vector>
287  Node>::
288 get_p_init(LocalOrdinal l) const
289 {
291  std::logic_error,
292  std::endl <<
293  "Error! twoD_diffusion_problem::get_p_init(): " <<
294  "Invalid parameter index l = " << l << std::endl);
295 
296  return p_init;
297 }
298 
299 template <typename Scalar, typename MeshScalar, typename BasisScalar,
300  typename LocalOrdinal, typename GlobalOrdinal, typename Node>
301 Teuchos::RCP<typename twoD_diffusion_problem<Scalar,MeshScalar,BasisScalar,
303  >::Tpetra_CrsMatrix>
305  Node>::
306 create_W() const
307 {
309  Teuchos::rcp(new Tpetra_CrsMatrix(graph));
310  AA->fillComplete();
311  return AA;
312 }
313 
314 template <typename Scalar, typename MeshScalar, typename BasisScalar,
315  typename LocalOrdinal, typename GlobalOrdinal, typename Node>
316 void
318  Node>::
319 computeResidual(const Tpetra_Vector& x,
320  const Tpetra_Vector& p,
321  Tpetra_Vector& f)
322 {
323  // f = A*x - b
324  if (log_normal)
325  computeA(*lnFunc, p, *A);
326  else
327  computeA(*klFunc, p, *A);
328  A->apply(x,f);
329  f.update(-1.0, *b, 1.0);
330 }
331 
332 template <typename Scalar, typename MeshScalar, typename BasisScalar,
333  typename LocalOrdinal, typename GlobalOrdinal, typename Node>
334 void
336  Node>::
337 computeJacobian(const Tpetra_Vector& x,
338  const Tpetra_Vector& p,
339  Tpetra_CrsMatrix& jac)
340 {
341  if (log_normal)
342  computeA(*lnFunc, p, jac);
343  else
344  computeA(*klFunc, p, jac);
345 }
346 
347 template <typename Scalar, typename MeshScalar, typename BasisScalar,
348  typename LocalOrdinal, typename GlobalOrdinal, typename Node>
349 void
351  Node>::
352 computeResponse(const Tpetra_Vector& x,
353  const Tpetra_Vector& p,
354  Tpetra_Vector& g)
355 {
356  // g = average of x
357  Teuchos::ArrayRCP<Scalar> g_view = g.get1dViewNonConst();
358  x.meanValue(g_view());
359  g_view[0] *= Scalar(x.getGlobalLength()) / Scalar(mesh.size());
360 }
361 
362 template <typename Scalar, typename MeshScalar, typename BasisScalar,
363  typename LocalOrdinal, typename GlobalOrdinal, typename Node>
364 template <typename FuncT>
365 void
367  Node>::
368 computeA(const FuncT& func, const Tpetra_Vector& p, Tpetra_CrsMatrix& jac)
369 {
370  using Teuchos::ArrayView;
371  using Teuchos::arrayView;
372 
373  jac.resumeFill();
374  jac.setAllToScalar(0.0);
375 
376  Teuchos::ArrayRCP<const Scalar> p_view = p.get1dView();
377  Teuchos::Array<Scalar> rv(p_view());
378  size_t NumMyElements = x_map->getNodeNumElements();
379  ArrayView<const GlobalOrdinal> MyGlobalElements =
380  x_map->getNodeElementList ();
381  MeshScalar h2 = h*h;
382  Scalar val;
383 
384  for(size_t i=0 ; i<NumMyElements; ++i ) {
385 
386  // Center
387  GlobalOrdinal global_idx = MyGlobalElements[i];
388  if (mesh[global_idx].boundary) {
389  val = 1.0;
390  jac.replaceGlobalValues(global_idx, arrayView(&global_idx,1),
391  arrayView(&val,1));
392  }
393  else {
394  Scalar a_down =
395  -func(mesh[global_idx].x, mesh[global_idx].y-h/2.0, rv)/h2;
396  Scalar a_left =
397  -func(mesh[global_idx].x-h/2.0, mesh[global_idx].y, rv)/h2;
398  Scalar a_right =
399  -func(mesh[global_idx].x+h/2.0, mesh[global_idx].y, rv)/h2;
400  Scalar a_up =
401  -func(mesh[global_idx].x, mesh[global_idx].y+h/2.0, rv)/h2;
402 
403  // Center
404  val = -(a_down + a_left + a_right + a_up);
405  jac.replaceGlobalValues(global_idx, arrayView(&global_idx,1),
406  arrayView(&val,1));
407 
408  // Down
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));
413 
414  // Left
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));
419 
420  // Right
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));
425 
426  // Up
427  if (!(eliminate_bcs && mesh[mesh[global_idx].up].boundary))
428  jac.replaceGlobalValues(global_idx,
429  arrayView(&mesh[global_idx].up,1),
430  arrayView(&a_up,1));
431  }
432  }
433  jac.fillComplete();
434 }
435 
436 template <typename Scalar, typename MeshScalar, typename BasisScalar,
437  typename LocalOrdinal, typename GlobalOrdinal, typename Node>
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)
443 {
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);
448  LocalOrdinal ndim = 2;
449  Teuchos::Array<MeshScalar> domain_upper(ndim), domain_lower(ndim),
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;
455  }
456  rfParams.set("Domain Upper Bounds", domain_upper);
457  rfParams.set("Domain Lower Bounds", domain_lower);
458  rfParams.set("Correlation Lengths", correlation_length);
459 
460  rf =
462 }
463 
464 template <typename Scalar, typename MeshScalar, typename BasisScalar,
465  typename LocalOrdinal, typename GlobalOrdinal, typename Node>
466 Scalar
468  Node>::KL_Diffusion_Func::
469 operator() (MeshScalar x, MeshScalar y, const Teuchos::Array<Scalar>& rv) const
470 {
471  point[0] = x;
472  point[1] = y;
473  return rf->evaluate(point, rv);
474 }
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
expr val()
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.