Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
epetra/SimpleME.cpp
Go to the documentation of this file.
1// $Id$
2// $Source$
3// @HEADER
4// ***********************************************************************
5//
6// Stokhos Package
7// Copyright (2009) Sandia Corporation
8//
9// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
10// license for use of this work by or on behalf of the U.S. Government.
11//
12// Redistribution and use in source and binary forms, with or without
13// modification, are permitted provided that the following conditions are
14// met:
15//
16// 1. Redistributions of source code must retain the above copyright
17// notice, this list of conditions and the following disclaimer.
18//
19// 2. Redistributions in binary form must reproduce the above copyright
20// notice, this list of conditions and the following disclaimer in the
21// documentation and/or other materials provided with the distribution.
22//
23// 3. Neither the name of the Corporation nor the names of the
24// contributors may be used to endorse or promote products derived from
25// this software without specific prior written permission.
26//
27// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38//
39// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
40//
41// ***********************************************************************
42// @HEADER
43
44#include "SimpleME.hpp"
45#include "Teuchos_Assert.hpp"
46#include "Stokhos_Epetra.hpp"
47
48SimpleME::SimpleME(const Teuchos::RCP<const Epetra_Comm>& comm)
49{
50 // Solution vector map
51 x_map = Teuchos::rcp(new Epetra_Map(2, 0, *comm));
52
53 // Overlapped solution vector map
54 x_overlapped_map = Teuchos::rcp(new Epetra_LocalMap(2, 0, *comm));
55
56 // Importer
57 importer = Teuchos::rcp(new Epetra_Import(*x_overlapped_map, *x_map));
58
59 // Initial guess, initialized to 1.5
60 x_init = Teuchos::rcp(new Epetra_Vector(*x_map));
61 x_init->PutScalar(1.5);
62
63 // Overlapped solution vector
64 x_overlapped = Teuchos::rcp(new Epetra_Vector(*x_overlapped_map));
65
66 // Parameter vector map
67 p_map = Teuchos::rcp(new Epetra_LocalMap(1, 0, *comm));
68
69 // Initial parameters
70 p_init = Teuchos::rcp(new Epetra_Vector(*p_map));
71 (*p_init)[0] = 2.0;
72
73 // Parameter names
74 p_names = Teuchos::rcp(new Teuchos::Array<std::string>(1));
75 (*p_names)[0] = "alpha";
76
77 // Jacobian graph (dense 2x2 matrix)
78 graph = Teuchos::rcp(new Epetra_CrsGraph(Copy, *x_map, 2));
79 int indices[2];
80 indices[0] = 0; indices[1] = 1;
81 if (x_map->MyGID(0))
82 graph->InsertGlobalIndices(0, 2, indices);
83 if (x_map->MyGID(1))
84 graph->InsertGlobalIndices(1, 2, indices);
85 graph->FillComplete();
86 graph->OptimizeStorage();
87}
88
89// Overridden from EpetraExt::ModelEvaluator
90
91Teuchos::RCP<const Epetra_Map>
93{
94 return x_map;
95}
96
97Teuchos::RCP<const Epetra_Map>
99{
100 return x_map;
101}
102
103Teuchos::RCP<const Epetra_Map>
105{
106 TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
107 std::logic_error,
108 std::endl <<
109 "Error! SimpleME::get_p_map(): " <<
110 "Invalid parameter index l = " << l << std::endl);
111
112 return p_map;
113}
114
115Teuchos::RCP<const Teuchos::Array<std::string> >
117{
118 TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
119 std::logic_error,
120 std::endl <<
121 "Error! SimpleME::get_p_names(): " <<
122 "Invalid parameter index l = " << l << std::endl);
123
124 return p_names;
125}
126
127Teuchos::RCP<const Epetra_Vector>
129{
130 return x_init;
131}
132
133Teuchos::RCP<const Epetra_Vector>
135{
136 TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
137 std::logic_error,
138 std::endl <<
139 "Error! SimpleME::get_p_init(): " <<
140 "Invalid parameter index l = " << l << std::endl);
141
142 return p_init;
143}
144
145Teuchos::RCP<Epetra_Operator>
147{
148 Teuchos::RCP<Epetra_CrsMatrix> A =
149 Teuchos::rcp(new Epetra_CrsMatrix(Copy, *graph));
150 A->FillComplete();
151 A->OptimizeStorage();
152 return A;
153}
154
155EpetraExt::ModelEvaluator::InArgs
157{
158 InArgsSetup inArgs;
159 inArgs.setModelEvalDescription("Simple Model Evaluator");
160
161 // Deterministic InArgs
162 inArgs.setSupports(IN_ARG_x,true);
163 inArgs.set_Np(1); // 1 parameter vector
164
165 // Stochastic InArgs
166 inArgs.setSupports(IN_ARG_x_sg,true);
167 inArgs.setSupports(IN_ARG_p_sg, 0, true); // 1 SG parameter vector
168 inArgs.setSupports(IN_ARG_sg_basis,true);
169 inArgs.setSupports(IN_ARG_sg_quadrature,true);
170 inArgs.setSupports(IN_ARG_sg_expansion,true);
171
172 return inArgs;
173}
174
175EpetraExt::ModelEvaluator::OutArgs
177{
178 OutArgsSetup outArgs;
179 outArgs.setModelEvalDescription("Simple Model Evaluator");
180
181 // Deterministic OutArgs
182 outArgs.set_Np_Ng(1, 0);
183 outArgs.setSupports(OUT_ARG_f,true);
184 outArgs.setSupports(OUT_ARG_W,true);
185
186 // Stochastic OutArgs
187 outArgs.setSupports(OUT_ARG_f_sg,true);
188 outArgs.setSupports(OUT_ARG_W_sg,true);
189
190 return outArgs;
191}
192
193void
194SimpleME::evalModel(const InArgs& inArgs, const OutArgs& outArgs) const
195{
196 //
197 // Determinisic calculation
198 //
199
200 // Solution vector
201 Teuchos::RCP<const Epetra_Vector> x = inArgs.get_x();
202 if (x != Teuchos::null)
203 x_overlapped->Import(*x, *importer, Insert);
204 double x0 = (*x_overlapped)[0];
205 double x1 = (*x_overlapped)[1];
206
207 // Parameters
208 Teuchos::RCP<const Epetra_Vector> p = inArgs.get_p(0);
209 if (p == Teuchos::null)
210 p = p_init;
211 double a = (*p)[0];
212
213 // Residual
214 // f = | a*a - x0 |
215 // | x1*x1 - x0 |
216 // where a = p[0].
217 Teuchos::RCP<Epetra_Vector> f = outArgs.get_f();
218 if (f != Teuchos::null) {
219 int row;
220 double v;
221 if (x_map->MyGID(0)) {
222 row = 0;
223 v = a*a - x0;
224 f->ReplaceGlobalValues(1, &v, &row);
225 }
226 if (x_map->MyGID(1)) {
227 row = 1;
228 v = x1*x1 - x0;
229 f->ReplaceGlobalValues(1, &v, &row);
230 }
231 }
232
233 // Jacobian
234 // J = | -1 0 |
235 // | -1 2*x1 |
236 Teuchos::RCP<Epetra_Operator> W = outArgs.get_W();
237 if (W != Teuchos::null) {
238 Teuchos::RCP<Epetra_CrsMatrix> jac =
239 Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(W, true);
240 int indices[2] = { 0, 1 };
241 double values[2];
242 if (x_map->MyGID(0)) {
243 values[0] = -1.0; values[1] = 0.0;
244 jac->ReplaceGlobalValues(0, 2, values, indices);
245 }
246 if (x_map->MyGID(1)) {
247 values[0] = -1.0; values[1] = 2.0*x1;
248 jac->ReplaceGlobalValues(1, 2, values, indices);
249 }
250 }
251
252 //
253 // Stochastic Galerkin calculation
254 //
255
256 // Get stochastic expansion data
257 Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> > basis =
258 inArgs.get_sg_basis();
259 Teuchos::RCP<Stokhos::OrthogPolyExpansion<int,double> > expn =
260 inArgs.get_sg_expansion();
261
262 // Stochastic solution vector
263 InArgs::sg_const_vector_t x_sg = inArgs.get_x_sg();
264 Stokhos::OrthogPolyApprox<int,double> x0_sg(basis), x1_sg(basis);
265 if (x_sg != Teuchos::null) {
266 for (int i=0; i<basis->size(); i++) {
267 x_overlapped->Import((*x_sg)[i], *importer, Insert);
268 x0_sg[i] = (*x_overlapped)[0];
269 x1_sg[i] = (*x_overlapped)[1];
270 }
271 }
272
273 // Stochastic parameters
274 InArgs::sg_const_vector_t p_sg = inArgs.get_p_sg(0);
276 if (p_sg != Teuchos::null) {
277 for (int i=0; i<basis->size(); i++) {
278 a_sg[i] = (*p_sg)[i][0];
279 }
280 }
281
282 // Stochastic residual
283 // f[i] = | <a*a - x0, psi_i>/<psi_i^2> |
284 // | <x1*x1 - x0, psi_i>/<psi_i^2> |
285 OutArgs::sg_vector_t f_sg = outArgs.get_f_sg();
286 Stokhos::OrthogPolyApprox<int,double> tmp0_sg(basis), tmp1_sg(basis);
287 if (f_sg != Teuchos::null) {
288 int row;
289 if (x_map->MyGID(0)) {
290 row = 0;
291 expn->times(tmp0_sg, a_sg, a_sg);
292 expn->minus(tmp1_sg, tmp0_sg, x0_sg);
293 for (int i=0; i<basis->size(); i++)
294 (*f_sg)[i].ReplaceGlobalValues(1, &tmp1_sg[i], &row);
295 }
296 if (x_map->MyGID(1)) {
297 row = 1;
298 expn->times(tmp0_sg, x1_sg, x1_sg);
299 expn->minus(tmp1_sg, tmp0_sg, x0_sg);
300 for (int i=0; i<basis->size(); i++)
301 (*f_sg)[i].ReplaceGlobalValues(1, &tmp1_sg[i], &row);
302 }
303 }
304
305 // Stochastic Jacobian
306 // J[0] = | -1 0 |, J[i] = | 0 0 |, i > 0
307 // | -1 2*x0[0] | | 0 2*x0[i] |
308 OutArgs::sg_operator_t W_sg = outArgs.get_W_sg();
309 if (W_sg != Teuchos::null) {
310 for (int i=0; i<basis->size(); i++) {
311 Teuchos::RCP<Epetra_CrsMatrix> jac =
312 Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(W_sg->getCoeffPtr(i), true);
313 int indices[2] = { 0, 1 };
314 double values[2];
315 if (x_map->MyGID(0)) {
316 if (i == 0)
317 values[0] = -1.0;
318 else
319 values[0] = 0.0;
320 values[1] = 0.0;
321 jac->ReplaceGlobalValues(0, 2, values, indices);
322 }
323 if (x_map->MyGID(1)) {
324 if (i == 0)
325 values[0] = -1.0;
326 else
327 values[0] = 0.0;
328 values[1] = 2.0*x1_sg[i];
329 jac->ReplaceGlobalValues(1, 2, values, indices);
330 }
331 }
332 }
333}
Insert
Copy
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.
Class to store coefficients of a projection onto an orthogonal polynomial basis.
ScalarType f(const Teuchos::Array< ScalarType > &x, double a, double b)