Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
uq_handbook/SimpleME.cpp
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#include "SimpleME.hpp"
43#include "Teuchos_Assert.hpp"
44#include "Stokhos_Epetra.hpp"
45#include "Stokhos_Sacado.hpp"
46#include "Sacado.hpp"
47
48namespace {
49
50 template <typename ScalarA, typename ScalarX>
51 void func(const ScalarA& a, const ScalarX x[2], ScalarX y[2]) {
52 y[0] = a*a - x[0];
53 y[1] = x[1]*x[1] - x[0];
54 }
55
56}
57
58SimpleME::SimpleME(const Teuchos::RCP<const Epetra_Comm>& comm)
59{
60 // Solution vector map
61 x_map = Teuchos::rcp(new Epetra_Map(2, 0, *comm));
62
63 // Overlapped solution vector map
64 x_overlapped_map = Teuchos::rcp(new Epetra_LocalMap(2, 0, *comm));
65
66 // Importer
67 importer = Teuchos::rcp(new Epetra_Import(*x_overlapped_map, *x_map));
68
69 // Initial guess, initialized to 1.5
70 x_init = Teuchos::rcp(new Epetra_Vector(*x_map));
71 x_init->PutScalar(1.5);
72
73 // Overlapped solution vector
74 x_overlapped = Teuchos::rcp(new Epetra_Vector(*x_overlapped_map));
75
76 // Parameter vector map
77 p_map = Teuchos::rcp(new Epetra_LocalMap(1, 0, *comm));
78
79 // Initial parameters
80 p_init = Teuchos::rcp(new Epetra_Vector(*p_map));
81 (*p_init)[0] = 2.0;
82
83 // Parameter names
84 p_names = Teuchos::rcp(new Teuchos::Array<std::string>(1));
85 (*p_names)[0] = "alpha";
86
87 // Jacobian graph (dense 2x2 matrix)
88 graph = Teuchos::rcp(new Epetra_CrsGraph(Copy, *x_map, 2));
89 int indices[2];
90 indices[0] = 0; indices[1] = 1;
91 if (x_map->MyGID(0))
92 graph->InsertGlobalIndices(0, 2, indices);
93 if (x_map->MyGID(1))
94 graph->InsertGlobalIndices(1, 2, indices);
95 graph->FillComplete();
96 graph->OptimizeStorage();
97}
98
99// Overridden from EpetraExt::ModelEvaluator
100
101Teuchos::RCP<const Epetra_Map>
103{
104 return x_map;
105}
106
107Teuchos::RCP<const Epetra_Map>
109{
110 return x_map;
111}
112
113Teuchos::RCP<const Epetra_Map>
114SimpleME::get_p_map(int l) const
115{
116 TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
117 std::logic_error,
118 std::endl <<
119 "Error! SimpleME::get_p_map(): " <<
120 "Invalid parameter index l = " << l << std::endl);
121
122 return p_map;
123}
124
125Teuchos::RCP<const Teuchos::Array<std::string> >
126SimpleME::get_p_names(int l) const
127{
128 TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
129 std::logic_error,
130 std::endl <<
131 "Error! SimpleME::get_p_names(): " <<
132 "Invalid parameter index l = " << l << std::endl);
133
134 return p_names;
135}
136
137Teuchos::RCP<const Epetra_Vector>
139{
140 return x_init;
141}
142
143Teuchos::RCP<const Epetra_Vector>
144SimpleME::get_p_init(int l) const
145{
146 TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
147 std::logic_error,
148 std::endl <<
149 "Error! SimpleME::get_p_init(): " <<
150 "Invalid parameter index l = " << l << std::endl);
151
152 return p_init;
153}
154
155Teuchos::RCP<Epetra_Operator>
156SimpleME::create_W() const
157{
158 Teuchos::RCP<Epetra_CrsMatrix> A =
159 Teuchos::rcp(new Epetra_CrsMatrix(Copy, *graph));
160 A->FillComplete();
161 A->OptimizeStorage();
162 return A;
163}
164
165EpetraExt::ModelEvaluator::InArgs
167{
168 InArgsSetup inArgs;
169 inArgs.setModelEvalDescription("Simple Model Evaluator");
170
171 // Deterministic InArgs
172 inArgs.setSupports(IN_ARG_x,true);
173 inArgs.set_Np(1); // 1 parameter vector
174
175 // Stochastic InArgs
176 inArgs.setSupports(IN_ARG_x_sg,true);
177 inArgs.setSupports(IN_ARG_p_sg, 0, true); // 1 SG parameter vector
178 inArgs.setSupports(IN_ARG_sg_basis,true);
179 inArgs.setSupports(IN_ARG_sg_quadrature,true);
180 inArgs.setSupports(IN_ARG_sg_expansion,true);
181
182 return inArgs;
183}
184
185EpetraExt::ModelEvaluator::OutArgs
187{
188 OutArgsSetup outArgs;
189 outArgs.setModelEvalDescription("Simple Model Evaluator");
190
191 // Deterministic OutArgs
192 outArgs.set_Np_Ng(1, 0);
193 outArgs.setSupports(OUT_ARG_f,true);
194 outArgs.setSupports(OUT_ARG_W,true);
195
196 // Stochastic OutArgs
197 outArgs.setSupports(OUT_ARG_f_sg,true);
198 outArgs.setSupports(OUT_ARG_W_sg,true);
199
200 return outArgs;
201}
202
203void
204SimpleME::evalModel(const InArgs& inArgs, const OutArgs& outArgs) const
205{
206 //
207 // Determinisic calculation
208 //
209
210 // Solution vector
211 Teuchos::RCP<const Epetra_Vector> x = inArgs.get_x();
212 if (x != Teuchos::null) {
213 x_overlapped->Import(*x, *importer, Insert);
214 double x0 = (*x_overlapped)[0];
215 double x1 = (*x_overlapped)[1];
216
217 // Parameters
218 Teuchos::RCP<const Epetra_Vector> p = inArgs.get_p(0);
219 if (p == Teuchos::null)
220 p = p_init;
221 double a = (*p)[0];
222
223 // Residual
224 // f = | a*a - x0 |
225 // | x1*x1 - x0 |
226 // where a = p[0].
227 Teuchos::RCP<Epetra_Vector> f = outArgs.get_f();
228 if (f != Teuchos::null) {
229 double x[2] = { x0, x1 };
230 double y[2];
231 func(a, x, y);
232
233 if (x_map->MyGID(0)) {
234 int row = 0;
235 f->ReplaceGlobalValues(1, &y[0], &row);
236 }
237 if (x_map->MyGID(1)) {
238 int row = 1;
239 f->ReplaceGlobalValues(1, &y[1], &row);
240 }
241 }
242
243 // Jacobian
244 // J = | -1 0 |
245 // | -1 2*x1 |
246 Teuchos::RCP<Epetra_Operator> W = outArgs.get_W();
247 if (W != Teuchos::null) {
248 typedef Sacado::Fad::SFad<double,2> fad_type;
249 fad_type x[2], y[2];
250 x[0] = fad_type(2, 0, x0);
251 x[1] = fad_type(2, 1, x1);
252 func(a, x, y);
253
254 Teuchos::RCP<Epetra_CrsMatrix> jac =
255 Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(W, true);
256 int indices[2] = { 0, 1 };
257 if (x_map->MyGID(0)) {
258 int row = 0;
259 jac->ReplaceGlobalValues(row, 2, y[0].dx(), indices);
260 }
261 if (x_map->MyGID(1)) {
262 int row = 1;
263 jac->ReplaceGlobalValues(row, 2, y[1].dx(), indices);
264 }
265 }
266 }
267
268 //
269 // Stochastic Galerkin calculation
270 //
271
272 // Stochastic solution vector
273 InArgs::sg_const_vector_t x_sg = inArgs.get_x_sg();
274 if (x_sg != Teuchos::null) {
275
276 // Get stochastic expansion data
277 Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> > basis =
278 inArgs.get_sg_basis();
279 Teuchos::RCP<Stokhos::OrthogPolyExpansion<int,double> > expn =
280 inArgs.get_sg_expansion();
283
284 pce_type x0(expn), x1(expn);
285 for (int i=0; i<basis->size(); i++) {
286 x_overlapped->Import((*x_sg)[i], *importer, Insert);
287 x0.fastAccessCoeff(i) = (*x_overlapped)[0];
288 x1.fastAccessCoeff(i) = (*x_overlapped)[1];
289 }
290
291 // Stochastic parameters
292 InArgs::sg_const_vector_t p_sg = inArgs.get_p_sg(0);
293 pce_type a(expn);
294 if (p_sg != Teuchos::null) {
295 for (int i=0; i<basis->size(); i++) {
296 a.fastAccessCoeff(i) = (*p_sg)[i][0];
297 }
298 }
299
300 // Stochastic residual
301 // f[i] = | <a*a - x0, psi_i>/<psi_i^2> |
302 // | <x1*x1 - x0, psi_i>/<psi_i^2> |
303 OutArgs::sg_vector_t f_sg = outArgs.get_f_sg();
304 if (f_sg != Teuchos::null) {
305 pce_type x[2] = { x0, x1 };
306 pce_type y[2];
307 func(a, x, y);
308
309 if (x_map->MyGID(0)) {
310 int row = 0;
311 for (int i=0; i<basis->size(); i++) {
312 double c = y[0].coeff(i);
313 (*f_sg)[i].ReplaceGlobalValues(1, &c, &row);
314 }
315 }
316 if (x_map->MyGID(1)) {
317 int row = 1;
318 for (int i=0; i<basis->size(); i++) {
319 double c = y[1].coeff(i);
320 (*f_sg)[i].ReplaceGlobalValues(1, &c, &row);
321 }
322 }
323 }
324
325 // Stochastic Jacobian
326 // J[0] = | -1 0 |, J[i] = | 0 0 |, i > 0
327 // | -1 2*x0[0] | | 0 2*x0[i] |
328 OutArgs::sg_operator_t W_sg = outArgs.get_W_sg();
329 if (W_sg != Teuchos::null) {
330 typedef Sacado::Fad::SFad<pce_type,2> fad_type;
331 fad_type x[2], y[2];
332 x[0] = fad_type(2, 0, x0);
333 x[1] = fad_type(2, 1, x1);
334 func(a, x, y);
335
336 for (int i=0; i<basis->size(); i++) {
337 Teuchos::RCP<Epetra_CrsMatrix> jac =
338 Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(W_sg->getCoeffPtr(i),
339 true);
340 int indices[2] = { 0, 1 };
341 if (x_map->MyGID(0)) {
342 int row = 0;
343 double values[2] = { y[0].dx(0).coeff(i), y[0].dx(1).coeff(i) };
344 jac->ReplaceGlobalValues(row, 2, values, indices);
345 }
346 if (x_map->MyGID(1)) {
347 int row = 1;
348 double values[2] = { y[1].dx(0).coeff(i), y[1].dx(1).coeff(i) };
349 jac->ReplaceGlobalValues(row, 2, values, indices);
350 }
351 }
352 }
353 }
354}
expr expr expr dx(i, j)
Stokhos::StandardStorage< int, double > storage_type
Sacado::Fad::DFad< double > fad_type
Teuchos::RCP< Epetra_Import > importer
Importer to overlapped distribution.
Teuchos::RCP< Teuchos::Array< std::string > > p_names
Parameter names.
OutArgs createOutArgs() const
Create OutArgs.
InArgs createInArgs() const
Create InArgs.
Teuchos::RCP< const Epetra_Vector > get_p_init(int l) const
Return initial parameters.
Teuchos::RCP< Epetra_Operator > create_W() const
Create W = alpha*M + beta*J matrix.
Teuchos::RCP< Epetra_CrsGraph > graph
Jacobian graph.
void evalModel(const InArgs &inArgs, const OutArgs &outArgs) const
Evaluate model on InArgs.
Teuchos::RCP< Epetra_Map > x_map
Solution vector map.
Teuchos::RCP< Epetra_Vector > p_init
Initial parameters.
Teuchos::RCP< Epetra_Vector > x_init
Initial guess.
Teuchos::RCP< Epetra_Map > x_overlapped_map
Overlapped solution vector map.
Teuchos::RCP< Epetra_Vector > x_overlapped
Overlapped solution vector.
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int l) const
Return array of parameter names.
Teuchos::RCP< const Epetra_Vector > get_x_init() const
Return initial solution.
Teuchos::RCP< const Epetra_Map > get_f_map() const
Return residual vector map.
Teuchos::RCP< Epetra_Map > p_map
Parameter vector map.
Teuchos::RCP< const Epetra_Map > get_p_map(int l) const
Return parameter vector map.
Teuchos::RCP< const Epetra_Map > get_x_map() const
Return solution vector map.
SimpleME(const Teuchos::RCP< const Epetra_Comm > &comm)
Constructor.
Sacado::ETPCE::OrthogPoly< double, Stokhos::StandardStorage< int, double > > pce_type
ScalarType f(const Teuchos::Array< ScalarType > &x, double a, double b)
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
Definition csr_vector.h:265
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType ValueType * y
Definition csr_vector.h:267