Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_SGQuadModelEvaluator.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
48#include "Epetra_Map.h"
49#include "Epetra_Vector.h"
50#include "Teuchos_TimeMonitor.hpp"
51#include "Teuchos_Assert.hpp"
52
55 const Teuchos::RCP<EpetraExt::ModelEvaluator>& me_) :
56 me(me_),
57 num_p(0),
58 num_g(0),
59 x_dot_qp(),
60 x_qp(),
61 p_qp(),
62 f_qp(),
63 W_qp(),
64 dfdp_qp(),
65 g_qp(),
66 dgdx_qp(),
67 dgdx_dot_qp(),
68 dgdp_qp()
69{
70 // Create storage for x_dot, x, and p at a quad point
71 InArgs me_inargs = me->createInArgs();
72 num_p = me_inargs.Np();
73 if (me_inargs.supports(IN_ARG_x_dot))
74 x_dot_qp = Teuchos::rcp(new Epetra_Vector(*(me->get_x_map())));
75 if (me_inargs.supports(IN_ARG_x))
76 x_qp = Teuchos::rcp(new Epetra_Vector((*me->get_x_map())));
77 p_qp.resize(num_p);
78 for (int i=0; i<num_p; i++)
79 p_qp[i] = Teuchos::rcp(new Epetra_Vector(*(me->get_p_map(i))));
80
81 // Create storage for f and W at a quad point
82 OutArgs me_outargs = me->createOutArgs();
83 num_g = me_outargs.Ng();
84
85 // f
86 if (me_outargs.supports(OUT_ARG_f))
87 f_qp = Teuchos::rcp(new Epetra_Vector(*(me->get_f_map())));
88
89 // W
90 if (me_outargs.supports(OUT_ARG_W))
91 W_qp = me->create_W();
92
93 // df/dp
94 dfdp_qp.resize(num_p);
95 for (int i=0; i<num_p; i++)
96 if (me_outargs.supports(OUT_ARG_DfDp,i).supports(DERIV_MV_BY_COL))
97 dfdp_qp[i] = EpetraExt::ModelEvaluator::Derivative(
98 Teuchos::rcp(new Epetra_MultiVector(
99 *(me->get_f_map()),
100 me->get_p_map(i)->NumGlobalElements())));
101 else if (me_outargs.supports(OUT_ARG_DfDp,i).supports(DERIV_TRANS_MV_BY_ROW))
102 dfdp_qp[i] = EpetraExt::ModelEvaluator::Derivative(
103 Teuchos::rcp(new Epetra_MultiVector(
104 *(me->get_p_map(i)),
105 me->get_f_map()->NumGlobalElements())));
106 else if (me_outargs.supports(OUT_ARG_DfDp,i).supports(DERIV_LINEAR_OP))
107 dfdp_qp[i] = EpetraExt::ModelEvaluator::Derivative(
108 me->create_DfDp_op(i));
109
110 g_qp.resize(num_g);
111 dgdx_qp.resize(num_g);
112 dgdx_dot_qp.resize(num_g);
113 dgdp_qp.resize(num_g);
114 for (int i=0; i<num_g; i++) {
115
116 // g
117 g_qp[i] =
118 Teuchos::rcp(new Epetra_Vector(*(me->get_g_map(i))));
119
120 // dg/dx
121 if (me_outargs.supports(OUT_ARG_DgDx, i).supports(DERIV_TRANS_MV_BY_ROW))
122 dgdx_qp[i] = EpetraExt::ModelEvaluator::Derivative(
123 Teuchos::rcp(new Epetra_MultiVector(
124 *(me->get_x_map()),
125 me->get_g_map(i)->NumGlobalElements())));
126 else if (me_outargs.supports(OUT_ARG_DgDx, i).supports(DERIV_MV_BY_COL))
127 dgdx_qp[i] = EpetraExt::ModelEvaluator::Derivative(
128 Teuchos::rcp(new Epetra_MultiVector(
129 *(me->get_g_map(i)),
130 me->get_x_map()->NumGlobalElements())));
131 else if (me_outargs.supports(OUT_ARG_DgDx, i).supports(DERIV_LINEAR_OP))
132 dgdx_qp[i] = EpetraExt::ModelEvaluator::Derivative(
133 me->create_DgDx_op(i));
134
135 // dg/dx_dot
136 if (me_outargs.supports(OUT_ARG_DgDx_dot, i).supports(DERIV_TRANS_MV_BY_ROW))
137 dgdx_dot_qp[i] = EpetraExt::ModelEvaluator::Derivative(
138 Teuchos::rcp(new Epetra_MultiVector(
139 *(me->get_x_map()),
140 me->get_g_map(i)->NumGlobalElements())));
141 else if (me_outargs.supports(OUT_ARG_DgDx_dot, i).supports(DERIV_MV_BY_COL))
142 dgdx_dot_qp[i] = EpetraExt::ModelEvaluator::Derivative(
143 Teuchos::rcp(new Epetra_MultiVector(
144 *(me->get_g_map(i)),
145 me->get_x_map()->NumGlobalElements())));
146 else if (me_outargs.supports(OUT_ARG_DgDx_dot, i).supports(DERIV_LINEAR_OP))
147 dgdx_dot_qp[i] = EpetraExt::ModelEvaluator::Derivative(
148 me->create_DgDx_dot_op(i));
149
150 // dg/dp
151 dgdp_qp[i].resize(num_p);
152 for (int j=0; j<num_p; j++)
153 if (me_outargs.supports(OUT_ARG_DgDp, i, j).supports(DERIV_TRANS_MV_BY_ROW))
154 dgdp_qp[i][j] = EpetraExt::ModelEvaluator::Derivative(
155 Teuchos::rcp(new Epetra_MultiVector(
156 *(me->get_p_map(j)),
157 me->get_g_map(i)->NumGlobalElements())));
158 else if (me_outargs.supports(OUT_ARG_DgDp, i, j).supports(DERIV_MV_BY_COL))
159 dgdp_qp[i][j] = EpetraExt::ModelEvaluator::Derivative(
160 Teuchos::rcp(new Epetra_MultiVector(
161 *(me->get_g_map(i)),
162 me->get_p_map(j)->NumGlobalElements())));
163 else if (me_outargs.supports(OUT_ARG_DgDp, i, j).supports(DERIV_LINEAR_OP))
164 dgdp_qp[i][j] = EpetraExt::ModelEvaluator::Derivative(
165 me->create_DgDp_op(i,j));
166 }
167}
168
169// Overridden from EpetraExt::ModelEvaluator
170
171Teuchos::RCP<const Epetra_Map>
173get_x_map() const
174{
175 return me->get_x_map();
176}
177
178Teuchos::RCP<const Epetra_Map>
180get_f_map() const
181{
182 return me->get_f_map();
183}
184
185Teuchos::RCP<const Epetra_Map>
187get_p_map(int l) const
188{
189 return me->get_p_map(l);
190}
191
192Teuchos::RCP<const Epetra_Map>
194get_g_map(int l) const
195{
196 return me->get_g_map(l);
197}
198
199Teuchos::RCP<const Teuchos::Array<std::string> >
201get_p_names(int l) const
202{
203 return me->get_p_names(l);
204}
205
206Teuchos::RCP<const Epetra_Vector>
208get_x_init() const
209{
210 return me->get_x_init();
211}
212
213Teuchos::RCP<const Epetra_Vector>
215get_p_init(int l) const
216{
217 return me->get_p_init(l);
218}
219
220Teuchos::RCP<Epetra_Operator>
222create_W() const
223{
224 return me->create_W();
225}
226
227EpetraExt::ModelEvaluator::InArgs
229createInArgs() const
230{
231 InArgsSetup inArgs;
232 InArgs me_inargs = me->createInArgs();
233
234 inArgs.setModelEvalDescription(this->description());
235 inArgs.set_Np(num_p);
236 inArgs.setSupports(IN_ARG_x_dot, me_inargs.supports(IN_ARG_x_dot));
237 inArgs.setSupports(IN_ARG_x, me_inargs.supports(IN_ARG_x));
238 inArgs.setSupports(IN_ARG_t, me_inargs.supports(IN_ARG_t));
239 inArgs.setSupports(IN_ARG_alpha, me_inargs.supports(IN_ARG_alpha));
240 inArgs.setSupports(IN_ARG_beta, me_inargs.supports(IN_ARG_beta));
241
242 for (int i=0; i<num_p; i++)
243 inArgs.setSupports(IN_ARG_p_sg, i, true);
244 inArgs.setSupports(IN_ARG_x_sg, me_inargs.supports(IN_ARG_x));
245 inArgs.setSupports(IN_ARG_x_dot_sg, me_inargs.supports(IN_ARG_x_dot));
246 inArgs.setSupports(IN_ARG_sg_basis, true);
247 inArgs.setSupports(IN_ARG_sg_quadrature, true);
248
249 return inArgs;
250}
251
252EpetraExt::ModelEvaluator::OutArgs
254createOutArgs() const
255{
256 OutArgsSetup outArgs;
257 OutArgs me_outargs = me->createOutArgs();
258
259 outArgs.setModelEvalDescription(this->description());
260 outArgs.set_Np_Ng(num_p, num_g);
261 outArgs.setSupports(OUT_ARG_f, me_outargs.supports(OUT_ARG_f));
262 outArgs.setSupports(OUT_ARG_W, me_outargs.supports(OUT_ARG_W));
263 for (int j=0; j<num_p; j++)
264 outArgs.setSupports(OUT_ARG_DfDp, j,
265 me_outargs.supports(OUT_ARG_DfDp, j));
266 for (int i=0; i<num_g; i++) {
267 outArgs.setSupports(OUT_ARG_DgDx, i,
268 me_outargs.supports(OUT_ARG_DgDx, i));
269 outArgs.setSupports(OUT_ARG_DgDx_dot, i,
270 me_outargs.supports(OUT_ARG_DgDx_dot, i));
271 for (int j=0; j<num_p; j++)
272 outArgs.setSupports(OUT_ARG_DgDp, i, j,
273 me_outargs.supports(OUT_ARG_DgDp, i, j));
274 }
275
276 outArgs.setSupports(OUT_ARG_f_sg, me_outargs.supports(OUT_ARG_f));
277 if (me_outargs.supports(OUT_ARG_W)) {
278 outArgs.set_W_properties(me_outargs.get_W_properties());
279 outArgs.setSupports(OUT_ARG_W_sg, true);
280 }
281 for (int j=0; j<num_p; j++) {
282 outArgs.setSupports(OUT_ARG_DfDp_sg, j,
283 me_outargs.supports(OUT_ARG_DfDp, j));
284 }
285 for (int i=0; i<num_g; i++) {
286 outArgs.setSupports(OUT_ARG_g_sg, i, true);
287 outArgs.setSupports(OUT_ARG_DgDx_sg, i,
288 me_outargs.supports(OUT_ARG_DgDx, i));
289 outArgs.setSupports(OUT_ARG_DgDx_dot_sg, i,
290 me_outargs.supports(OUT_ARG_DgDx_dot, i));
291 for (int j=0; j<num_p; j++)
292 outArgs.setSupports(OUT_ARG_DgDp_sg, i, j,
293 me_outargs.supports(OUT_ARG_DgDp, i, j));
294 }
295
296 return outArgs;
297}
298
299void
301evalModel(const InArgs& inArgs, const OutArgs& outArgs) const
302{
303 // Create underlying inargs
304 InArgs me_inargs = me->createInArgs();
305 if (me_inargs.supports(IN_ARG_x))
306 me_inargs.set_x(inArgs.get_x());
307 if (me_inargs.supports(IN_ARG_x_dot))
308 me_inargs.set_x_dot(inArgs.get_x_dot());
309 if (me_inargs.supports(IN_ARG_alpha))
310 me_inargs.set_alpha(inArgs.get_alpha());
311 if (me_inargs.supports(IN_ARG_beta))
312 me_inargs.set_beta(inArgs.get_beta());
313 if (me_inargs.supports(IN_ARG_t))
314 me_inargs.set_t(inArgs.get_t());
315 for (int i=0; i<num_p; i++)
316 me_inargs.set_p(i, inArgs.get_p(i));
317
318 // Create underlying outargs
319 OutArgs me_outargs = me->createOutArgs();
320 if (me_outargs.supports(OUT_ARG_f))
321 me_outargs.set_f(outArgs.get_f());
322 if (me_outargs.supports(OUT_ARG_W))
323 me_outargs.set_W(outArgs.get_W());
324 for (int j=0; j<num_p; j++)
325 if (!outArgs.supports(OUT_ARG_DfDp, j).none())
326 me_outargs.set_DfDp(j, outArgs.get_DfDp(j));
327 for (int i=0; i<num_g; i++) {
328 me_outargs.set_g(i, outArgs.get_g(i));
329 if (!outArgs.supports(OUT_ARG_DgDx, i).none())
330 me_outargs.set_DgDx(i, outArgs.get_DgDx(i));
331 if (!outArgs.supports(OUT_ARG_DgDx_dot, i).none())
332 me_outargs.set_DgDx(i, outArgs.get_DgDx_dot(i));
333 for (int j=0; j<num_p; j++)
334 if (!outArgs.supports(OUT_ARG_DgDp, i, j).none())
335 me_outargs.set_DgDp(i, j, outArgs.get_DgDp(i,j));
336 }
337
338 bool do_quad = false;
339 InArgs::sg_const_vector_t x_sg;
340 InArgs::sg_const_vector_t x_dot_sg;
341 Teuchos::Array<InArgs::sg_const_vector_t> p_sg(num_p);
342 OutArgs::sg_vector_t f_sg;
343 OutArgs::sg_operator_t W_sg;
344 Teuchos::Array<SGDerivative> dfdp_sg(num_p);
345 Teuchos::Array<OutArgs::sg_vector_t> g_sg(num_g);
346 Teuchos::Array<SGDerivative> dgdx_sg(num_g);
347 Teuchos::Array<SGDerivative> dgdx_dot_sg(num_g);
348 Teuchos::Array< Teuchos::Array<SGDerivative> > dgdp_sg(num_g);
349 TEUCHOS_TEST_FOR_EXCEPTION(inArgs.get_sg_basis() == Teuchos::null,
350 std::logic_error,
351 "Error! Stokhos::SGQuadModelEvaluator::evalModel(): " <<
352 "SG basis inArg cannot be null!");
353 TEUCHOS_TEST_FOR_EXCEPTION(inArgs.get_sg_quadrature() == Teuchos::null,
354 std::logic_error,
355 "Error! Stokhos::SGQuadModelEvaluator::evalModel(): " <<
356 "SG quadrature inArg cannot be null!");
357 Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> > basis =
358 inArgs.get_sg_basis();
359 Teuchos::RCP< const Stokhos::Quadrature<int,double> > quad =
360 inArgs.get_sg_quadrature();
361 if (inArgs.supports(IN_ARG_x_sg)) {
362 x_sg = inArgs.get_x_sg();
363 if (x_sg != Teuchos::null) {
364 do_quad = true;
365 }
366 }
367 if (inArgs.supports(IN_ARG_x_dot_sg)) {
368 x_dot_sg = inArgs.get_x_dot_sg();
369 if (x_dot_sg != Teuchos::null) {
370 do_quad = true;
371 }
372 }
373 for (int i=0; i<num_p; i++) {
374 p_sg[i] = inArgs.get_p_sg(i);
375 if (p_sg[i] != Teuchos::null) {
376 do_quad = true;
377 }
378 }
379 if (outArgs.supports(OUT_ARG_f_sg)) {
380 f_sg = outArgs.get_f_sg();
381 if (f_sg != Teuchos::null)
382 f_sg->init(0.0);
383 }
384 if (outArgs.supports(OUT_ARG_W_sg)) {
385 W_sg = outArgs.get_W_sg();
386 if (W_sg != Teuchos::null)
387 W_sg->init(0.0);
388 }
389 for (int i=0; i<num_p; i++) {
390 if (!outArgs.supports(OUT_ARG_DfDp_sg, i).none()) {
391 dfdp_sg[i] = outArgs.get_DfDp_sg(i);
392 if (dfdp_sg[i].getMultiVector() != Teuchos::null)
393 dfdp_sg[i].getMultiVector()->init(0.0);
394 else if (dfdp_sg[i].getLinearOp() != Teuchos::null)
395 dfdp_sg[i].getLinearOp()->init(0.0);
396 }
397 }
398
399 for (int i=0; i<num_g; i++) {
400 g_sg[i] = outArgs.get_g_sg(i);
401 if (g_sg[i] != Teuchos::null)
402 g_sg[i]->init(0.0);
403
404 if (!outArgs.supports(OUT_ARG_DgDx_sg, i).none()) {
405 dgdx_sg[i] = outArgs.get_DgDx_sg(i);
406 if (dgdx_sg[i].getMultiVector() != Teuchos::null)
407 dgdx_sg[i].getMultiVector()->init(0.0);
408 else if (dgdx_sg[i].getLinearOp() != Teuchos::null)
409 dgdx_sg[i].getLinearOp()->init(0.0);
410 }
411
412 if (!outArgs.supports(OUT_ARG_DgDx_dot_sg, i).none()) {
413 dgdx_dot_sg[i] = outArgs.get_DgDx_dot_sg(i);
414 if (dgdx_dot_sg[i].getMultiVector() != Teuchos::null)
415 dgdx_dot_sg[i].getMultiVector()->init(0.0);
416 else if (dgdx_dot_sg[i].getLinearOp() != Teuchos::null)
417 dgdx_dot_sg[i].getLinearOp()->init(0.0);
418 }
419
420 dgdp_sg[i].resize(num_p);
421 for (int j=0; j<num_p; j++) {
422 if (!outArgs.supports(OUT_ARG_DgDp_sg, i, j).none()) {
423 dgdp_sg[i][j] = outArgs.get_DgDp_sg(i,j);
424 if (dgdp_sg[i][j].getMultiVector() != Teuchos::null)
425 dgdp_sg[i][j].getMultiVector()->init(0.0);
426 else if (dgdp_sg[i][j].getLinearOp() != Teuchos::null)
427 dgdp_sg[i][j].getLinearOp()->init(0.0);
428 }
429 }
430 }
431
432 if (do_quad) {
433 // Get quadrature data
434 const Teuchos::Array< Teuchos::Array<double> >& quad_points =
435 quad->getQuadPoints();
436 const Teuchos::Array<double>& quad_weights =
437 quad->getQuadWeights();
438 const Teuchos::Array< Teuchos::Array<double> > & quad_values =
439 quad->getBasisAtQuadPoints();
440 const Teuchos::Array<double>& basis_norms = basis->norm_squared();
441
442 // Perform integrations
443 for (int qp=0; qp<quad_points.size(); qp++) {
444
445 // StieltjesPCEBasis can introduce quadrature points with zero weight
446 // Don't do those evaluations, since the model might not like the
447 // quadrature points (i.e., zero)
448 if (quad_weights[qp] == 0.0)
449 continue;
450
451 {
452#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
453 TEUCHOS_FUNC_TIME_MONITOR_DIFF("Stokhos: SGQuadModelEvaluator -- Polynomial Evaluation",
454 PolyEvaluation);
455#endif
456
457 // Evaluate inputs at quadrature points
458 if (x_sg != Teuchos::null) {
459#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
460 TEUCHOS_FUNC_TIME_MONITOR("Stokhos: SGQuadModelEvaluator -- X Evaluation");
461#endif
462 x_sg->evaluate(quad_values[qp], *x_qp);
463 me_inargs.set_x(x_qp);
464 }
465 if (x_dot_sg != Teuchos::null) {
466#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
467 TEUCHOS_FUNC_TIME_MONITOR("Stokhos: SGQuadModelEvaluator -- X_dot Evaluation");
468#endif
469 x_dot_sg->evaluate(quad_values[qp], *x_dot_qp);
470 me_inargs.set_x_dot(x_qp);
471 }
472 for (int i=0; i<num_p; i++) {
473 if (p_sg[i] != Teuchos::null) {
474#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
475 TEUCHOS_FUNC_TIME_MONITOR("Stokhos: SGQuadModelEvaluator -- P Evaluation");
476#endif
477 p_sg[i]->evaluate(quad_values[qp], *(p_qp[i]));
478 me_inargs.set_p(i, p_qp[i]);
479 }
480 }
481 if (f_sg != Teuchos::null)
482 me_outargs.set_f(f_qp);
483 if (W_sg != Teuchos::null)
484 me_outargs.set_W(W_qp);
485 for (int i=0; i<num_p; i++) {
486 if (!dfdp_sg[i].isEmpty())
487 me_outargs.set_DfDp(i, dfdp_qp[i]);
488 }
489 for (int i=0; i<num_g; i++) {
490 if (g_sg[i] != Teuchos::null)
491 me_outargs.set_g(i, g_qp[i]);
492 if (!dgdx_dot_sg[i].isEmpty())
493 me_outargs.set_DgDx_dot(i, dgdx_dot_qp[i]);
494 if (!dgdx_sg[i].isEmpty())
495 me_outargs.set_DgDx(i, dgdx_qp[i]);
496 for (int j=0; j<num_p; j++)
497 if (!dgdp_sg[i][j].isEmpty())
498 me_outargs.set_DgDp(i, j, dgdp_qp[i][j]);
499 }
500
501 }
502
503 {
504#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
505 TEUCHOS_FUNC_TIME_MONITOR("Stokhos: SGQuadModelEvaluator -- Model Evaluation");
506#endif
507
508 // Evaluate model at quadrature points
509 me->evalModel(me_inargs, me_outargs);
510
511 }
512
513 {
514#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
515 TEUCHOS_FUNC_TIME_MONITOR_DIFF(
516 "Stokhos: SGQuadModelEvaluator -- Polynomial Integration", Integration);
517#endif
518
519 // Sum in results
520 if (f_sg != Teuchos::null) {
521#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
522 TEUCHOS_FUNC_TIME_MONITOR("Stokhos: SGQuadModelEvaluator -- F Integration");
523#endif
524 f_sg->sumIntoAllTerms(quad_weights[qp], quad_values[qp], basis_norms,
525 *f_qp);
526 }
527 if (W_sg != Teuchos::null) {
528#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
529 TEUCHOS_FUNC_TIME_MONITOR("Stokhos: SGQuadModelEvaluator -- W Integration");
530#endif
531 W_sg->sumIntoAllTerms(quad_weights[qp], quad_values[qp], basis_norms,
532 *W_qp);
533 }
534 for (int j=0; j<num_p; j++) {
535 if (!dfdp_sg[j].isEmpty()) {
536#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
537 TEUCHOS_FUNC_TIME_MONITOR(
538 "Stokhos: SGQuadModelEvaluator -- df/dp Integration");
539#endif
540 if (dfdp_sg[j].getMultiVector() != Teuchos::null) {
541 dfdp_sg[j].getMultiVector()->sumIntoAllTerms(
542 quad_weights[qp], quad_values[qp], basis_norms,
543 *(dfdp_qp[j].getMultiVector()));
544 }
545 else if (dfdp_sg[j].getLinearOp() != Teuchos::null) {
546 dfdp_sg[j].getLinearOp()->sumIntoAllTerms(
547 quad_weights[qp], quad_values[qp], basis_norms,
548 *(dfdp_qp[j].getLinearOp()));
549 }
550 }
551 }
552 for (int i=0; i<num_g; i++) {
553 if (g_sg[i] != Teuchos::null) {
554#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
555 TEUCHOS_FUNC_TIME_MONITOR("Stokhos: SGQuadModelEvaluator -- G Integration");
556#endif
557 g_sg[i]->sumIntoAllTerms(quad_weights[qp], quad_values[qp],
558 basis_norms, *g_qp[i]);
559 }
560 if (!dgdx_dot_sg[i].isEmpty()) {
561#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
562 TEUCHOS_FUNC_TIME_MONITOR(
563 "Stokhos: SGQuadModelEvaluator -- dg/dx_dot Integration");
564#endif
565 if (dgdx_dot_sg[i].getMultiVector() != Teuchos::null) {
566 dgdx_dot_sg[i].getMultiVector()->sumIntoAllTerms(
567 quad_weights[qp], quad_values[qp], basis_norms,
568 *(dgdx_dot_qp[i].getMultiVector()));
569 }
570 else if (dgdx_dot_sg[i].getLinearOp() != Teuchos::null) {
571 dgdx_dot_sg[i].getLinearOp()->sumIntoAllTerms(
572 quad_weights[qp], quad_values[qp], basis_norms,
573 *(dgdx_dot_qp[i].getLinearOp()));
574 }
575 }
576 if (!dgdx_sg[i].isEmpty()) {
577#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
578 TEUCHOS_FUNC_TIME_MONITOR(
579 "Stokhos: SGQuadModelEvaluator -- dg/dx Integration");
580#endif
581 if (dgdx_sg[i].getMultiVector() != Teuchos::null) {
582 dgdx_sg[i].getMultiVector()->sumIntoAllTerms(
583 quad_weights[qp], quad_values[qp], basis_norms,
584 *(dgdx_qp[i].getMultiVector()));
585 }
586 else if (dgdx_sg[i].getLinearOp() != Teuchos::null) {
587 dgdx_sg[i].getLinearOp()->sumIntoAllTerms(
588 quad_weights[qp], quad_values[qp], basis_norms,
589 *(dgdx_qp[i].getLinearOp()));
590 }
591 }
592 for (int j=0; j<num_p; j++) {
593 if (!dgdp_sg[i][j].isEmpty()) {
594#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
595 TEUCHOS_FUNC_TIME_MONITOR(
596 "Stokhos: SGQuadModelEvaluator -- dg/dp Integration");
597#endif
598 if (dgdp_sg[i][j].getMultiVector() != Teuchos::null) {
599 dgdp_sg[i][j].getMultiVector()->sumIntoAllTerms(
600 quad_weights[qp], quad_values[qp], basis_norms,
601 *(dgdp_qp[i][j].getMultiVector()));
602 }
603 else if (dgdp_sg[i][j].getLinearOp() != Teuchos::null) {
604 dgdp_sg[i][j].getLinearOp()->sumIntoAllTerms(
605 quad_weights[qp], quad_values[qp], basis_norms,
606 *(dgdp_qp[i][j].getLinearOp()));
607 }
608 }
609 }
610 }
611
612 }
613 }
614 }
615 else {
616 // Compute the non-SG functions
617 me->evalModel(me_inargs, me_outargs);
618 }
619}
Teuchos::RCP< const Epetra_Map > get_f_map() const
Return residual vector map.
Teuchos::RCP< Epetra_Operator > W_qp
W operator.
Teuchos::RCP< const Epetra_Map > get_g_map(int l) const
Return observation vector map.
Teuchos::RCP< EpetraExt::ModelEvaluator > me
Underlying model evaluator.
Teuchos::RCP< Epetra_Vector > f_qp
Residual vector.
Teuchos::RCP< Epetra_Vector > x_dot_qp
Time derivative vector.
void evalModel(const InArgs &inArgs, const OutArgs &outArgs) const
Evaluate model on InArgs.
Teuchos::RCP< const Epetra_Vector > get_x_init() const
Return initial solution.
SGQuadModelEvaluator(const Teuchos::RCP< EpetraExt::ModelEvaluator > &me)
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int l) const
Return array of parameter names.
Teuchos::Array< Teuchos::RCP< Epetra_Vector > > p_qp
Parameter vectors.
Teuchos::Array< Teuchos::Array< EpetraExt::ModelEvaluator::Derivative > > dgdp_qp
Response sensitivities.
Teuchos::RCP< Epetra_Vector > x_qp
Solution vector.
Teuchos::Array< EpetraExt::ModelEvaluator::Derivative > dgdx_dot_qp
Response derivative.
Teuchos::RCP< const Epetra_Map > get_p_map(int l) const
Return parameter vector map.
int num_p
Number of parameter vectors.
OutArgs createOutArgs() const
Create OutArgs.
Teuchos::Array< Teuchos::RCP< Epetra_Vector > > g_qp
Response vectors.
Teuchos::Array< EpetraExt::ModelEvaluator::Derivative > dgdx_qp
Response derivative.
Teuchos::Array< EpetraExt::ModelEvaluator::Derivative > dfdp_qp
Residual derivatives.
Teuchos::RCP< const Epetra_Vector > get_p_init(int l) const
Return initial parameters.
Teuchos::RCP< const Epetra_Map > get_x_map() const
Return solution vector map.
Teuchos::RCP< Epetra_Operator > create_W() const
Create W = alpha*M + beta*J matrix.
int num_g
Number of response vectors.
InArgs createInArgs() const
Create InArgs.