Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
Stokhos_SGModelEvaluator_Adaptive.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<Stokhos::AdaptivityManager> & am,
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  bool onlyUseLinear_,int kExpOrder_,
63  const Teuchos::RCP<Teuchos::ParameterList>& params_)
64  : me(me_),
65  sg_basis(am->getMasterStochasticBasis()),
66  sg_row_dof_basis(am->getRowStochasticBasis()),
67  sg_quad(sg_quad_),
68  sg_exp(sg_exp_),
69  params(params_),
70  num_sg_blocks(sg_basis->size()),
71  num_W_blocks(sg_basis->size()),
72  num_p_blocks(sg_basis->size()),
73  supports_x(false),
74  x_map(me->get_x_map()),
75  f_map(me->get_f_map()),
76  sg_parallel_data(sg_parallel_data_),
77  sg_comm(sg_parallel_data->getMultiComm()),
78  epetraCijk(sg_parallel_data->getEpetraCijk()),
79  serialCijk(),
80  Cijk(epetraCijk->getParallelCijk()),
81  num_p(0),
82  num_p_sg(0),
83  sg_p_index_map(),
84  sg_p_map(),
85  sg_p_names(),
86  num_g(0),
87  num_g_sg(0),
88  sg_g_index_map(),
89  sg_g_map(),
90  x_dot_sg_blocks(),
91  x_sg_blocks(),
92  f_sg_blocks(),
93  W_sg_blocks(),
94  dgdx_dot_sg_blocks(),
95  dgdx_sg_blocks(),
96  sg_x_init(),
97  sg_p_init(),
98  eval_W_with_f(false),
99  kExpOrder(kExpOrder_),
100  onlyUseLinear(onlyUseLinear_),
101  scaleOP(am->isScaled()),
102  adaptMngr(am)
103 {
104  if (x_map != Teuchos::null)
105  supports_x = true;
106 
108  Teuchos::rcp(new Epetra_LocalMap(static_cast<int>(num_sg_blocks), 0, *(sg_parallel_data->getStochasticComm())));
110 
111  if (epetraCijk != Teuchos::null)
112  stoch_row_map = epetraCijk->getStochasticRowMap();
113 
114  if (supports_x) {
115 
116  // Create interlace SG x and f maps
117  adapted_x_map = adaptMngr->getAdaptedMap();
119 
120  adapted_f_map = adapted_x_map; // these are the same! No overlap possible in refinement!
122 
123  // Create importer/exporter from/to overlapped distribution
125  Teuchos::rcp(new Epetra_Import(*adapted_overlapped_x_map, *get_x_map()));
128 
129  // now we create the underlying Epetra block vectors
130  // that will be used by the model evaluator to construct
131  // the solution of the deterministic problem.
133 
134  // Create vector blocks
138 
139  // Create default sg_x_init
141  if (sg_x_init->myGID(0))
142  (*sg_x_init)[0] = *(me->get_x_init());
143 
144  // Preconditioner needs an x: This is adapted
145  my_x = Teuchos::rcp(new Epetra_Vector(*get_x_map()));
146 
147  // setup storage for W, these are blocked in Stokhos
148  // format
150 
151  // Determine W expansion type
152  std::string W_expansion_type =
153  params->get("Jacobian Expansion Type", "Full");
154  if (W_expansion_type == "Linear")
155  num_W_blocks = sg_basis->dimension() + 1;
156  else
158 
159  Teuchos::RCP<Epetra_BlockMap> W_overlap_map =
160  Teuchos::rcp(new Epetra_LocalMap(
161  static_cast<int>(num_W_blocks), 0,
162  *(sg_parallel_data->getStochasticComm())));
163  W_sg_blocks =
164  Teuchos::rcp(new Stokhos::EpetraOperatorOrthogPoly(
165  sg_basis, W_overlap_map, x_map, f_map, adapted_f_map,
166  sg_comm));
167  for (unsigned int i=0; i<num_W_blocks; i++)
168  W_sg_blocks->setCoeffPtr(i, me->create_W()); // allocate a bunch of matrices
169 
170  eval_W_with_f = params->get("Evaluate W with F", false);
171  }
172 
173  // Parameters -- The idea here is to add new parameter vectors
174  // for the stochastic Galerkin components of the parameters
175 
176  InArgs me_inargs = me->createInArgs();
177  OutArgs me_outargs = me->createOutArgs();
178  num_p = me_inargs.Np();
179 
180  // Get the p_sg's supported and build index map
181  for (int i=0; i<num_p; i++) {
182  if (me_inargs.supports(IN_ARG_p_sg, i))
183  sg_p_index_map.push_back(i);
184  }
185  num_p_sg = sg_p_index_map.size();
186 
187  sg_p_map.resize(num_p_sg);
188  sg_p_names.resize(num_p_sg);
189  sg_p_init.resize(num_p_sg);
190 
191  // Determine parameter expansion type
192  std::string p_expansion_type =
193  params->get("Parameter Expansion Type", "Full");
194  if (p_expansion_type == "Linear")
195  num_p_blocks = sg_basis->dimension() + 1;
196  else
198 
199  // Create parameter maps, names, and initial values
201  Teuchos::rcp(new Epetra_LocalMap(
202  static_cast<int>(num_p_blocks), 0,
203  *(sg_parallel_data->getStochasticComm())));
204  for (int i=0; i<num_p_sg; i++) {
205  Teuchos::RCP<const Epetra_Map> p_map = me->get_p_map(sg_p_index_map[i]);
206  sg_p_map[i] =
207  Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(
208  *p_map, *overlapped_stoch_p_map, *sg_comm));
209 
210  Teuchos::RCP<const Teuchos::Array<std::string> > p_names =
211  me->get_p_names(sg_p_index_map[i]);
212  if (p_names != Teuchos::null) {
213  sg_p_names[i] =
214  Teuchos::rcp(new Teuchos::Array<std::string>(num_sg_blocks*(p_names->size())));
215  for (int j=0; j<p_names->size(); j++) {
216  std::stringstream ss;
217  ss << (*p_names)[j] << " -- SG Coefficient " << i;
218  (*sg_p_names[i])[j] = ss.str();
219  }
220  }
221 
222  // Create default sg_p_init
224  sg_p_init[i]->init(0.0);
225  }
226 
227  // Responses -- The idea here is to add new parameter vectors
228  // for the stochastic Galerkin components of the respones
229  num_g = me_outargs.Ng();
230 
231  // Get the g_sg's supported and build index map
232  for (int i=0; i<num_g; i++) {
233  if (me_outargs.supports(OUT_ARG_g_sg, i))
234  sg_g_index_map.push_back(i);
235  }
236  num_g_sg = sg_g_index_map.size();
237 
238  sg_g_map.resize(num_g_sg);
240  dgdx_sg_blocks.resize(num_g_sg);
241 
242  // Create response maps
243  for (int i=0; i<num_g_sg; i++) {
244  Teuchos::RCP<const Epetra_Map> g_map = me->get_g_map(sg_g_index_map[i]);
245  sg_g_map[i] =
246  Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(
247  *g_map, *overlapped_stoch_row_map, *sg_comm));
248 
249  // Create dg/dxdot, dg/dx SG blocks
250  if (supports_x) {
251  dgdx_dot_sg_blocks[i] =
252  Teuchos::rcp(new Stokhos::EpetraMultiVectorOrthogPoly(
254  dgdx_sg_blocks[i] =
255  Teuchos::rcp(new Stokhos::EpetraMultiVectorOrthogPoly(
257  }
258  }
259 
260  // We don't support parallel for dgdx yet, so build a new EpetraCijk
261  if (supports_x) {
262  serialCijk =
263  Teuchos::rcp(new Stokhos::EpetraSparse3Tensor(sg_basis,
264  epetraCijk->getCijk(),
265  sg_comm,
267  }
268 
269 }
270 
272  const Teuchos::RCP<EpetraExt::ModelEvaluator>& me_,
273  const Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> >& sg_basis_,
274  const std::vector<Teuchos::RCP<const Stokhos::ProductBasis<int,double> > > & sg_row_dof_basis_,
275  const Teuchos::RCP<const Stokhos::Quadrature<int,double> >& sg_quad_,
276  const Teuchos::RCP<Stokhos::OrthogPolyExpansion<int,double> >& sg_exp_,
277  const Teuchos::RCP<const Stokhos::ParallelData>& sg_parallel_data_,
278  bool onlyUseLinear_,int kExpOrder_,
279  const Teuchos::RCP<Teuchos::ParameterList>& params_,
280  bool scaleOP_)
281  : me(me_),
282  sg_basis(sg_basis_),
283  sg_row_dof_basis(sg_row_dof_basis_),
284  sg_quad(sg_quad_),
285  sg_exp(sg_exp_),
286  params(params_),
287  num_sg_blocks(sg_basis->size()),
288  num_W_blocks(sg_basis->size()),
289  num_p_blocks(sg_basis->size()),
290  supports_x(false),
291  x_map(me->get_x_map()),
292  f_map(me->get_f_map()),
293  sg_parallel_data(sg_parallel_data_),
294  sg_comm(sg_parallel_data->getMultiComm()),
295  epetraCijk(sg_parallel_data->getEpetraCijk()),
296  serialCijk(),
297  Cijk(epetraCijk->getParallelCijk()),
298  num_p(0),
299  num_p_sg(0),
300  sg_p_index_map(),
301  sg_p_map(),
302  sg_p_names(),
303  num_g(0),
304  num_g_sg(0),
305  sg_g_index_map(),
306  sg_g_map(),
307  x_dot_sg_blocks(),
308  x_sg_blocks(),
309  f_sg_blocks(),
310  W_sg_blocks(),
311  dgdx_dot_sg_blocks(),
312  dgdx_sg_blocks(),
313  sg_x_init(),
314  sg_p_init(),
315  eval_W_with_f(false),
316  kExpOrder(kExpOrder_),
317  onlyUseLinear(onlyUseLinear_),
318  scaleOP(scaleOP_)
319 {
320  if (x_map != Teuchos::null)
321  supports_x = true;
322 
323  adaptMngr = Teuchos::rcp(new Stokhos::AdaptivityManager(Teuchos::rcp_dynamic_cast<const Stokhos::ProductBasis<int,double> >(sg_basis),
324  sg_row_dof_basis,*sg_parallel_data->getStochasticComm(),scaleOP));
325 
327  Teuchos::rcp(new Epetra_LocalMap(static_cast<int>(num_sg_blocks), 0, *(sg_parallel_data->getStochasticComm())));
329 
330  if (epetraCijk != Teuchos::null)
331  stoch_row_map = epetraCijk->getStochasticRowMap();
332 
333  if (supports_x) {
334 
335  // Create interlace SG x and f maps
336  adapted_x_map = adaptMngr->getAdaptedMap();
338 
339  adapted_f_map = adapted_x_map; // these are the same! No overlap possible in refinement!
341 
342  // Create importer/exporter from/to overlapped distribution
344  Teuchos::rcp(new Epetra_Import(*adapted_overlapped_x_map, *get_x_map()));
347 
348  // now we create the underlying Epetra block vectors
349  // that will be used by the model evaluator to construct
350  // the solution of the deterministic problem.
352 
353  // Create vector blocks
357 
358  // Create default sg_x_init
360  if (sg_x_init->myGID(0))
361  (*sg_x_init)[0] = *(me->get_x_init());
362 
363  // Preconditioner needs an x: This is adapted
364  my_x = Teuchos::rcp(new Epetra_Vector(*get_x_map()));
365 
366  // setup storage for W, these are blocked in Stokhos
367  // format
369 
370  // Determine W expansion type
371  std::string W_expansion_type =
372  params->get("Jacobian Expansion Type", "Full");
373  if (W_expansion_type == "Linear")
374  num_W_blocks = sg_basis->dimension() + 1;
375  else
377 
378  Teuchos::RCP<Epetra_BlockMap> W_overlap_map =
379  Teuchos::rcp(new Epetra_LocalMap(
380  static_cast<int>(num_W_blocks), 0,
381  *(sg_parallel_data->getStochasticComm())));
382  W_sg_blocks =
383  Teuchos::rcp(new Stokhos::EpetraOperatorOrthogPoly(
384  sg_basis, W_overlap_map, x_map, f_map, adapted_f_map,
385  sg_comm));
386  for (unsigned int i=0; i<num_W_blocks; i++)
387  W_sg_blocks->setCoeffPtr(i, me->create_W()); // allocate a bunch of matrices
388 
389  eval_W_with_f = params->get("Evaluate W with F", false);
390  }
391 
392  // Parameters -- The idea here is to add new parameter vectors
393  // for the stochastic Galerkin components of the parameters
394 
395  InArgs me_inargs = me->createInArgs();
396  OutArgs me_outargs = me->createOutArgs();
397  num_p = me_inargs.Np();
398 
399  // Get the p_sg's supported and build index map
400  for (int i=0; i<num_p; i++) {
401  if (me_inargs.supports(IN_ARG_p_sg, i))
402  sg_p_index_map.push_back(i);
403  }
404  num_p_sg = sg_p_index_map.size();
405 
406  sg_p_map.resize(num_p_sg);
407  sg_p_names.resize(num_p_sg);
408  sg_p_init.resize(num_p_sg);
409 
410  // Determine parameter expansion type
411  std::string p_expansion_type =
412  params->get("Parameter Expansion Type", "Full");
413  if (p_expansion_type == "Linear")
414  num_p_blocks = sg_basis->dimension() + 1;
415  else
417 
418  // Create parameter maps, names, and initial values
420  Teuchos::rcp(new Epetra_LocalMap(
421  static_cast<int>(num_p_blocks), 0,
422  *(sg_parallel_data->getStochasticComm())));
423  for (int i=0; i<num_p_sg; i++) {
424  Teuchos::RCP<const Epetra_Map> p_map = me->get_p_map(sg_p_index_map[i]);
425  sg_p_map[i] =
426  Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(
427  *p_map, *overlapped_stoch_p_map, *sg_comm));
428 
429  Teuchos::RCP<const Teuchos::Array<std::string> > p_names =
430  me->get_p_names(sg_p_index_map[i]);
431  if (p_names != Teuchos::null) {
432  sg_p_names[i] =
433  Teuchos::rcp(new Teuchos::Array<std::string>(num_sg_blocks*(p_names->size())));
434  for (int j=0; j<p_names->size(); j++) {
435  std::stringstream ss;
436  ss << (*p_names)[j] << " -- SG Coefficient " << i;
437  (*sg_p_names[i])[j] = ss.str();
438  }
439  }
440 
441  // Create default sg_p_init
443  sg_p_init[i]->init(0.0);
444  }
445 
446  // Responses -- The idea here is to add new parameter vectors
447  // for the stochastic Galerkin components of the respones
448  num_g = me_outargs.Ng();
449 
450  // Get the g_sg's supported and build index map
451  for (int i=0; i<num_g; i++) {
452  if (me_outargs.supports(OUT_ARG_g_sg, i))
453  sg_g_index_map.push_back(i);
454  }
455  num_g_sg = sg_g_index_map.size();
456 
457  sg_g_map.resize(num_g_sg);
459  dgdx_sg_blocks.resize(num_g_sg);
460 
461  // Create response maps
462  for (int i=0; i<num_g_sg; i++) {
463  Teuchos::RCP<const Epetra_Map> g_map = me->get_g_map(sg_g_index_map[i]);
464  sg_g_map[i] =
465  Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(
466  *g_map, *overlapped_stoch_row_map, *sg_comm));
467 
468  // Create dg/dxdot, dg/dx SG blocks
469  if (supports_x) {
470  dgdx_dot_sg_blocks[i] =
471  Teuchos::rcp(new Stokhos::EpetraMultiVectorOrthogPoly(
473  dgdx_sg_blocks[i] =
474  Teuchos::rcp(new Stokhos::EpetraMultiVectorOrthogPoly(
476  }
477  }
478 
479  // We don't support parallel for dgdx yet, so build a new EpetraCijk
480  if (supports_x) {
481  serialCijk =
482  Teuchos::rcp(new Stokhos::EpetraSparse3Tensor(sg_basis,
483  epetraCijk->getCijk(),
484  sg_comm,
486  }
487 
488 }
489 
490 // Overridden from EpetraExt::ModelEvaluator
491 
492 Teuchos::RCP<const Epetra_Map>
494 {
495  return adapted_x_map;
496 }
497 
498 Teuchos::RCP<const Epetra_Map>
500 {
501  return adapted_f_map;
502 }
503 
504 Teuchos::RCP<const Epetra_Map>
506 {
507  TEUCHOS_TEST_FOR_EXCEPTION(l < 0 || l >= num_p + num_p_sg, std::logic_error,
508  "Error! Invalid p map index " << l);
509  if (l < num_p)
510  return me->get_p_map(l);
511  else
512  return sg_p_map[l-num_p];
513 
514  return Teuchos::null;
515 }
516 
517 Teuchos::RCP<const Epetra_Map>
519 {
520  TEUCHOS_TEST_FOR_EXCEPTION(l < 0 || l >= num_g_sg, std::logic_error,
521  "Error! Invalid g map index " << l);
522  return sg_g_map[l];
523 }
524 
525 Teuchos::RCP<const Teuchos::Array<std::string> >
527 {
528  TEUCHOS_TEST_FOR_EXCEPTION(l < 0 || l >= num_p + num_p_sg, std::logic_error,
529  "Error! Invalid p map index " << l);
530  if (l < num_p)
531  return me->get_p_names(l);
532  else
533  return sg_p_names[l-num_p];
534 
535  return Teuchos::null;
536 }
537 
538 Teuchos::RCP<const Epetra_Vector>
540 {
541  // get stochastic galerking initial condition and write it out to x initial condition
542  Teuchos::RCP<Epetra_Vector> x_init = Teuchos::rcp(new Epetra_Vector(*get_x_map()));
543  adaptMngr->copyToAdaptiveVector(*get_x_sg_init(),*x_init);
544 
545  return x_init;
546 }
547 
548 Teuchos::RCP<const Epetra_Vector>
550 {
551  TEUCHOS_TEST_FOR_EXCEPTION(l < 0 || l >= num_p + num_p_sg, std::logic_error,
552  "Error! Invalid p map index " << l);
553  if (l < num_p)
554  return me->get_p_init(l);
555  else
556  return sg_p_init[l-num_p]->getBlockVector();
557 
558  return Teuchos::null;
559 }
560 
561 Teuchos::RCP<Epetra_Operator>
563 {
564  if (supports_x) {
565  Teuchos::RCP<Teuchos::ParameterList> sgOpParams =
566  Teuchos::rcp(&(params->sublist("SG Operator")), false);
567 
568  Teuchos::RCP<Epetra_CrsMatrix> W_crs
569  = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(me->create_W(), true);
570  const Epetra_CrsGraph & W_graph = W_crs->Graph();
571 
572  adaptMngr->setupWithGraph(W_graph,onlyUseLinear,kExpOrder);
573  my_W = adaptMngr->buildMatrixFromGraph();
574  adaptMngr->setupOperator(*my_W,*Cijk,*W_sg_blocks);
575 
576  return my_W;
577  }
578 
579  return Teuchos::null;
580 }
581 
582 EpetraExt::ModelEvaluator::InArgs
584 {
585  InArgsSetup inArgs;
586  InArgs me_inargs = me->createInArgs();
587 
588  inArgs.setModelEvalDescription(this->description());
589  inArgs.set_Np(num_p + num_p_sg);
590  inArgs.setSupports(IN_ARG_x_dot, me_inargs.supports(IN_ARG_x_dot_sg));
591  inArgs.setSupports(IN_ARG_x, me_inargs.supports(IN_ARG_x_sg));
592  inArgs.setSupports(IN_ARG_t, me_inargs.supports(IN_ARG_t));
593  inArgs.setSupports(IN_ARG_alpha, me_inargs.supports(IN_ARG_alpha));
594  inArgs.setSupports(IN_ARG_beta, me_inargs.supports(IN_ARG_beta));
595  inArgs.setSupports(IN_ARG_sg_basis, me_inargs.supports(IN_ARG_sg_basis));
596  inArgs.setSupports(IN_ARG_sg_quadrature,
597  me_inargs.supports(IN_ARG_sg_quadrature));
598  inArgs.setSupports(IN_ARG_sg_expansion,
599  me_inargs.supports(IN_ARG_sg_expansion));
600 
601  return inArgs;
602 }
603 
604 EpetraExt::ModelEvaluator::OutArgs
606 {
607  OutArgsSetup outArgs;
608  OutArgs me_outargs = me->createOutArgs();
609 
610  outArgs.setModelEvalDescription(this->description());
611  outArgs.set_Np_Ng(num_p+num_p_sg, num_g_sg);
612  outArgs.setSupports(OUT_ARG_f, me_outargs.supports(OUT_ARG_f_sg));
613  outArgs.setSupports(OUT_ARG_W, me_outargs.supports(OUT_ARG_W_sg));
614  outArgs.setSupports(OUT_ARG_WPrec, false);
615  for (int j=0; j<num_p; j++)
616  outArgs.setSupports(OUT_ARG_DfDp, j,
617  me_outargs.supports(OUT_ARG_DfDp_sg, j));
618  for (int i=0; i<num_g_sg; i++) {
619  int ii = sg_g_index_map[i];
620 // if (!me_outargs.supports(OUT_ARG_DgDx_dot_sg, ii).none())
621 // outArgs.setSupports(OUT_ARG_DgDx_dot, i, DERIV_LINEAR_OP);
622 // if (!me_outargs.supports(OUT_ARG_DgDx_sg, ii).none())
623 // outArgs.setSupports(OUT_ARG_DgDx, i, DERIV_LINEAR_OP);
624  for (int j=0; j<num_p; j++)
625  outArgs.setSupports(OUT_ARG_DgDp, i, j,
626  me_outargs.supports(OUT_ARG_DgDp_sg, ii, j));
627  }
628 
629  // We do not support derivatives w.r.t. the new SG parameters, so their
630  // support defaults to none.
631 
632  return outArgs;
633 }
634 
635 void
637  const OutArgs& outArgs) const
638 {
639  // Get the input arguments
640  Teuchos::RCP<const Epetra_Vector> x;
641  if (inArgs.supports(IN_ARG_x)) {
642  x = inArgs.get_x();
643  if (x != Teuchos::null)
644  *my_x = *x;
645  }
646  Teuchos::RCP<const Epetra_Vector> x_dot;
647  if (inArgs.supports(IN_ARG_x_dot))
648  x_dot = inArgs.get_x_dot();
649 
650  // Get the output arguments
651  EpetraExt::ModelEvaluator::Evaluation<Epetra_Vector> f_out;
652  if (outArgs.supports(OUT_ARG_f))
653  f_out = outArgs.get_f();
654  Teuchos::RCP<Epetra_Operator> W_out;
655  if (outArgs.supports(OUT_ARG_W))
656  W_out = outArgs.get_W();
657 
658  // Create underlying inargs
659  InArgs me_inargs = me->createInArgs();
660  if (x != Teuchos::null) {
661  // copy to a poly orthog vector
662  adaptMngr->copyFromAdaptiveVector(*x,*x_sg_blocks);
663  me_inargs.set_x_sg(x_sg_blocks);
664  }
665  if (x_dot != Teuchos::null) {
666  // copy to a poly orthog vector
667  adaptMngr->copyFromAdaptiveVector(*x_dot,*x_dot_sg_blocks);
668  me_inargs.set_x_dot_sg(x_dot_sg_blocks);
669  }
670  if (me_inargs.supports(IN_ARG_alpha))
671  me_inargs.set_alpha(inArgs.get_alpha());
672  if (me_inargs.supports(IN_ARG_beta))
673  me_inargs.set_beta(inArgs.get_beta());
674  if (me_inargs.supports(IN_ARG_t))
675  me_inargs.set_t(inArgs.get_t());
676  if (me_inargs.supports(IN_ARG_sg_basis)) {
677  if (inArgs.get_sg_basis() != Teuchos::null)
678  me_inargs.set_sg_basis(inArgs.get_sg_basis());
679  else
680  me_inargs.set_sg_basis(sg_basis);
681  }
682  if (me_inargs.supports(IN_ARG_sg_quadrature)) {
683  if (inArgs.get_sg_quadrature() != Teuchos::null)
684  me_inargs.set_sg_quadrature(inArgs.get_sg_quadrature());
685  else
686  me_inargs.set_sg_quadrature(sg_quad);
687  }
688  if (me_inargs.supports(IN_ARG_sg_expansion)) {
689  if (inArgs.get_sg_expansion() != Teuchos::null)
690  me_inargs.set_sg_expansion(inArgs.get_sg_expansion());
691  else
692  me_inargs.set_sg_expansion(sg_exp);
693  }
694 
695  // Pass parameters
696  for (int i=0; i<num_p; i++)
697  me_inargs.set_p(i, inArgs.get_p(i));
698  for (int i=0; i<num_p_sg; i++) {
699  Teuchos::RCP<const Epetra_Vector> p = inArgs.get_p(i+num_p);
700 
701  // We always need to pass in the SG parameters, so just use
702  // the initial parameters if it is NULL
703  if (p == Teuchos::null)
704  p = sg_p_init[i]->getBlockVector();
705 
706  // Convert block p to SG polynomial
707  Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> p_sg =
708  create_p_sg(sg_p_index_map[i], View, p.get());
709  me_inargs.set_p_sg(sg_p_index_map[i], p_sg);
710  }
711 
712  // Create underlying outargs
713  OutArgs me_outargs = me->createOutArgs();
714 
715  // f
716  if (f_out != Teuchos::null) {
717  me_outargs.set_f_sg(f_sg_blocks);
718  if (eval_W_with_f)
719  me_outargs.set_W_sg(W_sg_blocks);
720  }
721 
722  // W
723  if (W_out != Teuchos::null && !eval_W_with_f )
724  me_outargs.set_W_sg(W_sg_blocks);
725 
726  // df/dp
727  for (int i=0; i<num_p; i++) {
728  if (!outArgs.supports(OUT_ARG_DfDp, i).none()) {
729  Derivative dfdp = outArgs.get_DfDp(i);
730  if (dfdp.getMultiVector() != Teuchos::null) {
731  Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> dfdp_sg;
732  if (dfdp.getMultiVectorOrientation() == DERIV_MV_BY_COL)
733  dfdp_sg =
734  Teuchos::rcp(new Stokhos::EpetraMultiVectorOrthogPoly(
735  sg_basis, overlapped_stoch_row_map,
736  me->get_f_map(), sg_comm,
737  me->get_p_map(i)->NumMyElements()));
738  else if (dfdp.getMultiVectorOrientation() == DERIV_TRANS_MV_BY_ROW)
739  dfdp_sg =
740  Teuchos::rcp(new Stokhos::EpetraMultiVectorOrthogPoly(
741  sg_basis, overlapped_stoch_row_map,
742  me->get_p_map(i), sg_comm,
743  me->get_f_map()->NumMyElements()));
744  me_outargs.set_DfDp_sg(i,
745  SGDerivative(dfdp_sg,
746  dfdp.getMultiVectorOrientation()));
747  }
748  TEUCHOS_TEST_FOR_EXCEPTION(dfdp.getLinearOp() != Teuchos::null, std::logic_error,
749  "Error! Stokhos::SGModelEvaluator_Adaptive::evalModel " <<
750  "cannot handle operator form of df/dp!");
751  }
752  }
753 
754  // Responses (g, dg/dx, dg/dp, ...)
755  for (int i=0; i<num_g_sg; i++) {
756  int ii = sg_g_index_map[i];
757 
758  // g
759  Teuchos::RCP<Epetra_Vector> g = outArgs.get_g(i);
760  if (g != Teuchos::null) {
761  Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> g_sg =
762  create_g_sg(sg_g_index_map[i], View, g.get());
763  me_outargs.set_g_sg(ii, g_sg);
764  }
765 
766  // dg/dxdot
767  if (outArgs.supports(OUT_ARG_DgDx_dot, i).supports(DERIV_LINEAR_OP)) {
768  Derivative dgdx_dot = outArgs.get_DgDx_dot(i);
769  if (dgdx_dot.getLinearOp() != Teuchos::null) {
770  Teuchos::RCP<Stokhos::SGOperator> op =
771  Teuchos::rcp_dynamic_cast<Stokhos::SGOperator>(
772  dgdx_dot.getLinearOp(), true);
773  Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > sg_blocks =
774  op->getSGPolynomial();
775  if (me_outargs.supports(OUT_ARG_DgDx, ii).supports(DERIV_LINEAR_OP))
776  me_outargs.set_DgDx_dot_sg(ii, sg_blocks);
777  else {
778  for (unsigned int k=0; k<num_sg_blocks; k++) {
779  Teuchos::RCP<Epetra_MultiVector> mv =
780  Teuchos::rcp_dynamic_cast<Stokhos::EpetraMultiVectorOperator>(
781  sg_blocks->getCoeffPtr(k), true)->getMultiVector();
782  dgdx_dot_sg_blocks[i]->setCoeffPtr(k, mv);
783  }
784  if (me_outargs.supports(OUT_ARG_DgDx_dot_sg, ii).supports(DERIV_MV_BY_COL))
785  me_outargs.set_DgDx_dot_sg(ii, SGDerivative(dgdx_dot_sg_blocks[i],
786  DERIV_MV_BY_COL));
787  else
788  me_outargs.set_DgDx_dot_sg(ii, SGDerivative(dgdx_dot_sg_blocks[i],
789  DERIV_TRANS_MV_BY_ROW));
790  }
791  }
792  TEUCHOS_TEST_FOR_EXCEPTION(dgdx_dot.getLinearOp() == Teuchos::null &&
793  dgdx_dot.isEmpty() == false,
794  std::logic_error,
795  "Error! Stokhos::SGModelEvaluator_Adaptive::evalModel: " <<
796  "Operator form of dg/dxdot is required!");
797  }
798 
799  // dg/dx
800  if (outArgs.supports(OUT_ARG_DgDx, i).supports(DERIV_LINEAR_OP)) {
801  Derivative dgdx = outArgs.get_DgDx(i);
802  if (dgdx.getLinearOp() != Teuchos::null) {
803  Teuchos::RCP<Stokhos::SGOperator> op =
804  Teuchos::rcp_dynamic_cast<Stokhos::SGOperator>(
805  dgdx.getLinearOp(), true);
806  Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > sg_blocks =
807  op->getSGPolynomial();
808  if (me_outargs.supports(OUT_ARG_DgDx, ii).supports(DERIV_LINEAR_OP))
809  me_outargs.set_DgDx_sg(ii, sg_blocks);
810  else {
811  for (unsigned int k=0; k<num_sg_blocks; k++) {
812  Teuchos::RCP<Epetra_MultiVector> mv =
813  Teuchos::rcp_dynamic_cast<Stokhos::EpetraMultiVectorOperator>(
814  sg_blocks->getCoeffPtr(k), true)->getMultiVector();
815  dgdx_sg_blocks[i]->setCoeffPtr(k, mv);
816  }
817  if (me_outargs.supports(OUT_ARG_DgDx_sg, ii).supports(DERIV_MV_BY_COL))
818  me_outargs.set_DgDx_sg(ii, SGDerivative(dgdx_sg_blocks[i],
819  DERIV_MV_BY_COL));
820  else
821  me_outargs.set_DgDx_sg(ii, SGDerivative(dgdx_sg_blocks[i],
822  DERIV_TRANS_MV_BY_ROW));
823  }
824  }
825  TEUCHOS_TEST_FOR_EXCEPTION(dgdx.getLinearOp() == Teuchos::null &&
826  dgdx.isEmpty() == false,
827  std::logic_error,
828  "Error! Stokhos::SGModelEvaluator_Adaptive::evalModel: " <<
829  "Operator form of dg/dxdot is required!");
830  }
831 
832  // dg/dp
833  // Rembember, no derivatives w.r.t. sg parameters
834  for (int j=0; j<num_p; j++) {
835  if (!outArgs.supports(OUT_ARG_DgDp, i, j).none()) {
836  Derivative dgdp = outArgs.get_DgDp(i,j);
837  if (dgdp.getMultiVector() != Teuchos::null) {
838  Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> dgdp_sg;
839  if (dgdp.getMultiVectorOrientation() == DERIV_MV_BY_COL)
840  dgdp_sg =
841  Teuchos::rcp(new Stokhos::EpetraMultiVectorOrthogPoly(
842  sg_basis, overlapped_stoch_row_map,
843  me->get_g_map(ii), sg_g_map[i], sg_comm,
844  View, *(dgdp.getMultiVector())));
845  else if (dgdp.getMultiVectorOrientation() == DERIV_TRANS_MV_BY_ROW) {
846  Teuchos::RCP<const Epetra_BlockMap> product_map =
847  Teuchos::rcp(&(dgdp.getMultiVector()->Map()),false);
848  dgdp_sg =
849  Teuchos::rcp(new Stokhos::EpetraMultiVectorOrthogPoly(
850  sg_basis, overlapped_stoch_row_map,
851  me->get_p_map(j), product_map, sg_comm,
852  View, *(dgdp.getMultiVector())));
853  }
854  me_outargs.set_DgDp_sg(ii, j,
855  SGDerivative(dgdp_sg,
856  dgdp.getMultiVectorOrientation()));
857  }
858  TEUCHOS_TEST_FOR_EXCEPTION(dgdp.getLinearOp() != Teuchos::null,
859  std::logic_error,
860  "Error! Stokhos::SGModelEvaluator_Adaptive::evalModel " <<
861  "cannot handle operator form of dg/dp!");
862  }
863  }
864 
865  }
866 
867  // Compute the functions
868  me->evalModel(me_inargs, me_outargs);
869 
870  // Copy block SG components for W
871  if ((W_out != Teuchos::null || (eval_W_with_f && f_out != Teuchos::null)) ) {
872 
873  Teuchos::RCP<Epetra_Operator> W;
874  if (W_out != Teuchos::null)
875  W = W_out;
876  else
877  W = my_W;
878 
879  Teuchos::RCP<Epetra_CrsMatrix> W_sg = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(W, true);
880  adaptMngr->setupOperator(*W_sg,*Cijk,*W_sg_blocks);
881  }
882 
883  // f
884  if (f_out!=Teuchos::null){
885  if (!scaleOP)
886  for (int i=0; i<sg_basis->size(); i++)
887  (*f_sg_blocks)[i].Scale(sg_basis->norm_squared(i));
888 
889  adaptMngr->copyToAdaptiveVector(*f_sg_blocks,*f_out);
890  }
891 
892  // df/dp
893  for (int i=0; i<num_p; i++) {
894  if (!outArgs.supports(OUT_ARG_DfDp, i).none()) {
895  Derivative dfdp = outArgs.get_DfDp(i);
896  SGDerivative dfdp_sg = me_outargs.get_DfDp_sg(i);
897  if (dfdp.getMultiVector() != Teuchos::null) {
898 
899  dfdp.getMultiVector()->Export(
900  *(dfdp_sg.getMultiVector()->getBlockMultiVector()),
901  *adapted_overlapped_f_exporter, Insert);
902  }
903  }
904  }
905 }
906 
907 void
909  const Stokhos::EpetraVectorOrthogPoly& x_sg_in)
910 {
911  *sg_x_init = x_sg_in;
912 }
913 
914 Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly>
916 {
917  return sg_x_init;
918 }
919 
920 void
922  int i, const Stokhos::EpetraVectorOrthogPoly& p_sg_in)
923 {
924  *sg_p_init[i] = p_sg_in;
925 }
926 
927 Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly>
929 {
930  return sg_p_init[l];
931 }
932 
933 Teuchos::Array<int>
935 {
936  return sg_p_index_map;
937 }
938 
939 Teuchos::Array<int>
941 {
942  return sg_g_index_map;
943 }
944 
945 Teuchos::Array< Teuchos::RCP<const Epetra_Map> >
947 {
948  Teuchos::Array< Teuchos::RCP<const Epetra_Map> > base_maps(num_g);
949  for (int i=0; i<num_g; i++)
950  base_maps[i] = me->get_g_map(i);
951  return base_maps;
952  }
953 
954 Teuchos::RCP<const Epetra_BlockMap>
956 {
957  return overlapped_stoch_row_map;
958 }
959 
960 Teuchos::RCP<const Epetra_BlockMap>
962 {
963  return adapted_overlapped_x_map;
964 }
965 
966 Teuchos::RCP<const Epetra_Import>
968 {
969  return adapted_overlapped_x_importer;
970 }
971 
972 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
974 {
975  Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_x;
976  sg_x = Teuchos::rcp(new Stokhos::EpetraVectorOrthogPoly(
977  sg_basis, stoch_row_map, x_map, sg_comm));
978  return sg_x;
979 }
980 
981 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
983 {
984  return create_x_sg();
985 }
986 
987 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly>
989 {
990  Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> sg_x;
991  sg_x = Teuchos::rcp(new Stokhos::EpetraMultiVectorOrthogPoly(
992  sg_basis, stoch_row_map, x_map, sg_comm, num_vecs));
993  return sg_x;
994 }
995 
996 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly>
998  int num_vecs) const
999 {
1000  return create_x_mv_sg(num_vecs);
1001 }
1002 
1003 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
1005  const Epetra_Vector* v) const
1006 {
1007  Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_p;
1008  Teuchos::Array<int>::const_iterator it = std::find(sg_p_index_map.begin(),
1009  sg_p_index_map.end(),
1010  l);
1011  TEUCHOS_TEST_FOR_EXCEPTION(it == sg_p_index_map.end(), std::logic_error,
1012  "Error! Invalid p map index " << l);
1013  int ll = it - sg_p_index_map.begin();
1014  if (v == NULL)
1015  sg_p = Teuchos::rcp(new Stokhos::EpetraVectorOrthogPoly(
1016  sg_basis, overlapped_stoch_p_map, me->get_p_map(l),
1017  sg_p_map[ll], sg_comm));
1018  else
1019  sg_p = Teuchos::rcp(new Stokhos::EpetraVectorOrthogPoly(
1020  sg_basis, overlapped_stoch_p_map, me->get_p_map(l),
1021  sg_p_map[ll], sg_comm, CV, *v));
1022  return sg_p;
1023 }
1024 
1025 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
1027 {
1028  Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_f;
1029  sg_f = Teuchos::rcp(new Stokhos::EpetraVectorOrthogPoly(
1030  sg_basis, stoch_row_map, f_map, sg_comm));
1031  return sg_f;
1032 }
1033 
1034 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
1036 {
1037  return create_f_sg();
1038 }
1039 
1040 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly>
1042  int num_vecs) const
1043 {
1044  Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> sg_f;
1045  sg_f = Teuchos::rcp(new Stokhos::EpetraMultiVectorOrthogPoly(
1046  sg_basis, stoch_row_map, f_map, sg_comm, num_vecs));
1047  return sg_f;
1048 }
1049 
1050 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly>
1052  int num_vecs) const
1053 {
1054  return create_f_mv_sg_overlap(num_vecs);
1055 }
1056 
1057 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
1059  const Epetra_Vector* v) const
1060 {
1061  Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_g;
1062  Teuchos::Array<int>::const_iterator it = std::find(sg_g_index_map.begin(),
1063  sg_g_index_map.end(),
1064  l);
1065  TEUCHOS_TEST_FOR_EXCEPTION(it == sg_g_index_map.end(), std::logic_error,
1066  "Error! Invalid g map index " << l);
1067  int ll = it - sg_g_index_map.begin();
1068  if (v == NULL)
1069  sg_g = Teuchos::rcp(new Stokhos::EpetraVectorOrthogPoly(
1070  sg_basis, overlapped_stoch_row_map,
1071  me->get_g_map(l),
1072  sg_g_map[ll], sg_comm));
1073  else
1074  sg_g = Teuchos::rcp(new Stokhos::EpetraVectorOrthogPoly(
1075  sg_basis, overlapped_stoch_row_map,
1076  me->get_g_map(l),
1077  sg_g_map[ll], sg_comm, CV, *v));
1078  return sg_g;
1079 }
1080 
1081 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly>
1083  Epetra_DataAccess CV,
1084  const Epetra_MultiVector* v) const
1085 {
1086  Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> sg_g;
1087  Teuchos::Array<int>::const_iterator it = std::find(sg_g_index_map.begin(),
1088  sg_g_index_map.end(),
1089  l);
1090  TEUCHOS_TEST_FOR_EXCEPTION(it == sg_g_index_map.end(), std::logic_error,
1091  "Error! Invalid g map index " << l);
1092  int ll = it - sg_g_index_map.begin();
1093  if (v == NULL)
1094  sg_g = Teuchos::rcp(new Stokhos::EpetraMultiVectorOrthogPoly(
1095  sg_basis, overlapped_stoch_row_map,
1096  me->get_g_map(l),
1097  sg_g_map[ll], sg_comm, num_vecs));
1098  else
1099  sg_g = Teuchos::rcp(new Stokhos::EpetraMultiVectorOrthogPoly(
1100  sg_basis, overlapped_stoch_row_map,
1101  me->get_g_map(l),
1102  sg_g_map[ll], sg_comm, CV, *v));
1103  return sg_g;
1104 }
void evalModel(const InArgs &inArgs, const OutArgs &outArgs) const
Evaluate model on InArgs.
Teuchos::RCP< const Epetra_Map > get_g_map(int l) const
Return response map.
Teuchos::Array< int > get_p_sg_map_indices() const
Get indices of SG parameters.
bool eval_W_with_f
Whether to always evaluate W with f.
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > x_sg_blocks
x stochastic Galerkin components
A container class storing an orthogonal polynomial whose coefficients are vectors, operators, or in general any type that would have an expensive copy constructor.
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > create_x_sg() const
Create vector orthog poly using x map and owned sg map.
int num_g_sg
Number of stochastic response vectors.
Teuchos::RCP< Epetra_Operator > create_W() const
Create W = alpha*M + beta*J matrix.
Teuchos::Array< int > get_g_sg_map_indices() const
Get indices of SG responses.
Teuchos::RCP< const Epetra_Map > get_x_map() const
Return solution vector map.
Teuchos::RCP< const Epetra_BlockMap > stoch_row_map
Map for stochastic blocks.
Teuchos::RCP< const Stokhos::EpetraSparse3Tensor > serialCijk
Serial Epetra Cijk for dgdx*.
Teuchos::Array< Teuchos::RCP< Stokhos::EpetraMultiVectorOrthogPoly > > dgdx_sg_blocks
dg/dx stochastic Galerkin components
Teuchos::RCP< const Epetra_Map > get_f_map() const
Return residual vector map.
Teuchos::RCP< Epetra_Export > adapted_overlapped_f_exporter
Exporter from SG-overlapped to SG maps.
Teuchos::Array< Teuchos::RCP< Teuchos::Array< std::string > > > sg_p_names
SG coefficient parameter names.
Teuchos::Array< Teuchos::RCP< Stokhos::EpetraMultiVectorOrthogPoly > > dgdx_dot_sg_blocks
dg/dxdot stochastic Galerkin components
void set_p_sg_init(int i, const Stokhos::EpetraVectorOrthogPoly &p_sg_in)
Set initial parameter polynomial.
Insert
Teuchos::RCP< const Epetra_Map > x_map
Underlying unknown map.
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::ParallelData > sg_parallel_data
Parallel SG data.
Teuchos::RCP< const Epetra_BlockMap > get_overlap_stochastic_map() const
Return overlap stochastic map.
Teuchos::RCP< const Epetra_Map > adapted_overlapped_x_map
Block SG overlapped unknown map.
int num_g
Number of response vectors of underlying model evaluator.
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > sg_x_init
SG initial x.
Teuchos::Array< Teuchos::RCP< const Epetra_Map > > sg_g_map
Block SG response map.
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > create_f_sg_overlap() const
Create vector orthog poly using f map and overlap sg map.
Teuchos::RCP< const Epetra_Map > get_p_map(int l) const
Return parameter vector map.
Teuchos::RCP< const Stokhos::EpetraSparse3Tensor > epetraCijk
Epetra Cijk.
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int l) const
Return array of parameter names.
Abstract base class for orthogonal polynomial-based expansions.
unsigned int num_W_blocks
Number of W stochastic blocks (may be smaller than num_sg_blocks)
std::vector< Teuchos::RCP< const Stokhos::ProductBasis< int, double > > > sg_row_dof_basis
Teuchos::RCP< EpetraExt::ModelEvaluator > me
Underlying model evaluator.
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
Definition: csr_vector.h:260
An adaptor that supplies the operator interface to a multi-vector.
An abstract class to represent a generic stochastic Galerkin operator as an Epetra_Operator.
Teuchos::RCP< const Stokhos::EpetraVectorOrthogPoly > get_p_sg_init(int l) const
Return initial SG parameters.
bool supports_x
Whether we support x (and thus f and W)
Abstract base class for quadrature methods.
int num_p
Number of parameter vectors of underlying model evaluator.
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > create_p_sg(int l, Epetra_DataAccess CV=Copy, const Epetra_Vector *v=0) const
Create vector orthog poly using p map.
A container class storing an orthogonal polynomial whose coefficients are vectors, operators, or in general any type that would have an expensive copy constructor.
virtual Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > getSGPolynomial()=0
Get SG polynomial.
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > create_x_sg_overlap() const
Create vector orthog poly using x map and overlap sg map.
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > x_dot_sg_blocks
x_dot stochastic Galerkin components
Teuchos::RCP< const Epetra_Map > adapted_overlapped_f_map
Block SG overlapped residual map.
Teuchos::RCP< const Epetra_BlockMap > overlapped_stoch_row_map
Overlapped map for stochastic blocks (local map)
Teuchos::Array< Teuchos::RCP< const Epetra_Map > > get_g_sg_base_maps() const
Get base maps of SG responses.
Abstract base class for multivariate orthogonal polynomials generated from tensor products of univari...
expr expr expr expr j
Teuchos::RCP< Stokhos::EpetraMultiVectorOrthogPoly > create_x_mv_sg(int num_vecs) const
Create vector orthog poly using x map and owned sg map.
Teuchos::RCP< const Epetra_Map > adapted_f_map
Block SG residual map.
Teuchos::RCP< const Epetra_Map > adapted_x_map
Block SG unknown map.
Teuchos::RCP< const Stokhos::EpetraVectorOrthogPoly > get_x_sg_init() const
Return initial SG x.
Teuchos::Array< Teuchos::RCP< const Epetra_Map > > sg_p_map
Block SG parameter map.
Teuchos::RCP< Stokhos::EpetraMultiVectorOrthogPoly > create_f_mv_sg_overlap(int num_vecs) const
Create multi-vector orthog poly using f map and overlap sg map.
Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > W_sg_blocks
W stochastic Galerkin components.
Teuchos::RCP< Epetra_Import > adapted_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::EpetraMultiVectorOrthogPoly > create_x_mv_sg_overlap(int num_vecs) const
Create vector orthog poly using x map and overlap sg map.
Teuchos::RCP< const Epetra_Vector > get_p_init(int l) const
Return initial parameters.
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::EpetraMultiVectorOrthogPoly > create_f_mv_sg(int num_vecs) const
Create multi-vector orthog poly using f map and owned sg map.
Teuchos::Array< int > sg_p_index_map
Index map between block-p and p_sg maps.
Teuchos::RCP< const Epetra_Vector > get_x_init() const
Return initial solution.
int num_p_sg
Number of stochastic parameter vectors.
void set_x_sg_init(const Stokhos::EpetraVectorOrthogPoly &x_sg_in)
Set initial solution polynomial.
unsigned int num_sg_blocks
Number of stochastic blocks.
Teuchos::Array< int > sg_g_index_map
Index map between block-g and g_sg maps.
int Scale(double ScalarConstant)
Epetra_DataAccess
Teuchos::RCP< const EpetraExt::MultiComm > sg_comm
Parallel SG communicator.
Teuchos::RCP< const Epetra_Import > get_x_sg_importer() const
Return x sg importer.
Teuchos::RCP< Epetra_Vector > my_x
x pointer for evaluating preconditioner
Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > sg_basis
Stochastic Galerkin basis.
Teuchos::RCP< Teuchos::ParameterList > params
Algorithmic parameters.
unsigned int num_p_blocks
Number of p stochastic blocks (may be smaller than num_sg_blocks)
Teuchos::RCP< Stokhos::AdaptivityManager > adaptMngr
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > create_f_sg() const
Create vector orthog poly using f map and owned sg map.
ScalarType g(const Teuchos::Array< ScalarType > &x, const ScalarType &y)
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > f_sg_blocks
f stochastic Galerkin components
Teuchos::RCP< const Epetra_BlockMap > overlapped_stoch_p_map
Overlapped map for p stochastic blocks (local map)
SGModelEvaluator_Adaptive(const Teuchos::RCP< EpetraExt::ModelEvaluator > &me_, const Teuchos::RCP< Stokhos::AdaptivityManager > &am, 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_, bool onlyUseLinear_, int kExpOrder_, const Teuchos::RCP< Teuchos::ParameterList > &params_)
A container class storing an orthogonal polynomial whose coefficients are vectors, operators, or in general any type that would have an expensive copy constructor.
Teuchos::Array< Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > > sg_p_init
SG initial p.
Teuchos::RCP< const Epetra_Map > f_map
Underlying residual map.