Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_SGModelEvaluator_Interlaced.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
43
44#include <algorithm>
45#include "Teuchos_Assert.hpp"
46#include "EpetraExt_BlockUtility.h"
47#include "EpetraExt_BlockMultiVector.h"
52#include "Epetra_LocalMap.h"
53#include "Epetra_Export.h"
54#include "Epetra_Import.h"
55
57 const Teuchos::RCP<EpetraExt::ModelEvaluator>& me_,
58 const Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> >& sg_basis_,
59 const Teuchos::RCP<const Stokhos::Quadrature<int,double> >& sg_quad_,
60 const Teuchos::RCP<Stokhos::OrthogPolyExpansion<int,double> >& sg_exp_,
61 const Teuchos::RCP<const Stokhos::ParallelData>& sg_parallel_data_,
62 const Teuchos::RCP<Teuchos::ParameterList>& params_,
63 bool scaleOP_)
64 : me(me_),
65 sg_basis(sg_basis_),
66 sg_quad(sg_quad_),
67 sg_exp(sg_exp_),
68 params(params_),
69 num_sg_blocks(sg_basis->size()),
70 num_W_blocks(sg_basis->size()),
71 num_p_blocks(sg_basis->size()),
72 supports_x(false),
73 x_map(me->get_x_map()),
74 f_map(me->get_f_map()),
75 sg_parallel_data(sg_parallel_data_),
76 sg_comm(sg_parallel_data->getMultiComm()),
77 epetraCijk(sg_parallel_data->getEpetraCijk()),
78 serialCijk(),
79 num_p(0),
80 num_p_sg(0),
81 sg_p_index_map(),
82 sg_p_map(),
83 sg_p_names(),
84 num_g(0),
85 num_g_sg(0),
86 sg_g_index_map(),
87 sg_g_map(),
88 x_dot_sg_blocks(),
89 x_sg_blocks(),
90 f_sg_blocks(),
91 W_sg_blocks(),
92 dgdx_dot_sg_blocks(),
93 dgdx_sg_blocks(),
94 sg_x_init(),
95 sg_p_init(),
96 eval_W_with_f(false),
97 scaleOP(scaleOP_)
98{
99 if (x_map != Teuchos::null)
100 supports_x = true;
101
103 Teuchos::rcp(new Epetra_LocalMap(
104 static_cast<int>(num_sg_blocks), 0, *(sg_parallel_data->getStochasticComm())));
106 if (epetraCijk != Teuchos::null)
107 stoch_row_map = epetraCijk->getStochasticRowMap();
108
109 if (supports_x) {
110
111 // Create interlace SG x and f maps
114
117
118 // Create importer/exporter from/to overlapped distribution
120 Teuchos::rcp(new Epetra_Import(*interlace_overlapped_x_map, *get_x_map()));
123
124 // now we create the underlying Epetra block vectors
125 // that will be used by the model evaluator to construct
126 // the solution of the deterministic problem.
128
129 // Create vector blocks
133
134 // Create default sg_x_init
136 if (sg_x_init->myGID(0))
137 (*sg_x_init)[0] = *(me->get_x_init());
138
139 // Preconditioner needs an x: This is interlaced
140 my_x = Teuchos::rcp(new Epetra_Vector(*get_x_map()));
141
142
143 // setup storage for W, these are blocked in Stokhos
144 // format
146
147 // Determine W expansion type
148 std::string W_expansion_type =
149 params->get("Jacobian Expansion Type", "Full");
150 if (W_expansion_type == "Linear")
151 num_W_blocks = sg_basis->dimension() + 1;
152 else
154
155 Teuchos::RCP<Epetra_BlockMap> W_overlap_map =
156 Teuchos::rcp(new Epetra_LocalMap(
157 static_cast<int>(num_W_blocks), 0,
158 *(sg_parallel_data->getStochasticComm())));
160 Teuchos::rcp(new Stokhos::EpetraOperatorOrthogPoly(
161 sg_basis, W_overlap_map, x_map, f_map, interlace_f_map,
162 sg_comm));
163 for (unsigned int i=0; i<num_W_blocks; i++)
164 W_sg_blocks->setCoeffPtr(i, me->create_W()); // allocate a bunch of matrices
165
166 eval_W_with_f = params->get("Evaluate W with F", false);
167 }
168
169 // Parameters -- The idea here is to add new parameter vectors
170 // for the stochastic Galerkin components of the parameters
171
172 InArgs me_inargs = me->createInArgs();
173 OutArgs me_outargs = me->createOutArgs();
174 num_p = me_inargs.Np();
175
176 // Get the p_sg's supported and build index map
177 for (int i=0; i<num_p; i++) {
178 if (me_inargs.supports(IN_ARG_p_sg, i))
179 sg_p_index_map.push_back(i);
180 }
181 num_p_sg = sg_p_index_map.size();
182
183 sg_p_map.resize(num_p_sg);
184 sg_p_names.resize(num_p_sg);
185 sg_p_init.resize(num_p_sg);
186
187 // Determine parameter expansion type
188 std::string p_expansion_type =
189 params->get("Parameter Expansion Type", "Full");
190 if (p_expansion_type == "Linear")
191 num_p_blocks = sg_basis->dimension() + 1;
192 else
194
195 // Create parameter maps, names, and initial values
197 Teuchos::rcp(new Epetra_LocalMap(
198 static_cast<int>(num_p_blocks), 0,
199 *(sg_parallel_data->getStochasticComm())));
200 for (int i=0; i<num_p_sg; i++) {
201 Teuchos::RCP<const Epetra_Map> p_map = me->get_p_map(sg_p_index_map[i]);
202 sg_p_map[i] =
203 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(
204 *p_map, *overlapped_stoch_p_map, *sg_comm));
205
206 Teuchos::RCP<const Teuchos::Array<std::string> > p_names =
207 me->get_p_names(sg_p_index_map[i]);
208 if (p_names != Teuchos::null) {
209 sg_p_names[i] =
210 Teuchos::rcp(new Teuchos::Array<std::string>(num_sg_blocks*(p_names->size())));
211 for (int j=0; j<p_names->size(); j++) {
212 std::stringstream ss;
213 ss << (*p_names)[j] << " -- SG Coefficient " << i;
214 (*sg_p_names[i])[j] = ss.str();
215 }
216 }
217
218 // Create default sg_p_init
220 sg_p_init[i]->init(0.0);
221 }
222
223 // Responses -- The idea here is to add new parameter vectors
224 // for the stochastic Galerkin components of the respones
225
226 // Get number of SG parameters model supports derivatives w.r.t.
227 num_g = me_outargs.Ng();
228
229 // Get the g_sg's supported and build index map
230 for (int i=0; i<num_g; i++) {
231 if (me_outargs.supports(OUT_ARG_g_sg, i))
232 sg_g_index_map.push_back(i);
233 }
234 num_g_sg = sg_g_index_map.size();
235
236 sg_g_map.resize(num_g_sg);
238 dgdx_sg_blocks.resize(num_g_sg);
239
240 // Create response maps
241 for (int i=0; i<num_g_sg; i++) {
242 Teuchos::RCP<const Epetra_Map> g_map = me->get_g_map(sg_g_index_map[i]);
243 sg_g_map[i] =
244 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(
246
247 // Create dg/dxdot, dg/dx SG blocks
248 if (supports_x) {
252 dgdx_sg_blocks[i] =
255 }
256 }
257
258 // We don't support parallel for dgdx yet, so build a new EpetraCijk
259 if (supports_x) {
260 serialCijk =
261 Teuchos::rcp(new Stokhos::EpetraSparse3Tensor(sg_basis,
262 epetraCijk->getCijk(),
263 sg_comm,
265 }
266
267}
268
269// Overridden from EpetraExt::ModelEvaluator
270
271Teuchos::RCP<const Epetra_Map>
273{
274 return interlace_x_map;
275}
276
277Teuchos::RCP<const Epetra_Map>
279{
280 return interlace_f_map;
281}
282
283Teuchos::RCP<const Epetra_Map>
285{
286 TEUCHOS_TEST_FOR_EXCEPTION(l < 0 || l >= num_p + num_p_sg, std::logic_error,
287 "Error! Invalid p map index " << l);
288 if (l < num_p)
289 return me->get_p_map(l);
290 else
291 return sg_p_map[l-num_p];
292
293 return Teuchos::null;
294}
295
296Teuchos::RCP<const Epetra_Map>
298{
299 TEUCHOS_TEST_FOR_EXCEPTION(l < 0 || l >= num_g_sg, std::logic_error,
300 "Error! Invalid g map index " << l);
301 return sg_g_map[l];
302}
303
304Teuchos::RCP<const Teuchos::Array<std::string> >
306{
307 TEUCHOS_TEST_FOR_EXCEPTION(l < 0 || l >= num_p + num_p_sg, std::logic_error,
308 "Error! Invalid p map index " << l);
309 if (l < num_p)
310 return me->get_p_names(l);
311 else
312 return sg_p_names[l-num_p];
313
314 return Teuchos::null;
315}
316
317Teuchos::RCP<const Epetra_Vector>
319{
320 return sg_x_init->getBlockVector();
321}
322
323Teuchos::RCP<const Epetra_Vector>
325{
326 TEUCHOS_TEST_FOR_EXCEPTION(l < 0 || l >= num_p + num_p_sg, std::logic_error,
327 "Error! Invalid p map index " << l);
328 if (l < num_p)
329 return me->get_p_init(l);
330 else
331 return sg_p_init[l-num_p]->getBlockVector();
332
333 return Teuchos::null;
334}
335
336Teuchos::RCP<Epetra_Operator>
338{
339 if (supports_x) {
340 Teuchos::RCP<Teuchos::ParameterList> sgOpParams =
341 Teuchos::rcp(&(params->sublist("SG Operator")), false);
342 Teuchos::RCP<Epetra_CrsMatrix> W_crs;
343 W_crs = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(me->create_W(), true);
344 Teuchos::RCP<const Epetra_CrsGraph> W_graph =
345 Teuchos::rcp(&(W_crs->Graph()), false);
346
347 my_W = Teuchos::rcp(new Stokhos::InterlacedOperator(sg_comm, sg_basis, epetraCijk, W_graph,sgOpParams));
348 my_W->setupOperator(W_sg_blocks);
349
350 return my_W;
351 }
352
353 return Teuchos::null;
354}
355
356EpetraExt::ModelEvaluator::InArgs
358{
359 InArgsSetup inArgs;
360 InArgs me_inargs = me->createInArgs();
361
362 inArgs.setModelEvalDescription(this->description());
363 inArgs.set_Np(num_p + num_p_sg);
364 inArgs.setSupports(IN_ARG_x_dot, me_inargs.supports(IN_ARG_x_dot_sg));
365 inArgs.setSupports(IN_ARG_x, me_inargs.supports(IN_ARG_x_sg));
366 inArgs.setSupports(IN_ARG_t, me_inargs.supports(IN_ARG_t));
367 inArgs.setSupports(IN_ARG_alpha, me_inargs.supports(IN_ARG_alpha));
368 inArgs.setSupports(IN_ARG_beta, me_inargs.supports(IN_ARG_beta));
369 inArgs.setSupports(IN_ARG_sg_basis, me_inargs.supports(IN_ARG_sg_basis));
370 inArgs.setSupports(IN_ARG_sg_quadrature,
371 me_inargs.supports(IN_ARG_sg_quadrature));
372 inArgs.setSupports(IN_ARG_sg_expansion,
373 me_inargs.supports(IN_ARG_sg_expansion));
374
375 return inArgs;
376}
377
378EpetraExt::ModelEvaluator::OutArgs
380{
381 OutArgsSetup outArgs;
382 OutArgs me_outargs = me->createOutArgs();
383
384 outArgs.setModelEvalDescription(this->description());
385 outArgs.set_Np_Ng(num_p+num_p_sg, num_g_sg);
386 outArgs.setSupports(OUT_ARG_f, me_outargs.supports(OUT_ARG_f_sg));
387 outArgs.setSupports(OUT_ARG_W, me_outargs.supports(OUT_ARG_W_sg));
388 outArgs.setSupports(OUT_ARG_WPrec, false);
389 for (int j=0; j<num_p; j++)
390 outArgs.setSupports(OUT_ARG_DfDp, j,
391 me_outargs.supports(OUT_ARG_DfDp_sg, j));
392 for (int i=0; i<num_g_sg; i++) {
393 int ii = sg_g_index_map[i];
394// if (!me_outargs.supports(OUT_ARG_DgDx_dot_sg, ii).none())
395// outArgs.setSupports(OUT_ARG_DgDx_dot, i, DERIV_LINEAR_OP);
396// if (!me_outargs.supports(OUT_ARG_DgDx_sg, i).none())
397// outArgs.setSupports(OUT_ARG_DgDx, i, DERIV_LINEAR_OP);
398 for (int j=0; j<num_p; j++)
399 outArgs.setSupports(OUT_ARG_DgDp, i, j,
400 me_outargs.supports(OUT_ARG_DgDp_sg, ii, j));
401 }
402
403 // We do not support derivatives w.r.t. the new SG parameters, so their
404 // support defaults to none.
405
406 return outArgs;
407}
408
409void
411 const OutArgs& outArgs) const
412{
413 // Get the input arguments
414 Teuchos::RCP<const Epetra_Vector> x;
415 if (inArgs.supports(IN_ARG_x)) {
416 x = inArgs.get_x();
417 if (x != Teuchos::null)
418 *my_x = *x;
419 }
420 Teuchos::RCP<const Epetra_Vector> x_dot;
421 if (inArgs.supports(IN_ARG_x_dot))
422 x_dot = inArgs.get_x_dot();
423
424 // Get the output arguments
425 EpetraExt::ModelEvaluator::Evaluation<Epetra_Vector> f_out;
426 if (outArgs.supports(OUT_ARG_f))
427 f_out = outArgs.get_f();
428 Teuchos::RCP<Epetra_Operator> W_out;
429 if (outArgs.supports(OUT_ARG_W))
430 W_out = outArgs.get_W();
431
432 // Create underlying inargs
433 InArgs me_inargs = me->createInArgs();
434 if (x != Teuchos::null) {
435 Teuchos::RCP<Epetra_Vector> overlapped_x
436 = Teuchos::rcp(new Epetra_Vector(*interlace_overlapped_x_map));
437 overlapped_x->Import(*x,*interlace_overlapped_x_importer,Insert);
438
439 // x_sg_blocks->getBlockVector()->Import(*x, *interlace_overlapped_x_importer,
440 // Insert);
441
442 copyToPolyOrthogVector(*overlapped_x,*x_sg_blocks);
443 me_inargs.set_x_sg(x_sg_blocks);
444 }
445 if (x_dot != Teuchos::null) {
446 Teuchos::RCP<Epetra_Vector> overlapped_x_dot
447 = Teuchos::rcp(new Epetra_Vector(*interlace_overlapped_x_map));
448 overlapped_x_dot->Import(*x_dot,*interlace_overlapped_x_importer,Insert);
449
450 // x_dot_sg_blocks->getBlockVector()->Import(*x_dot, *interlace_overlapped_x_importer,
451 // Insert);
452
453 copyToPolyOrthogVector(*overlapped_x_dot,*x_dot_sg_blocks);
454 me_inargs.set_x_dot_sg(x_dot_sg_blocks);
455 }
456 if (me_inargs.supports(IN_ARG_alpha))
457 me_inargs.set_alpha(inArgs.get_alpha());
458 if (me_inargs.supports(IN_ARG_beta))
459 me_inargs.set_beta(inArgs.get_beta());
460 if (me_inargs.supports(IN_ARG_t))
461 me_inargs.set_t(inArgs.get_t());
462 if (me_inargs.supports(IN_ARG_sg_basis)) {
463 if (inArgs.get_sg_basis() != Teuchos::null)
464 me_inargs.set_sg_basis(inArgs.get_sg_basis());
465 else
466 me_inargs.set_sg_basis(sg_basis);
467 }
468 if (me_inargs.supports(IN_ARG_sg_quadrature)) {
469 if (inArgs.get_sg_quadrature() != Teuchos::null)
470 me_inargs.set_sg_quadrature(inArgs.get_sg_quadrature());
471 else
472 me_inargs.set_sg_quadrature(sg_quad);
473 }
474 if (me_inargs.supports(IN_ARG_sg_expansion)) {
475 if (inArgs.get_sg_expansion() != Teuchos::null)
476 me_inargs.set_sg_expansion(inArgs.get_sg_expansion());
477 else
478 me_inargs.set_sg_expansion(sg_exp);
479 }
480
481 // Pass parameters
482 for (int i=0; i<num_p; i++)
483 me_inargs.set_p(i, inArgs.get_p(i));
484 for (int i=0; i<num_p_sg; i++) {
485 Teuchos::RCP<const Epetra_Vector> p = inArgs.get_p(i+num_p);
486
487 // We always need to pass in the SG parameters, so just use
488 // the initial parameters if it is NULL
489 if (p == Teuchos::null)
490 p = sg_p_init[i]->getBlockVector();
491
492 // Convert block p to SG polynomial
493 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> p_sg =
494 create_p_sg(sg_p_index_map[i], View, p.get());
495 me_inargs.set_p_sg(sg_p_index_map[i], p_sg);
496 }
497
498 // Create underlying outargs
499 OutArgs me_outargs = me->createOutArgs();
500
501 // f
502 if (f_out != Teuchos::null) {
503 me_outargs.set_f_sg(f_sg_blocks);
504 if (eval_W_with_f)
505 me_outargs.set_W_sg(W_sg_blocks);
506 }
507
508 // W
509 if (W_out != Teuchos::null && !eval_W_with_f )
510 me_outargs.set_W_sg(W_sg_blocks);
511
512 // df/dp -- deterministic p
513 for (int i=0; i<outArgs.Np(); i++) {
514 if (!outArgs.supports(OUT_ARG_DfDp, i).none()) {
515 Derivative dfdp = outArgs.get_DfDp(i);
516 if (dfdp.getMultiVector() != Teuchos::null) {
517 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> dfdp_sg;
518 if (dfdp.getMultiVectorOrientation() == DERIV_MV_BY_COL)
519 dfdp_sg =
521 sg_basis, overlapped_stoch_row_map,
522 me->get_f_map(), interlace_overlapped_f_map, sg_comm,
523 me->get_p_map(i)->NumMyElements()));
524 else if (dfdp.getMultiVectorOrientation() == DERIV_TRANS_MV_BY_ROW)
525 dfdp_sg =
527 sg_basis, overlapped_stoch_row_map,
528 me->get_p_map(i), sg_comm,
529 me->get_f_map()->NumMyElements()));
530 me_outargs.set_DfDp_sg(i,
531 SGDerivative(dfdp_sg,
532 dfdp.getMultiVectorOrientation()));
533 }
534 TEUCHOS_TEST_FOR_EXCEPTION(dfdp.getLinearOp() != Teuchos::null, std::logic_error,
535 "Error! Stokhos::SGModelEvaluator_Interlaced::evalModel " <<
536 "cannot handle operator form of df/dp!");
537 }
538 }
539
540 // Responses (g, dg/dx, dg/dp, ...)
541 for (int i=0; i<num_g_sg; i++) {
542 int ii = sg_g_index_map[i];
543
544 // g
545 Teuchos::RCP<Epetra_Vector> g = outArgs.get_g(i);
546 if (g != Teuchos::null) {
547 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> g_sg =
548 create_g_sg(sg_g_index_map[i], View, g.get());
549 me_outargs.set_g_sg(i, g_sg);
550 }
551
552 // dg/dxdot
553 if (outArgs.supports(OUT_ARG_DgDx_dot, i).supports(DERIV_LINEAR_OP)) {
554 Derivative dgdx_dot = outArgs.get_DgDx_dot(i);
555 if (dgdx_dot.getLinearOp() != Teuchos::null) {
556 Teuchos::RCP<Stokhos::SGOperator> op =
557 Teuchos::rcp_dynamic_cast<Stokhos::SGOperator>(
558 dgdx_dot.getLinearOp(), true);
559 Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > sg_blocks =
560 op->getSGPolynomial();
561 if (me_outargs.supports(OUT_ARG_DgDx, ii).supports(DERIV_LINEAR_OP))
562 me_outargs.set_DgDx_dot_sg(ii, sg_blocks);
563 else {
564 for (unsigned int k=0; k<num_sg_blocks; k++) {
565 Teuchos::RCP<Epetra_MultiVector> mv =
566 Teuchos::rcp_dynamic_cast<Stokhos::EpetraMultiVectorOperator>(
567 sg_blocks->getCoeffPtr(k), true)->getMultiVector();
568 dgdx_dot_sg_blocks[i]->setCoeffPtr(k, mv);
569 }
570 if (me_outargs.supports(OUT_ARG_DgDx_dot_sg, ii).supports(DERIV_MV_BY_COL))
571 me_outargs.set_DgDx_dot_sg(ii, SGDerivative(dgdx_dot_sg_blocks[i],
572 DERIV_MV_BY_COL));
573 else
574 me_outargs.set_DgDx_dot_sg(ii, SGDerivative(dgdx_dot_sg_blocks[i],
575 DERIV_TRANS_MV_BY_ROW));
576 }
577 }
578 TEUCHOS_TEST_FOR_EXCEPTION(dgdx_dot.getLinearOp() == Teuchos::null &&
579 dgdx_dot.isEmpty() == false,
580 std::logic_error,
581 "Error! Stokhos::SGModelEvaluator_Interlaced::evalModel: " <<
582 "Operator form of dg/dxdot is required!");
583 }
584
585 // dg/dx
586 if (outArgs.supports(OUT_ARG_DgDx, i).supports(DERIV_LINEAR_OP)) {
587 Derivative dgdx = outArgs.get_DgDx(i);
588 if (dgdx.getLinearOp() != Teuchos::null) {
589 Teuchos::RCP<Stokhos::SGOperator> op =
590 Teuchos::rcp_dynamic_cast<Stokhos::SGOperator>(
591 dgdx.getLinearOp(), true);
592 Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > sg_blocks =
593 op->getSGPolynomial();
594 if (me_outargs.supports(OUT_ARG_DgDx, ii).supports(DERIV_LINEAR_OP))
595 me_outargs.set_DgDx_sg(i, sg_blocks);
596 else {
597 for (unsigned int k=0; k<num_sg_blocks; k++) {
598 Teuchos::RCP<Epetra_MultiVector> mv =
599 Teuchos::rcp_dynamic_cast<Stokhos::EpetraMultiVectorOperator>(
600 sg_blocks->getCoeffPtr(k), true)->getMultiVector();
601 dgdx_sg_blocks[i]->setCoeffPtr(k, mv);
602 }
603 if (me_outargs.supports(OUT_ARG_DgDx_sg, ii).supports(DERIV_MV_BY_COL))
604 me_outargs.set_DgDx_sg(ii, SGDerivative(dgdx_sg_blocks[i],
605 DERIV_MV_BY_COL));
606 else
607 me_outargs.set_DgDx_sg(ii, SGDerivative(dgdx_sg_blocks[i],
608 DERIV_TRANS_MV_BY_ROW));
609 }
610 }
611 TEUCHOS_TEST_FOR_EXCEPTION(dgdx.getLinearOp() == Teuchos::null &&
612 dgdx.isEmpty() == false,
613 std::logic_error,
614 "Error! Stokhos::SGModelEvaluator_Interlaced::evalModel: " <<
615 "Operator form of dg/dxdot is required!");
616 }
617
618 // dg/dp -- deterministic p
619 // Rembember, no derivatives w.r.t. sg parameters
620 for (int j=0; j<num_p; j++) {
621 if (!outArgs.supports(OUT_ARG_DgDp, i, j).none()) {
622 Derivative dgdp = outArgs.get_DgDp(i,j);
623 if (dgdp.getMultiVector() != Teuchos::null) {
624 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> dgdp_sg;
625 if (dgdp.getMultiVectorOrientation() == DERIV_MV_BY_COL)
626 dgdp_sg =
628 sg_basis, overlapped_stoch_row_map,
629 me->get_g_map(ii), sg_g_map[i], sg_comm,
630 View, *(dgdp.getMultiVector())));
631 else if (dgdp.getMultiVectorOrientation() == DERIV_TRANS_MV_BY_ROW) {
632 Teuchos::RCP<const Epetra_BlockMap> product_map =
633 Teuchos::rcp(&(dgdp.getMultiVector()->Map()),false);
634 dgdp_sg =
636 sg_basis, overlapped_stoch_row_map,
637 me->get_p_map(j), product_map, sg_comm,
638 View, *(dgdp.getMultiVector())));
639 }
640 me_outargs.set_DgDp_sg(ii, j,
641 SGDerivative(dgdp_sg,
642 dgdp.getMultiVectorOrientation()));
643 }
644 TEUCHOS_TEST_FOR_EXCEPTION(dgdp.getLinearOp() != Teuchos::null,
645 std::logic_error,
646 "Error! Stokhos::SGModelEvaluator_Interlaced::evalModel " <<
647 "cannot handle operator form of dg/dp!");
648 }
649 }
650
651 }
652
653 // Compute the functions
654 me->evalModel(me_inargs, me_outargs);
655
656 // Copy block SG components for W
657 if ((W_out != Teuchos::null || (eval_W_with_f && f_out != Teuchos::null)) ) {
658
659 Teuchos::RCP<Epetra_Operator> W;
660 if (W_out != Teuchos::null)
661 W = W_out;
662 else
663 W = my_W;
664 Teuchos::RCP<Stokhos::SGOperator> W_sg =
665 Teuchos::rcp_dynamic_cast<Stokhos::SGOperator>(W, true);
666 W_sg->setupOperator(W_sg_blocks);
667 }
668
669 // f
670 if (f_out!=Teuchos::null){
671 if (!scaleOP)
672 for (int i=0; i<sg_basis->size(); i++)
673 (*f_sg_blocks)[i].Scale(sg_basis->norm_squared(i));
674
675 Teuchos::RCP<Epetra_Vector> overlapped_f
676 = Teuchos::rcp(new Epetra_Vector(*interlace_overlapped_f_map));
677 copyToInterlacedVector(*f_sg_blocks,*overlapped_f);
678 f_out->Export(*overlapped_f,*interlace_overlapped_f_exporter,Insert);
679 }
680
681 // df/dp -- determinsistic p
682 for (int i=0; i<num_p; i++) {
683 if (!outArgs.supports(OUT_ARG_DfDp, i).none()) {
684 Derivative dfdp = outArgs.get_DfDp(i);
685 SGDerivative dfdp_sg = me_outargs.get_DfDp_sg(i);
686 if (dfdp.getMultiVector() != Teuchos::null) {
687 dfdp.getMultiVector()->Export(
688 *(dfdp_sg.getMultiVector()->getBlockMultiVector()),
689 *interlace_overlapped_f_exporter, Insert);
690 }
691 }
692 }
693}
694
695void
697 const Stokhos::EpetraVectorOrthogPoly& x_sg_in)
698{
699 *sg_x_init = x_sg_in;
700}
701
702Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly>
704{
705 return sg_x_init;
706}
707
708void
710 int i, const Stokhos::EpetraVectorOrthogPoly& p_sg_in)
711{
712 *sg_p_init[i] = p_sg_in;
713}
714
715Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly>
717{
718 return sg_p_init[l];
719}
720
721Teuchos::Array<int>
723{
724 return sg_p_index_map;
725}
726
727Teuchos::Array<int>
729{
730 return sg_g_index_map;
731}
732
733Teuchos::Array< Teuchos::RCP<const Epetra_Map> >
735{
736 Teuchos::Array< Teuchos::RCP<const Epetra_Map> > base_maps(num_g);
737 for (int i=0; i<num_g; i++)
738 base_maps[i] = me->get_g_map(i);
739 return base_maps;
740 }
741
742Teuchos::RCP<const Epetra_BlockMap>
744{
745 return overlapped_stoch_row_map;
746}
747
748Teuchos::RCP<const Epetra_BlockMap>
750{
751 return interlace_overlapped_x_map;
752}
753
754Teuchos::RCP<const Epetra_Import>
756{
757 return interlace_overlapped_x_importer;
758}
759
760Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
762 const Epetra_Vector* v) const
763{
764 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_x;
765 if (v == NULL)
766 sg_x = Teuchos::rcp(new Stokhos::EpetraVectorOrthogPoly(
767 sg_basis, stoch_row_map, x_map, get_x_map(), sg_comm));
768 else
769 sg_x = Teuchos::rcp(new Stokhos::EpetraVectorOrthogPoly(
770 sg_basis, stoch_row_map, x_map, get_x_map(), sg_comm,
771 CV, *v));
772 return sg_x;
773}
774
775Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
777 const Epetra_Vector* v) const
778{
779 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_x;
780 if (v == NULL)
781 sg_x = Teuchos::rcp(new Stokhos::EpetraVectorOrthogPoly(
782 sg_basis, overlapped_stoch_row_map, x_map,
783 get_x_sg_overlap_map(), sg_comm));
784 else
785 sg_x = Teuchos::rcp(new Stokhos::EpetraVectorOrthogPoly(
786 sg_basis, overlapped_stoch_row_map, x_map,
787 get_x_sg_overlap_map(), sg_comm, CV, *v));
788 return sg_x;
789}
790
791Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly>
793 const Epetra_MultiVector* v) const
794{
795 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> sg_x;
796 if (v == NULL)
797 sg_x = Teuchos::rcp(new Stokhos::EpetraMultiVectorOrthogPoly(
798 sg_basis, stoch_row_map, x_map, get_x_map(), sg_comm,
799 num_vecs));
800 else
801 sg_x = Teuchos::rcp(new Stokhos::EpetraMultiVectorOrthogPoly(
802 sg_basis, stoch_row_map, x_map, get_x_map(), sg_comm,
803 CV, *v));
804 return sg_x;
805}
806
807Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly>
809 int num_vecs,
811 const Epetra_MultiVector* v) const
812{
813 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> sg_x;
814 if (v == NULL)
815 sg_x = Teuchos::rcp(new Stokhos::EpetraMultiVectorOrthogPoly(
816 sg_basis, overlapped_stoch_row_map, x_map,
817 get_x_sg_overlap_map(), sg_comm, num_vecs));
818 else
819 sg_x = Teuchos::rcp(new Stokhos::EpetraMultiVectorOrthogPoly(
820 sg_basis, overlapped_stoch_row_map, x_map,
821 get_x_sg_overlap_map(), sg_comm, CV, *v));
822 return sg_x;
823}
824
825Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
827 const Epetra_Vector* v) const
828{
829 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_p;
830 Teuchos::Array<int>::const_iterator it = std::find(sg_p_index_map.begin(),
831 sg_p_index_map.end(),
832 l);
833 TEUCHOS_TEST_FOR_EXCEPTION(it == sg_p_index_map.end(), std::logic_error,
834 "Error! Invalid p map index " << l);
835 int ll = it - sg_p_index_map.begin();
836 if (v == NULL)
837 sg_p = Teuchos::rcp(new Stokhos::EpetraVectorOrthogPoly(
838 sg_basis, overlapped_stoch_p_map, me->get_p_map(l),
839 sg_p_map[ll], sg_comm));
840 else
841 sg_p = Teuchos::rcp(new Stokhos::EpetraVectorOrthogPoly(
842 sg_basis, overlapped_stoch_p_map, me->get_p_map(l),
843 sg_p_map[ll], sg_comm, CV, *v));
844 return sg_p;
845}
846
847Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
849 const Epetra_Vector* v) const
850{
851 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_f;
852 if (v == NULL)
853 sg_f = Teuchos::rcp(new Stokhos::EpetraVectorOrthogPoly(
854 sg_basis, stoch_row_map, f_map, interlace_f_map, sg_comm));
855 else
856 sg_f = Teuchos::rcp(new Stokhos::EpetraVectorOrthogPoly(
857 sg_basis, stoch_row_map, f_map, interlace_f_map, sg_comm,
858 CV, *v));
859 return sg_f;
860}
861
862Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
864 const Epetra_Vector* v) const
865{
866 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_f;
867 if (v == NULL)
868 sg_f = Teuchos::rcp(new Stokhos::EpetraVectorOrthogPoly(
869 sg_basis, overlapped_stoch_row_map, f_map,
870 interlace_overlapped_f_map, sg_comm));
871 else
872 sg_f = Teuchos::rcp(new Stokhos::EpetraVectorOrthogPoly(
873 sg_basis, overlapped_stoch_row_map, f_map,
874 interlace_overlapped_f_map, sg_comm, CV, *v));
875 return sg_f;
876}
877
878Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly>
880 int num_vecs,
882 const Epetra_MultiVector* v) const
883{
884 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> sg_f;
885 if (v == NULL)
886 sg_f = Teuchos::rcp(new Stokhos::EpetraMultiVectorOrthogPoly(
887 sg_basis, stoch_row_map, f_map, interlace_f_map, sg_comm,
888 num_vecs));
889 else
890 sg_f = Teuchos::rcp(new Stokhos::EpetraMultiVectorOrthogPoly(
891 sg_basis, stoch_row_map, f_map, interlace_f_map, sg_comm,
892 CV, *v));
893 return sg_f;
894}
895
896Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly>
898 int num_vecs,
900 const Epetra_MultiVector* v) const
901{
902 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> sg_f;
903 if (v == NULL)
904 sg_f = Teuchos::rcp(new Stokhos::EpetraMultiVectorOrthogPoly(
905 sg_basis, overlapped_stoch_row_map, f_map,
906 interlace_overlapped_f_map, sg_comm, num_vecs));
907 else
908 sg_f = Teuchos::rcp(new Stokhos::EpetraMultiVectorOrthogPoly(
909 sg_basis, overlapped_stoch_row_map, f_map,
910 interlace_overlapped_f_map, sg_comm, CV, *v));
911 return sg_f;
912}
913
914Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
916 const Epetra_Vector* v) const
917{
918 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_g;
919 Teuchos::Array<int>::const_iterator it = std::find(sg_g_index_map.begin(),
920 sg_g_index_map.end(),
921 l);
922 TEUCHOS_TEST_FOR_EXCEPTION(it == sg_g_index_map.end(), std::logic_error,
923 "Error! Invalid g map index " << l);
924 int ll = it - sg_g_index_map.begin();
925 if (v == NULL)
926 sg_g = Teuchos::rcp(new Stokhos::EpetraVectorOrthogPoly(
927 sg_basis, overlapped_stoch_row_map,
928 me->get_g_map(l),
929 sg_g_map[ll], sg_comm));
930 else
931 sg_g = Teuchos::rcp(new Stokhos::EpetraVectorOrthogPoly(
932 sg_basis, overlapped_stoch_row_map,
933 me->get_g_map(l),
934 sg_g_map[ll], sg_comm, CV, *v));
935 return sg_g;
936}
937
938Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly>
941 const Epetra_MultiVector* v) const
942{
943 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> sg_g;
944 Teuchos::Array<int>::const_iterator it = std::find(sg_g_index_map.begin(),
945 sg_g_index_map.end(),
946 l);
947 TEUCHOS_TEST_FOR_EXCEPTION(it == sg_g_index_map.end(), std::logic_error,
948 "Error! Invalid g map index " << l);
949 int ll = it - sg_g_index_map.begin();
950 if (v == NULL)
951 sg_g = Teuchos::rcp(new Stokhos::EpetraMultiVectorOrthogPoly(
952 sg_basis, overlapped_stoch_row_map,
953 me->get_g_map(l),
954 sg_g_map[ll], sg_comm, num_vecs));
955 else
956 sg_g = Teuchos::rcp(new Stokhos::EpetraMultiVectorOrthogPoly(
957 sg_basis, overlapped_stoch_row_map,
958 me->get_g_map(l),
959 sg_g_map[ll], sg_comm, CV, *v));
960 return sg_g;
961}
962
963Teuchos::RCP<Epetra_Map>
965{
966 int stochaUnks = stocha_map.NumMyElements();
967 int determUnks = determ_map.NumMyElements();
968
969 // these must be equal for the "adapt" model evaluator
970 TEUCHOS_ASSERT(stocha_map.NumGlobalElements()==stochaUnks);
971
972 // build interlaced unknowns
973 std::vector<int> interlacedUnks(stochaUnks*determUnks);
974 std::size_t i=0;
975 for(int d=0;d<determUnks;d++)
976 for(int s=0;s<stochaUnks;s++,i++)
977 interlacedUnks[i] = determ_map.GID(d)*stochaUnks+s;
978
979 return Teuchos::rcp(new Epetra_Map(-1,interlacedUnks.size(),&interlacedUnks[0],0,determ_map.Comm()));
980}
981
983 Epetra_Vector & x)
984{
985 std::size_t numBlocks = x_sg.size();
986 Teuchos::RCP<const EpetraExt::BlockVector> bv_x = x_sg.getBlockVector();
987
988 // loop over all blocks
989 for(std::size_t blk=0;blk<numBlocks;blk++) {
990 const Epetra_Vector & v = *bv_x->GetBlock(blk);
991
992 for(int dof=0;dof<v.MyLength();dof++)
993 x[dof*numBlocks+blk] = v[dof];
994 }
995}
996
999{
1000 std::size_t numBlocks = x_sg.size();
1001 Teuchos::RCP<EpetraExt::BlockVector> bv_x = x_sg.getBlockVector();
1002
1003 // loop over all blocks
1004 for(std::size_t blk=0;blk<numBlocks;blk++) {
1005 Epetra_Vector & v = *bv_x->GetBlock(blk);
1006
1007 for(int dof=0;dof<v.MyLength();dof++)
1008 v[dof] = x[dof*numBlocks+blk];
1009 }
1010}
Insert
Epetra_DataAccess
int GID(int LID) const
int NumGlobalElements() const
const Epetra_Comm & Comm() const
int NumMyElements() const
int MyLength() const
A container class storing an orthogonal polynomial whose coefficients are vectors,...
A container class storing an orthogonal polynomial whose coefficients are vectors,...
A container class storing an orthogonal polynomial whose coefficients are vectors,...
An Epetra operator representing the block stochastic Galerkin operator generated by fully assembling ...
Abstract base class for multivariate orthogonal polynomials.
Abstract base class for orthogonal polynomial-based expansions.
ordinal_type size() const
Return size.
Teuchos::RCP< EpetraExt::BlockVector > getBlockVector()
Get block vector.
Abstract base class for quadrature methods.
Teuchos::RCP< const Epetra_Map > get_p_map(int l) const
Return parameter vector map.
Teuchos::RCP< const Epetra_Vector > get_p_init(int l) const
Return initial parameters.
Teuchos::RCP< const Stokhos::EpetraVectorOrthogPoly > get_p_sg_init(int l) const
Return initial SG parameters.
Teuchos::RCP< Epetra_Operator > create_W() const
Create W = alpha*M + beta*J matrix.
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int l) const
Return array of parameter names.
Teuchos::RCP< Epetra_Import > interlace_overlapped_x_importer
Importer from SG to SG-overlapped maps.
Teuchos::RCP< const Epetra_BlockMap > get_x_sg_overlap_map() const
Return x sg overlap map.
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > create_x_sg_overlap(Epetra_DataAccess CV=Copy, const Epetra_Vector *v=NULL) const
Create vector orthog poly using x map and overlap sg map.
static Teuchos::RCP< Epetra_Map > buildInterlaceMap(const Epetra_BlockMap &determ_map, const Epetra_BlockMap &stocha_map)
Teuchos::RCP< const Epetra_Map > get_g_map(int l) const
Return response map.
Teuchos::RCP< Stokhos::EpetraMultiVectorOrthogPoly > create_x_mv_sg(int num_vecs, Epetra_DataAccess CV=Copy, const Epetra_MultiVector *v=NULL) const
Create vector orthog poly using x map and owned sg map.
static void copyToPolyOrthogVector(const Epetra_Vector &x, Stokhos::EpetraVectorOrthogPoly &x_sg)
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > create_g_sg(int l, Epetra_DataAccess CV=Copy, const Epetra_Vector *v=NULL) const
Create vector orthog poly using g map.
Teuchos::RCP< const Stokhos::EpetraSparse3Tensor > serialCijk
Serial Epetra Cijk for dgdx*.
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > x_dot_sg_blocks
x_dot stochastic Galerkin components
Teuchos::Array< Teuchos::RCP< const Epetra_Map > > get_g_sg_base_maps() const
Get base maps of SG responses.
Teuchos::Array< Teuchos::RCP< Teuchos::Array< std::string > > > sg_p_names
SG coefficient parameter names.
Teuchos::Array< Teuchos::RCP< Stokhos::EpetraMultiVectorOrthogPoly > > dgdx_sg_blocks
dg/dx stochastic Galerkin components
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > create_f_sg(Epetra_DataAccess CV=Copy, const Epetra_Vector *v=NULL) const
Create vector orthog poly using f map and owned sg map.
bool eval_W_with_f
Whether to always evaluate W with f.
static void copyToInterlacedVector(const Stokhos::EpetraVectorOrthogPoly &x_sg, Epetra_Vector &x)
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > sg_x_init
SG initial x.
Teuchos::Array< int > sg_p_index_map
Index map between block-p and p_sg maps.
Teuchos::RCP< const Epetra_BlockMap > overlapped_stoch_row_map
Overlapped map for stochastic blocks (local map)
unsigned int num_p_blocks
Number of p stochastic blocks (may be smaller than num_sg_blocks)
Teuchos::RCP< const Stokhos::EpetraSparse3Tensor > epetraCijk
Epetra Cijk.
Teuchos::RCP< Stokhos::EpetraMultiVectorOrthogPoly > create_f_mv_sg(int num_vecs, Epetra_DataAccess CV=Copy, const Epetra_MultiVector *v=NULL) const
Create multi-vector orthog poly using f map and owned sg map.
void evalModel(const InArgs &inArgs, const OutArgs &outArgs) const
Evaluate model on InArgs.
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > create_x_sg(Epetra_DataAccess CV=Copy, const Epetra_Vector *v=NULL) const
Create vector orthog poly using x map and owned sg map.
Teuchos::RCP< Stokhos::EpetraMultiVectorOrthogPoly > create_f_mv_sg_overlap(int num_vecs, Epetra_DataAccess CV=Copy, const Epetra_MultiVector *v=NULL) const
Create multi-vector orthog poly using f map and overlap sg map.
bool supports_x
Whether we support x (and thus f and W)
Teuchos::RCP< const Epetra_BlockMap > overlapped_stoch_p_map
Overlapped map for p stochastic blocks (local map)
int num_p_sg
Number of stochastic parameter vectors.
Teuchos::RCP< const Epetra_BlockMap > stoch_row_map
Map for stochastic blocks.
Teuchos::RCP< EpetraExt::ModelEvaluator > me
Underlying model evaluator.
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > f_sg_blocks
f stochastic Galerkin components
int num_p
Number of parameter vectors of underlying model evaluator.
Teuchos::RCP< Epetra_Vector > my_x
x pointer for evaluating preconditioner
unsigned int num_sg_blocks
Number of stochastic blocks.
Teuchos::Array< int > get_p_sg_map_indices() const
Get indices of SG parameters.
Teuchos::RCP< const Epetra_Map > interlace_f_map
Block SG residual map.
Teuchos::RCP< const Epetra_Vector > get_x_init() const
Return initial solution.
Teuchos::Array< Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > > sg_p_init
SG initial p.
int num_g_sg
Number of stochastic response vectors.
Teuchos::RCP< Stokhos::EpetraMultiVectorOrthogPoly > create_x_mv_sg_overlap(int num_vecs, Epetra_DataAccess CV=Copy, const Epetra_MultiVector *v=NULL) const
Create vector orthog poly using x map and overlap sg map.
Teuchos::RCP< Epetra_Export > interlace_overlapped_f_exporter
Exporter from SG-overlapped to SG maps.
Teuchos::RCP< const Epetra_Import > get_x_sg_importer() const
Return x sg importer.
Teuchos::RCP< const Epetra_Map > f_map
Underlying residual map.
Teuchos::RCP< const EpetraExt::MultiComm > sg_comm
Parallel SG communicator.
Teuchos::RCP< const Epetra_Map > get_f_map() const
Return residual vector map.
Teuchos::Array< Teuchos::RCP< Stokhos::EpetraMultiVectorOrthogPoly > > dgdx_dot_sg_blocks
dg/dxdot stochastic Galerkin components
Teuchos::RCP< const Stokhos::EpetraVectorOrthogPoly > get_x_sg_init() const
Return initial SG x.
Teuchos::Array< int > sg_g_index_map
Index map between block-g and g_sg maps.
Teuchos::RCP< const Epetra_Map > interlace_overlapped_f_map
Block SG overlapped residual map.
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > create_p_sg(int l, Epetra_DataAccess CV=Copy, const Epetra_Vector *v=NULL) const
Create vector orthog poly using p map.
void set_x_sg_init(const Stokhos::EpetraVectorOrthogPoly &x_sg_in)
Set initial solution polynomial.
Teuchos::Array< Teuchos::RCP< const Epetra_Map > > sg_g_map
Block SG response map.
Teuchos::Array< Teuchos::RCP< const Epetra_Map > > sg_p_map
Block SG parameter map.
Teuchos::RCP< Teuchos::ParameterList > params
Algorithmic parameters.
void set_p_sg_init(int i, const Stokhos::EpetraVectorOrthogPoly &p_sg_in)
Set initial parameter polynomial.
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > create_f_sg_overlap(Epetra_DataAccess CV=Copy, const Epetra_Vector *v=NULL) const
Create vector orthog poly using f map and overlap sg map.
Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > sg_basis
Stochastic Galerkin basis.
Teuchos::RCP< const Epetra_Map > interlace_overlapped_x_map
Block SG overlapped unknown map.
SGModelEvaluator_Interlaced(const Teuchos::RCP< EpetraExt::ModelEvaluator > &me, const Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > &sg_basis, const Teuchos::RCP< const Stokhos::Quadrature< int, double > > &sg_quad, const Teuchos::RCP< Stokhos::OrthogPolyExpansion< int, double > > &sg_exp, const Teuchos::RCP< const Stokhos::ParallelData > &sg_parallel_data, const Teuchos::RCP< Teuchos::ParameterList > &params, bool scaleOP=true)
Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > W_sg_blocks
W stochastic Galerkin components.
Teuchos::RCP< const Epetra_BlockMap > get_overlap_stochastic_map() const
Return overlap stochastic map.
Teuchos::RCP< const Epetra_Map > x_map
Underlying unknown map.
Teuchos::RCP< const Stokhos::ParallelData > sg_parallel_data
Parallel SG data.
int num_g
Number of response vectors of underlying model evaluator.
Teuchos::Array< int > get_g_sg_map_indices() const
Get indices of SG responses.
Teuchos::RCP< Stokhos::EpetraMultiVectorOrthogPoly > create_g_mv_sg(int l, int num_vecs, Epetra_DataAccess CV=Copy, const Epetra_MultiVector *v=NULL) const
Create multi-vector orthog poly using g map.
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > x_sg_blocks
x stochastic Galerkin components
unsigned int num_W_blocks
Number of W stochastic blocks (may be smaller than num_sg_blocks)
Teuchos::RCP< const Epetra_Map > interlace_x_map
Block SG unknown map.
Teuchos::RCP< const Epetra_Map > get_x_map() const
Return solution vector map.
ScalarType g(const Teuchos::Array< ScalarType > &x, const ScalarType &y)