Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_SGQuadMPModelEvaluator.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
51#include "Epetra_Map.h"
52#include "Epetra_Vector.h"
53#include "Teuchos_TimeMonitor.hpp"
54#include "Teuchos_Assert.hpp"
55
58 const Teuchos::RCP<EpetraExt::ModelEvaluator>& me_,
59 const Teuchos::RCP<const EpetraExt::MultiComm>& mp_comm_,
60 const Teuchos::RCP<const Epetra_Map>& mp_block_map_) :
61 me(me_),
62 mp_comm(mp_comm_),
63 mp_block_map(mp_block_map_),
64 num_p(0),
65 num_g(0),
66 num_p_mp(0),
67 num_g_mp(0),
68 mp_p_index_map(),
69 mp_g_index_map(),
70 x_dot_mp(),
71 x_mp(),
72 p_mp(),
73 f_mp(),
74 W_mp(),
75 dfdp_mp(),
76 g_mp(),
77 dgdx_mp(),
78 dgdx_dot_mp(),
79 dgdp_mp()
80{
81 int num_qp = mp_block_map->NumMyElements();
82
83 Teuchos::RCP<const Epetra_Map> x_map;
84 Teuchos::RCP<const Epetra_Map> f_map;
85
86 // Create storage for x_dot, x, and p at a quad point
87 InArgs me_inargs = me->createInArgs();
88 if (me_inargs.supports(IN_ARG_x_dot)) {
89 x_map = me->get_x_map();
90 x_dot_mp = Teuchos::rcp(new ProductEpetraVector(
91 mp_block_map, x_map, mp_comm));
92 }
93 if (me_inargs.supports(IN_ARG_x)) {
94 x_map = me->get_x_map();
95 x_mp = Teuchos::rcp(new ProductEpetraVector(
96 mp_block_map, me->get_x_map(), mp_comm));
97 }
98
99 // Get the p_mp's supported and build index map
100 num_p = me_inargs.Np();
101 for (int i=0; i<num_p; i++) {
102 if (me_inargs.supports(IN_ARG_p_mp, i))
103 mp_p_index_map.push_back(i);
104 }
105 num_p_mp = mp_p_index_map.size();
106
107 p_mp.resize(num_p_mp);
108 for (int i=0; i<num_p_mp; i++)
109 p_mp[i] = Teuchos::rcp(new ProductEpetraVector(
110 mp_block_map, me->get_p_map(mp_p_index_map[i]),
111 mp_comm));
112
113 // Create storage for f and W at a quad point
114 OutArgs me_outargs = me->createOutArgs();
115
116 // f
117 if (me_outargs.supports(OUT_ARG_f)) {
118 f_map = me->get_f_map();
119 f_mp = Teuchos::rcp(new ProductEpetraVector(mp_block_map, f_map, mp_comm));
120 }
121
122 // W
123 if (me_outargs.supports(OUT_ARG_W)) {
124 W_mp = Teuchos::rcp(new ProductEpetraOperator(mp_block_map, x_map, f_map,
125 mp_comm));
126 for (int i=0; i<num_qp; i++)
127 W_mp->setCoeffPtr(i, me->create_W());
128 }
129
130 // df/dp -- note we potentially support derivatives w.r.t. all parameters,
131 // not just those for which p_mp is supported
132 dfdp_mp.resize(num_p);
133 for (int i=0; i<num_p; i++) {
134 Teuchos::RCP<const Epetra_Map> p_map = me->get_p_map(i);
135 DerivativeSupport ds = me_outargs.supports(OUT_ARG_DfDp_mp,i);
136 if (ds.supports(DERIV_MV_BY_COL))
137 dfdp_mp[i] = EpetraExt::ModelEvaluator::MPDerivative(
138 Teuchos::rcp(new ProductEpetraMultiVector(
139 mp_block_map, f_map, mp_comm,
140 p_map->NumGlobalElements())));
141 else if (ds.supports(DERIV_TRANS_MV_BY_ROW))
142 dfdp_mp[i] = EpetraExt::ModelEvaluator::MPDerivative(
143 Teuchos::rcp(new ProductEpetraMultiVector(
144 mp_block_map, p_map, mp_comm,
145 f_map->NumGlobalElements())));
146 else if (ds.supports(DERIV_LINEAR_OP)) {
147 Teuchos::RCP< ProductEpetraOperator > dfdp_mp_op =
148 Teuchos::rcp(new ProductEpetraOperator(
149 mp_block_map, p_map, f_map,
150 mp_comm));
151 for (int j=0; j<num_qp; j++)
152 dfdp_mp_op->setCoeffPtr(j, me->create_DfDp_op(i));
153 dfdp_mp[i] = EpetraExt::ModelEvaluator::MPDerivative(dfdp_mp_op);
154 }
155
156 }
157
158 // Get the g_mp's supported and build index map
159 num_g = me_outargs.Ng();
160 for (int i=0; i<num_g; i++) {
161 if (me_outargs.supports(OUT_ARG_g_mp, i))
162 mp_g_index_map.push_back(i);
163 }
164 num_g_mp = mp_g_index_map.size();
165
166 g_mp.resize(num_g_mp);
167 dgdx_mp.resize(num_g_mp);
168 dgdx_dot_mp.resize(num_g_mp);
169 dgdp_mp.resize(num_g_mp);
170 for (int i=0; i<num_g_mp; i++) {
171 int ii = mp_g_index_map[i];
172 Teuchos::RCP<const Epetra_Map> g_map = me->get_g_map(ii);
173
174 // g
175 g_mp[i] =
176 Teuchos::rcp(new ProductEpetraVector(mp_block_map, g_map, mp_comm));
177
178 // dg/dx
179 DerivativeSupport ds_x = me_outargs.supports(OUT_ARG_DgDx_mp, ii);
180 if (ds_x.supports(DERIV_TRANS_MV_BY_ROW))
181 dgdx_mp[i] = EpetraExt::ModelEvaluator::MPDerivative(
182 Teuchos::rcp(new ProductEpetraMultiVector(mp_block_map, x_map, mp_comm,
183 g_map->NumGlobalElements())));
184 else if (ds_x.supports(DERIV_MV_BY_COL))
185 dgdx_mp[i] = EpetraExt::ModelEvaluator::MPDerivative(
186 Teuchos::rcp(new ProductEpetraMultiVector(mp_block_map, g_map, mp_comm,
187 x_map->NumGlobalElements())));
188 else if (ds_x.supports(DERIV_LINEAR_OP)) {
189 Teuchos::RCP<ProductEpetraOperator> dgdx_mp_op =
190 Teuchos::rcp(new ProductEpetraOperator(mp_block_map, x_map, g_map,
191 mp_comm));
192 for (int j=0; j<num_qp; j++)
193 dgdx_mp_op->setCoeffPtr(j, me->create_DgDx_op(ii));
194 dgdx_mp[i] = EpetraExt::ModelEvaluator::MPDerivative(dgdx_mp_op);
195 }
196
197 // dg/dx_dot
198 DerivativeSupport ds_xdot = me_outargs.supports(OUT_ARG_DgDx_dot_mp, ii);
199 if (ds_xdot.supports(DERIV_TRANS_MV_BY_ROW))
200 dgdx_dot_mp[i] = EpetraExt::ModelEvaluator::MPDerivative(
201 Teuchos::rcp(new ProductEpetraMultiVector(mp_block_map, x_map, mp_comm,
202 g_map->NumGlobalElements())));
203 else if (ds_xdot.supports(DERIV_MV_BY_COL))
204 dgdx_dot_mp[i] = EpetraExt::ModelEvaluator::MPDerivative(
205 Teuchos::rcp(new ProductEpetraMultiVector(mp_block_map, g_map, mp_comm,
206 x_map->NumGlobalElements())));
207 else if (ds_xdot.supports(DERIV_LINEAR_OP)) {
208 Teuchos::RCP< ProductEpetraOperator > dgdx_dot_mp_op =
209 Teuchos::rcp(new ProductEpetraOperator(mp_block_map, x_map, g_map,
210 mp_comm));
211 for (int j=0; j<num_qp; j++)
212 dgdx_dot_mp_op->setCoeffPtr(j, me->create_DgDx_dot_op(ii));
213 dgdx_dot_mp[i] = EpetraExt::ModelEvaluator::MPDerivative(dgdx_dot_mp_op);
214 }
215
216 // dg/dp -- note we potentially support derivatives w.r.t. all parameters,
217 // not just those for which p_mp is supported
218 dgdp_mp[i].resize(num_p);
219 for (int j=0; j<num_p; j++) {
220 Teuchos::RCP<const Epetra_Map> p_map = me->get_p_map(j);
221 DerivativeSupport ds_p = me_outargs.supports(OUT_ARG_DgDp_mp, ii, j);
222 if (ds_p.supports(DERIV_TRANS_MV_BY_ROW))
223 dgdp_mp[i][j] = EpetraExt::ModelEvaluator::MPDerivative(
224 Teuchos::rcp(new ProductEpetraMultiVector(
225 mp_block_map, p_map, mp_comm,
226 g_map->NumGlobalElements())));
227 else if (ds_p.supports(DERIV_MV_BY_COL))
228 dgdp_mp[i][j] = EpetraExt::ModelEvaluator::MPDerivative(
229 Teuchos::rcp(new ProductEpetraMultiVector(
230 mp_block_map, g_map, mp_comm,
231 p_map->NumGlobalElements())));
232 else if (ds_p.supports(DERIV_LINEAR_OP)) {
233 Teuchos::RCP<ProductEpetraOperator> dgdp_mp_op =
234 Teuchos::rcp(new ProductEpetraOperator(mp_block_map, p_map, g_map,
235 mp_comm));
236 for (int k=0; k<num_qp; k++)
237 dgdp_mp_op->setCoeffPtr(k, me->create_DgDp_op(ii,j));
238 dgdp_mp[i][j] = EpetraExt::ModelEvaluator::MPDerivative(dgdp_mp_op);
239 }
240 }
241 }
242}
243
244// Overridden from EpetraExt::ModelEvaluator
245
246Teuchos::RCP<const Epetra_Map>
248get_x_map() const
249{
250 return me->get_x_map();
251}
252
253Teuchos::RCP<const Epetra_Map>
255get_f_map() const
256{
257 return me->get_f_map();
258}
259
260Teuchos::RCP<const Epetra_Map>
262get_p_map(int l) const
263{
264 return me->get_p_map(l);
265}
266
267Teuchos::RCP<const Epetra_Map>
269get_g_map(int l) const
270{
271 return me->get_g_map(l);
272}
273
274Teuchos::RCP<const Teuchos::Array<std::string> >
276get_p_names(int l) const
277{
278 return me->get_p_names(l);
279}
280
281Teuchos::RCP<const Epetra_Vector>
283get_x_init() const
284{
285 return me->get_x_init();
286}
287
288Teuchos::RCP<const Epetra_Vector>
290get_p_init(int l) const
291{
292 return me->get_p_init(l);
293}
294
295Teuchos::RCP<Epetra_Operator>
297create_W() const
298{
299 return me->create_W();
300}
301
302EpetraExt::ModelEvaluator::InArgs
304createInArgs() const
305{
306 InArgsSetup inArgs;
307 InArgs me_inargs = me->createInArgs();
308
309 inArgs.setModelEvalDescription(this->description());
310 inArgs.set_Np(num_p);
311 inArgs.setSupports(IN_ARG_x_dot, me_inargs.supports(IN_ARG_x_dot));
312 inArgs.setSupports(IN_ARG_x, me_inargs.supports(IN_ARG_x));
313 inArgs.setSupports(IN_ARG_t, me_inargs.supports(IN_ARG_t));
314 inArgs.setSupports(IN_ARG_alpha, me_inargs.supports(IN_ARG_alpha));
315 inArgs.setSupports(IN_ARG_beta, me_inargs.supports(IN_ARG_beta));
316
317 for (int i=0; i<num_p_mp; i++)
318 inArgs.setSupports(IN_ARG_p_sg, mp_p_index_map[i], true);
319 inArgs.setSupports(IN_ARG_x_sg, me_inargs.supports(IN_ARG_x));
320 inArgs.setSupports(IN_ARG_x_dot_sg, me_inargs.supports(IN_ARG_x_dot));
321 inArgs.setSupports(IN_ARG_sg_basis, true);
322 inArgs.setSupports(IN_ARG_sg_quadrature, true);
323
324 return inArgs;
325}
326
327EpetraExt::ModelEvaluator::OutArgs
329createOutArgs() const
330{
331 OutArgsSetup outArgs;
332 OutArgs me_outargs = me->createOutArgs();
333
334 outArgs.setModelEvalDescription(this->description());
335 outArgs.set_Np_Ng(num_p, num_g);
336 outArgs.setSupports(OUT_ARG_f, me_outargs.supports(OUT_ARG_f));
337 outArgs.setSupports(OUT_ARG_W, me_outargs.supports(OUT_ARG_W));
338 for (int j=0; j<num_p; j++)
339 outArgs.setSupports(OUT_ARG_DfDp, j,
340 me_outargs.supports(OUT_ARG_DfDp, j));
341 for (int i=0; i<num_g; i++) {
342 outArgs.setSupports(OUT_ARG_DgDx, i,
343 me_outargs.supports(OUT_ARG_DgDx, i));
344 outArgs.setSupports(OUT_ARG_DgDx_dot, i,
345 me_outargs.supports(OUT_ARG_DgDx_dot, i));
346 for (int j=0; j<num_p; j++)
347 outArgs.setSupports(OUT_ARG_DgDp, i, j,
348 me_outargs.supports(OUT_ARG_DgDp, i, j));
349 }
350
351 outArgs.setSupports(OUT_ARG_f_sg, me_outargs.supports(OUT_ARG_f_mp));
352 if (me_outargs.supports(OUT_ARG_W_mp)) {
353 outArgs.set_W_properties(me_outargs.get_W_properties());
354 outArgs.setSupports(OUT_ARG_W_sg, true);
355 }
356 for (int i=0; i<num_g_mp; i++)
357 outArgs.setSupports(OUT_ARG_g_sg, mp_g_index_map[i], true);
358 for (int j=0; j<num_p; j++)
359 outArgs.setSupports(OUT_ARG_DfDp_sg, j,
360 me_outargs.supports(OUT_ARG_DfDp_mp, j));
361 for (int i=0; i<num_g_mp; i++) {
362 int ii = mp_g_index_map[i];
363 outArgs.setSupports(OUT_ARG_DgDx_sg, ii,
364 me_outargs.supports(OUT_ARG_DgDx_mp, ii));
365 outArgs.setSupports(OUT_ARG_DgDx_dot_sg, ii,
366 me_outargs.supports(OUT_ARG_DgDx_dot_mp, ii));
367 for (int j=0; j<num_p; j++)
368 outArgs.setSupports(OUT_ARG_DgDp_sg, ii, j,
369 me_outargs.supports(OUT_ARG_DgDp_mp, ii, j));
370 }
371
372 return outArgs;
373}
374
375void
377evalModel(const InArgs& inArgs, const OutArgs& outArgs) const
378{
379 // Create underlying inargs
380 InArgs me_inargs = me->createInArgs();
381 if (me_inargs.supports(IN_ARG_x))
382 me_inargs.set_x(inArgs.get_x());
383 if (me_inargs.supports(IN_ARG_x_dot))
384 me_inargs.set_x_dot(inArgs.get_x_dot());
385 if (me_inargs.supports(IN_ARG_alpha))
386 me_inargs.set_alpha(inArgs.get_alpha());
387 if (me_inargs.supports(IN_ARG_beta))
388 me_inargs.set_beta(inArgs.get_beta());
389 if (me_inargs.supports(IN_ARG_t))
390 me_inargs.set_t(inArgs.get_t());
391 for (int i=0; i<num_p; i++)
392 me_inargs.set_p(i, inArgs.get_p(i));
393
394 // Create underlying outargs
395 OutArgs me_outargs = me->createOutArgs();
396 if (me_outargs.supports(OUT_ARG_f))
397 me_outargs.set_f(outArgs.get_f());
398 if (me_outargs.supports(OUT_ARG_W))
399 me_outargs.set_W(outArgs.get_W());
400 for (int j=0; j<num_p; j++)
401 if (!outArgs.supports(OUT_ARG_DfDp, j).none())
402 me_outargs.set_DfDp(j, outArgs.get_DfDp(j));
403 for (int i=0; i<num_g; i++) {
404 me_outargs.set_g(i, outArgs.get_g(i));
405 if (!outArgs.supports(OUT_ARG_DgDx, i).none())
406 me_outargs.set_DgDx(i, outArgs.get_DgDx(i));
407 if (!outArgs.supports(OUT_ARG_DgDx_dot, i).none())
408 me_outargs.set_DgDx(i, outArgs.get_DgDx_dot(i));
409 for (int j=0; j<outArgs.Np(); j++)
410 if (!outArgs.supports(OUT_ARG_DgDp, i, j).none())
411 me_outargs.set_DgDp(i, j, outArgs.get_DgDp(i,j));
412 }
413
414 bool do_quad = false;
415 InArgs::sg_const_vector_t x_sg;
416 InArgs::sg_const_vector_t x_dot_sg;
417 Teuchos::Array<InArgs::sg_const_vector_t> p_sg(num_p_mp);
418 OutArgs::sg_vector_t f_sg;
419 OutArgs::sg_operator_t W_sg;
420 Teuchos::Array<SGDerivative> dfdp_sg(num_p);
421 Teuchos::Array<OutArgs::sg_vector_t> g_sg(num_g_mp);
422 Teuchos::Array<SGDerivative> dgdx_sg(num_g_mp);
423 Teuchos::Array<SGDerivative> dgdx_dot_sg(num_g_mp);
424 Teuchos::Array< Teuchos::Array<SGDerivative> > dgdp_sg(num_g_mp);
425 TEUCHOS_TEST_FOR_EXCEPTION(inArgs.get_sg_basis() == Teuchos::null,
426 std::logic_error,
427 "Error! Stokhos::SGQuadModelEvaluator::evalModel(): " <<
428 "SG basis inArg cannot be null!");
429 TEUCHOS_TEST_FOR_EXCEPTION(inArgs.get_sg_quadrature() == Teuchos::null,
430 std::logic_error,
431 "Error! Stokhos::SGQuadModelEvaluator::evalModel(): " <<
432 "SG quadrature inArg cannot be null!");
433 Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> > basis =
434 inArgs.get_sg_basis();
435 Teuchos::RCP< const Stokhos::Quadrature<int,double> > quad =
436 inArgs.get_sg_quadrature();
437
438 if (inArgs.supports(IN_ARG_x_sg)) {
439 x_sg = inArgs.get_x_sg();
440 if (x_sg != Teuchos::null) {
441 do_quad = true;
442 }
443 }
444 if (inArgs.supports(IN_ARG_x_dot_sg)) {
445 x_dot_sg = inArgs.get_x_dot_sg();
446 if (x_dot_sg != Teuchos::null) {
447 do_quad = true;
448 }
449 }
450 for (int i=0; i<num_p_mp; i++) {
451 p_sg[i] = inArgs.get_p_sg(mp_p_index_map[i]);
452 if (p_sg[i] != Teuchos::null) {
453 do_quad = true;
454 }
455 }
456 if (outArgs.supports(OUT_ARG_f_sg)) {
457 f_sg = outArgs.get_f_sg();
458 if (f_sg != Teuchos::null)
459 f_sg->init(0.0);
460 }
461 if (outArgs.supports(OUT_ARG_W_sg)) {
462 W_sg = outArgs.get_W_sg();
463 if (W_sg != Teuchos::null)
464 W_sg->init(0.0);
465 }
466 for (int i=0; i<num_p; i++) {
467 if (!outArgs.supports(OUT_ARG_DfDp_sg, i).none()) {
468 dfdp_sg[i] = outArgs.get_DfDp_sg(i);
469 if (dfdp_sg[i].getMultiVector() != Teuchos::null)
470 dfdp_sg[i].getMultiVector()->init(0.0);
471 else if (dfdp_sg[i].getLinearOp() != Teuchos::null)
472 dfdp_sg[i].getLinearOp()->init(0.0);
473 }
474 }
475
476 for (int i=0; i<num_g_mp; i++) {
477 int ii = mp_g_index_map[i];
478 g_sg[i] = outArgs.get_g_sg(ii);
479 if (g_sg[i] != Teuchos::null)
480 g_sg[i]->init(0.0);
481
482 if (!outArgs.supports(OUT_ARG_DgDx_sg, ii).none()) {
483 dgdx_sg[i] = outArgs.get_DgDx_sg(ii);
484 if (dgdx_sg[i].getMultiVector() != Teuchos::null)
485 dgdx_sg[i].getMultiVector()->init(0.0);
486 else if (dgdx_sg[i].getLinearOp() != Teuchos::null)
487 dgdx_sg[i].getLinearOp()->init(0.0);
488 }
489
490 if (!outArgs.supports(OUT_ARG_DgDx_dot_sg, ii).none()) {
491 dgdx_dot_sg[i] = outArgs.get_DgDx_dot_sg(ii);
492 if (dgdx_dot_sg[i].getMultiVector() != Teuchos::null)
493 dgdx_dot_sg[i].getMultiVector()->init(0.0);
494 else if (dgdx_dot_sg[i].getLinearOp() != Teuchos::null)
495 dgdx_dot_sg[i].getLinearOp()->init(0.0);
496 }
497
498 dgdp_sg[i].resize(num_p);
499 for (int j=0; j<num_p; j++) {
500 if (!outArgs.supports(OUT_ARG_DgDp_sg, ii, j).none()) {
501 dgdp_sg[i][j] = outArgs.get_DgDp_sg(ii,j);
502 if (dgdp_sg[i][j].getMultiVector() != Teuchos::null)
503 dgdp_sg[i][j].getMultiVector()->init(0.0);
504 else if (dgdp_sg[i][j].getLinearOp() != Teuchos::null)
505 dgdp_sg[i][j].getLinearOp()->init(0.0);
506 }
507 }
508 }
509
510 if (do_quad) {
511 // Get quadrature data
512 const Teuchos::Array<double>& quad_weights =
513 quad->getQuadWeights();
514 const Teuchos::Array< Teuchos::Array<double> > & quad_values =
515 quad->getBasisAtQuadPoints();
516 const Teuchos::Array<double>& basis_norms = basis->norm_squared();
517
518 // Evaluate inputs at quadrature points
519 int nqp = mp_block_map->NumMyElements();
520 for (int qp=0; qp<nqp; qp++) {
521#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
522 TEUCHOS_FUNC_TIME_MONITOR_DIFF(
523 "Stokhos: SGQuadMPModelEvaluator -- Polynomial Evaluation",
524 PolyEvaluation);
525#endif
526
527 int gqp = mp_block_map->GID(qp);
528
529 if (x_sg != Teuchos::null) {
530 x_sg->evaluate(quad_values[gqp], (*x_mp)[qp]);
531 me_inargs.set_x_mp(x_mp);
532 }
533 if (x_dot_sg != Teuchos::null) {
534 x_dot_sg->evaluate(quad_values[gqp], (*x_dot_mp)[qp]);
535 me_inargs.set_x_dot_mp(x_mp);
536 }
537 for (int i=0; i<num_p_mp; i++) {
538 if (p_sg[i] != Teuchos::null) {
539 p_sg[i]->evaluate(quad_values[gqp], (*(p_mp[i]))[qp]);
540 me_inargs.set_p_mp(mp_p_index_map[i], p_mp[i]);
541 }
542 }
543
544 }
545
546 // Set OutArgs
547 if (f_sg != Teuchos::null)
548 me_outargs.set_f_mp(f_mp);
549 if (W_sg != Teuchos::null)
550 me_outargs.set_W_mp(W_mp);
551 for (int i=0; i<num_p_mp; i++) {
552 if (!dfdp_sg[i].isEmpty())
553 me_outargs.set_DfDp_mp(i, dfdp_mp[i]);
554 }
555 for (int i=0; i<num_g_mp; i++) {
556 int ii = mp_g_index_map[i];
557 if (g_sg[i] != Teuchos::null)
558 me_outargs.set_g_mp(ii, g_mp[i]);
559 if (!dgdx_dot_sg[i].isEmpty())
560 me_outargs.set_DgDx_dot_mp(ii, dgdx_dot_mp[i]);
561 if (!dgdx_sg[i].isEmpty())
562 me_outargs.set_DgDx_mp(ii, dgdx_mp[i]);
563 for (int j=0; j<num_p; j++)
564 if (!dgdp_sg[i][j].isEmpty())
565 me_outargs.set_DgDp_mp(ii, j, dgdp_mp[i][j]);
566 }
567
568
569 {
570#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
571 TEUCHOS_FUNC_TIME_MONITOR("Stokhos: SGQuadMPModelEvaluator -- Model Evaluation");
572#endif
573
574 // Evaluate multi-point model at quadrature points
575 me->evalModel(me_inargs, me_outargs);
576
577 }
578
579 // Perform integrations
580 for (int qp=0; qp<nqp; qp++) {
581#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
582 TEUCHOS_FUNC_TIME_MONITOR_DIFF(
583 "Stokhos: SGQuadMPModelEvaluator -- Polynomial Integration", Integration);
584#endif
585
586 int gqp = mp_block_map->GID(qp);
587
588 // Sum in results
589 if (f_sg != Teuchos::null) {
590 f_sg->sumIntoAllTerms(quad_weights[gqp], quad_values[gqp], basis_norms,
591 (*f_mp)[qp]);
592 }
593 if (W_sg != Teuchos::null) {
594 W_sg->sumIntoAllTerms(quad_weights[gqp], quad_values[gqp], basis_norms,
595 (*W_mp)[qp]);
596 }
597 for (int j=0; j<num_p; j++) {
598 if (!dfdp_sg[j].isEmpty()) {
599 if (dfdp_sg[j].getMultiVector() != Teuchos::null) {
600 dfdp_sg[j].getMultiVector()->sumIntoAllTerms(
601 quad_weights[gqp], quad_values[gqp], basis_norms,
602 (*(dfdp_mp[j].getMultiVector()))[qp]);
603 }
604 else if (dfdp_sg[j].getLinearOp() != Teuchos::null) {
605 dfdp_sg[j].getLinearOp()->sumIntoAllTerms(
606 quad_weights[gqp], quad_values[gqp], basis_norms,
607 (*(dfdp_mp[j].getLinearOp()))[qp]);
608 }
609 }
610 }
611 for (int i=0; i<num_g_mp; i++) {
612 if (g_sg[i] != Teuchos::null) {
613 g_sg[i]->sumIntoAllTerms(quad_weights[gqp], quad_values[gqp],
614 basis_norms, (*g_mp[i])[qp]);
615 }
616 if (!dgdx_dot_sg[i].isEmpty()) {
617 if (dgdx_dot_sg[i].getMultiVector() != Teuchos::null) {
618 dgdx_dot_sg[i].getMultiVector()->sumIntoAllTerms(
619 quad_weights[gqp], quad_values[gqp], basis_norms,
620 (*(dgdx_dot_mp[i].getMultiVector()))[qp]);
621 }
622 else if (dgdx_dot_sg[i].getLinearOp() != Teuchos::null) {
623 dgdx_dot_sg[i].getLinearOp()->sumIntoAllTerms(
624 quad_weights[gqp], quad_values[gqp], basis_norms,
625 (*(dgdx_dot_mp[i].getLinearOp()))[qp]);
626 }
627 }
628 if (!dgdx_sg[i].isEmpty()) {
629 if (dgdx_sg[i].getMultiVector() != Teuchos::null) {
630 dgdx_sg[i].getMultiVector()->sumIntoAllTerms(
631 quad_weights[gqp], quad_values[gqp], basis_norms,
632 (*(dgdx_mp[i].getMultiVector()))[qp]);
633 }
634 else if (dgdx_sg[i].getLinearOp() != Teuchos::null) {
635 dgdx_sg[i].getLinearOp()->sumIntoAllTerms(
636 quad_weights[gqp], quad_values[gqp], basis_norms,
637 (*(dgdx_mp[i].getLinearOp()))[qp]);
638 }
639 }
640 for (int j=0; j<num_p; j++) {
641 if (!dgdp_sg[i][j].isEmpty()) {
642 if (dgdp_sg[i][j].getMultiVector() != Teuchos::null) {
643 dgdp_sg[i][j].getMultiVector()->sumIntoAllTerms(
644 quad_weights[gqp], quad_values[gqp], basis_norms,
645 (*(dgdp_mp[i][j].getMultiVector()))[qp]);
646 }
647 else if (dgdp_sg[i][j].getLinearOp() != Teuchos::null) {
648 dgdp_sg[i][j].getLinearOp()->sumIntoAllTerms(
649 quad_weights[gqp], quad_values[gqp], basis_norms,
650 (*(dgdp_mp[i][j].getLinearOp()))[qp]);
651 }
652 }
653 }
654 }
655
656 }
657
658 // Now communicate partial sums across processors
659 if (mp_block_map->DistributedGlobal()) {
660 for (int i=0; i<num_g_mp; i++) {
661 if (g_sg[i] != Teuchos::null) {
662 g_sg[i]->sumAll();
663
664 // Need to do the same for all of the other out args --
665 // function needs to be added to multi-vectors and operators
666 }
667 }
668 }
669 }
670 else {
671 // Compute the non-SG functions
672 me->evalModel(me_inargs, me_outargs);
673 }
674}
A container class storing products of Epetra_MultiVector's.
A container class for products of Epetra_Vector's.
A container class for products of Epetra_Vector's.
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int l) const
Return array of parameter names.
SGQuadMPModelEvaluator(const Teuchos::RCP< EpetraExt::ModelEvaluator > &me, const Teuchos::RCP< const EpetraExt::MultiComm > &mp_comm, const Teuchos::RCP< const Epetra_Map > &mp_block_map)
Teuchos::Array< mp_vector_t > p_mp
Parameter vectors.
Teuchos::RCP< const EpetraExt::MultiComm > mp_comm
Parallel MP communicator.
void evalModel(const InArgs &inArgs, const OutArgs &outArgs) const
Evaluate model on InArgs.
Teuchos::Array< mp_vector_t > g_mp
Response vectors.
int num_p_mp
Number of multipoint parameter vectors.
Teuchos::RCP< const Epetra_Map > get_g_map(int l) const
Return observation vector map.
Teuchos::Array< MPDerivative > dgdx_mp
Response derivative.
Teuchos::Array< int > mp_g_index_map
Index map between block-g and g_mp maps.
Teuchos::Array< Teuchos::Array< MPDerivative > > dgdp_mp
Response sensitivities.
Teuchos::RCP< const Epetra_Map > get_p_map(int l) const
Return parameter vector map.
mp_vector_t x_dot_mp
Time derivative vector.
Teuchos::RCP< const Epetra_Vector > get_x_init() const
Return initial solution.
Teuchos::RCP< const Epetra_Map > mp_block_map
Map for layout of parallel MP blocks.
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.
OutArgs createOutArgs() const
Create OutArgs.
Teuchos::Array< int > mp_p_index_map
Index map between block-p and p_mp maps.
Teuchos::Array< MPDerivative > dfdp_mp
Residual derivatives.
Teuchos::RCP< const Epetra_Map > get_f_map() const
Return residual vector map.
int num_g_mp
Number of multipoint response vectors.
Teuchos::Array< MPDerivative > dgdx_dot_mp
Response derivative.
Teuchos::RCP< EpetraExt::ModelEvaluator > me
Underlying model evaluator.
Teuchos::RCP< Epetra_Operator > create_W() const
Create W = alpha*M + beta*J matrix.