Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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
42template <typename Scalar, typename MeshScalar, typename BasisScalar,
43 typename LocalOrdinal, typename GlobalOrdinal, typename Node>
46twoD_diffusion_problem(
47 const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
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->getLocalNumElements();
134 ArrayView<const GlobalOrdinal> MyGlobalElements =
135 x_map->getLocalElementList ();
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));
183 lnFunc = rcp(new LogNormal_Diffusion_Func<KL_Diffusion_Func>(*klFunc));
184}
185
186// Overridden from EpetraExt::ModelEvaluator
187template <typename Scalar, typename MeshScalar, typename BasisScalar,
188 typename LocalOrdinal, typename GlobalOrdinal, typename Node>
189Teuchos::RCP<const
190 typename twoD_diffusion_problem<Scalar,MeshScalar,BasisScalar,
192 >::Tpetra_Map>
194 Node>::
195get_x_map() const
196{
197 return x_map;
198}
199
200template <typename Scalar, typename MeshScalar, typename BasisScalar,
201 typename LocalOrdinal, typename GlobalOrdinal, typename Node>
202Teuchos::RCP<const
203 typename twoD_diffusion_problem<Scalar,MeshScalar,BasisScalar,
205 >::Tpetra_Map>
207 Node>::
208get_f_map() const
209{
210 return x_map;
211}
212
213template <typename Scalar, typename MeshScalar, typename BasisScalar,
214 typename LocalOrdinal, typename GlobalOrdinal, typename Node>
215Teuchos::RCP<const
216 typename twoD_diffusion_problem<Scalar,MeshScalar,BasisScalar,
218 >::Tpetra_Map>
220 Node>::
221get_p_map(LocalOrdinal l) const
222{
223 TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
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
232template <typename Scalar, typename MeshScalar, typename BasisScalar,
233 typename LocalOrdinal, typename GlobalOrdinal, typename Node>
234Teuchos::RCP<const
235 typename twoD_diffusion_problem<Scalar,MeshScalar,BasisScalar,
237 >::Tpetra_Map>
239 Node>::
240get_g_map(LocalOrdinal l) const
241{
242 TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
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
251template <typename Scalar, typename MeshScalar, typename BasisScalar,
252 typename LocalOrdinal, typename GlobalOrdinal, typename Node>
253Teuchos::RCP<const Teuchos::Array<std::string> >
255 Node>::
256get_p_names(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_names(): " <<
262 "Invalid parameter index l = " << l << std::endl);
263
264 return p_names;
265}
266
267template <typename Scalar, typename MeshScalar, typename BasisScalar,
268 typename LocalOrdinal, typename GlobalOrdinal, typename Node>
269Teuchos::RCP<const
270 typename twoD_diffusion_problem<Scalar,MeshScalar,BasisScalar,
272 >::Tpetra_Vector>
274 Node>::
275get_x_init() const
276{
277 return x_init;
278}
279
280template <typename Scalar, typename MeshScalar, typename BasisScalar,
281 typename LocalOrdinal, typename GlobalOrdinal, typename Node>
282Teuchos::RCP<const
283 typename twoD_diffusion_problem<Scalar,MeshScalar,BasisScalar,
285 >::Tpetra_Vector>
287 Node>::
288get_p_init(LocalOrdinal l) const
289{
290 TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
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
299template <typename Scalar, typename MeshScalar, typename BasisScalar,
300 typename LocalOrdinal, typename GlobalOrdinal, typename Node>
301Teuchos::RCP<typename twoD_diffusion_problem<Scalar,MeshScalar,BasisScalar,
303 >::Tpetra_CrsMatrix>
305 Node>::
306create_W() const
307{
308 Teuchos::RCP<Tpetra_CrsMatrix> AA =
309 Teuchos::rcp(new Tpetra_CrsMatrix(graph));
310 AA->fillComplete();
311 return AA;
312}
313
314template <typename Scalar, typename MeshScalar, typename BasisScalar,
315 typename LocalOrdinal, typename GlobalOrdinal, typename Node>
316void
319computeResidual(const Tpetra_Vector& x,
320 const Tpetra_Vector& p,
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
332template <typename Scalar, typename MeshScalar, typename BasisScalar,
333 typename LocalOrdinal, typename GlobalOrdinal, typename Node>
334void
337computeJacobian(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
347template <typename Scalar, typename MeshScalar, typename BasisScalar,
348 typename LocalOrdinal, typename GlobalOrdinal, typename Node>
349void
352computeResponse(const Tpetra_Vector& x,
353 const Tpetra_Vector& p,
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
362template <typename Scalar, typename MeshScalar, typename BasisScalar,
363 typename LocalOrdinal, typename GlobalOrdinal, typename Node>
364template <typename FuncT>
365void
368computeA(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->getLocalNumElements();
379 ArrayView<const GlobalOrdinal> MyGlobalElements =
380 x_map->getLocalElementList ();
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
436template <typename Scalar, typename MeshScalar, typename BasisScalar,
437 typename LocalOrdinal, typename GlobalOrdinal, typename Node>
439 Node>::KL_Diffusion_Func::
440KL_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 =
461 Teuchos::rcp(new Stokhos::KL::ExponentialRandomField<MeshScalar>(rfParams));
462}
463
464template <typename Scalar, typename MeshScalar, typename BasisScalar,
465 typename LocalOrdinal, typename GlobalOrdinal, typename Node>
466Scalar
468 Node>::KL_Diffusion_Func::
469operator() (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}
expr val()
Class representing a KL expansion of an exponential random field.
A linear 2-D diffusion problem.
Teuchos::RCP< Epetra_Vector > p_init
Initial parameters.
Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > Tpetra_CrsMatrix
Teuchos::RCP< Epetra_Map > x_map
Solution vector map.
Teuchos::RCP< Epetra_Vector > x_init
Initial guess.
Teuchos::RCP< Epetra_CrsMatrix > A
Matrix to store deterministic operator.
Teuchos::RCP< Teuchos::Array< std::string > > p_names
Parameter names.
Teuchos::RCP< Epetra_Map > p_map
Parameter vector map.
Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > Tpetra_CrsGraph
Teuchos::RCP< LogNormal_Diffusion_Func< KL_Diffusion_Func > > lnFunc
Teuchos::Array< MeshPoint > mesh
Teuchos::Array< double > point
Array to store a point for basis evaluation.
Tpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > Tpetra_Vector
Teuchos::RCP< KL_Diffusion_Func > klFunc
Teuchos::RCP< Epetra_Vector > b
Deterministic RHS.
Teuchos::RCP< Epetra_CrsGraph > graph
Jacobian graph.
Teuchos::RCP< Epetra_Map > g_map
Response vector map.
KokkosClassic::DefaultNode::DefaultNodeType Node
ScalarType g(const Teuchos::Array< ScalarType > &x, const ScalarType &y)
ScalarType f(const Teuchos::Array< ScalarType > &x, double a, double b)
Teuchos::RCP< Stokhos::KL::ExponentialRandomField< MeshScalar > > rf