45#include "Teuchos_Assert.hpp"
46#include "EpetraExt_BlockUtility.h"
47#include "EpetraExt_BlockMultiVector.h"
57 const Teuchos::RCP<EpetraExt::ModelEvaluator>& me_,
61 const Teuchos::RCP<const Stokhos::ParallelData>& sg_parallel_data_,
62 const Teuchos::RCP<Teuchos::ParameterList>& params_,
69 num_sg_blocks(sg_basis->size()),
70 num_W_blocks(sg_basis->size()),
71 num_p_blocks(sg_basis->size()),
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()),
99 if (
x_map != Teuchos::null)
137 (*sg_x_init)[0] = *(
me->get_x_init());
148 std::string W_expansion_type =
149 params->get(
"Jacobian Expansion Type",
"Full");
150 if (W_expansion_type ==
"Linear")
155 Teuchos::RCP<Epetra_BlockMap> W_overlap_map =
172 InArgs me_inargs =
me->createInArgs();
173 OutArgs me_outargs =
me->createOutArgs();
174 num_p = me_inargs.Np();
177 for (
int i=0; i<
num_p; i++) {
178 if (me_inargs.supports(IN_ARG_p_sg, i))
188 std::string p_expansion_type =
189 params->get(
"Parameter Expansion Type",
"Full");
190 if (p_expansion_type ==
"Linear")
203 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(
206 Teuchos::RCP<const Teuchos::Array<std::string> > p_names =
208 if (p_names != Teuchos::null) {
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;
227 num_g = me_outargs.Ng();
230 for (
int i=0; i<
num_g; i++) {
231 if (me_outargs.supports(OUT_ARG_g_sg, i))
244 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(
271Teuchos::RCP<const Epetra_Map>
274 return interlace_x_map;
277Teuchos::RCP<const Epetra_Map>
280 return interlace_f_map;
283Teuchos::RCP<const Epetra_Map>
286 TEUCHOS_TEST_FOR_EXCEPTION(l < 0 || l >= num_p + num_p_sg, std::logic_error,
287 "Error! Invalid p map index " << l);
289 return me->get_p_map(l);
291 return sg_p_map[l-num_p];
293 return Teuchos::null;
296Teuchos::RCP<const Epetra_Map>
299 TEUCHOS_TEST_FOR_EXCEPTION(l < 0 || l >= num_g_sg, std::logic_error,
300 "Error! Invalid g map index " << l);
304Teuchos::RCP<const Teuchos::Array<std::string> >
307 TEUCHOS_TEST_FOR_EXCEPTION(l < 0 || l >= num_p + num_p_sg, std::logic_error,
308 "Error! Invalid p map index " << l);
310 return me->get_p_names(l);
312 return sg_p_names[l-num_p];
314 return Teuchos::null;
317Teuchos::RCP<const Epetra_Vector>
320 return sg_x_init->getBlockVector();
323Teuchos::RCP<const Epetra_Vector>
326 TEUCHOS_TEST_FOR_EXCEPTION(l < 0 || l >= num_p + num_p_sg, std::logic_error,
327 "Error! Invalid p map index " << l);
329 return me->get_p_init(l);
331 return sg_p_init[l-num_p]->getBlockVector();
333 return Teuchos::null;
336Teuchos::RCP<Epetra_Operator>
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);
348 my_W->setupOperator(W_sg_blocks);
353 return Teuchos::null;
356EpetraExt::ModelEvaluator::InArgs
360 InArgs me_inargs = me->createInArgs();
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));
378EpetraExt::ModelEvaluator::OutArgs
381 OutArgsSetup outArgs;
382 OutArgs me_outargs = me->createOutArgs();
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];
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));
411 const OutArgs& outArgs)
const
414 Teuchos::RCP<const Epetra_Vector> x;
415 if (inArgs.supports(IN_ARG_x)) {
417 if (x != Teuchos::null)
420 Teuchos::RCP<const Epetra_Vector> x_dot;
421 if (inArgs.supports(IN_ARG_x_dot))
422 x_dot = inArgs.get_x_dot();
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();
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);
442 copyToPolyOrthogVector(*overlapped_x,*x_sg_blocks);
443 me_inargs.set_x_sg(x_sg_blocks);
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);
453 copyToPolyOrthogVector(*overlapped_x_dot,*x_dot_sg_blocks);
454 me_inargs.set_x_dot_sg(x_dot_sg_blocks);
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());
466 me_inargs.set_sg_basis(sg_basis);
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());
472 me_inargs.set_sg_quadrature(sg_quad);
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());
478 me_inargs.set_sg_expansion(sg_exp);
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);
489 if (p == Teuchos::null)
490 p = sg_p_init[i]->getBlockVector();
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);
499 OutArgs me_outargs = me->createOutArgs();
502 if (f_out != Teuchos::null) {
503 me_outargs.set_f_sg(f_sg_blocks);
505 me_outargs.set_W_sg(W_sg_blocks);
509 if (W_out != Teuchos::null && !eval_W_with_f )
510 me_outargs.set_W_sg(W_sg_blocks);
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)
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)
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()));
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!");
541 for (
int i=0; i<num_g_sg; i++) {
542 int ii = sg_g_index_map[i];
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);
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);
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);
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],
574 me_outargs.set_DgDx_dot_sg(ii, SGDerivative(dgdx_dot_sg_blocks[i],
575 DERIV_TRANS_MV_BY_ROW));
578 TEUCHOS_TEST_FOR_EXCEPTION(dgdx_dot.getLinearOp() == Teuchos::null &&
579 dgdx_dot.isEmpty() ==
false,
581 "Error! Stokhos::SGModelEvaluator_Interlaced::evalModel: " <<
582 "Operator form of dg/dxdot is required!");
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);
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);
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],
607 me_outargs.set_DgDx_sg(ii, SGDerivative(dgdx_sg_blocks[i],
608 DERIV_TRANS_MV_BY_ROW));
611 TEUCHOS_TEST_FOR_EXCEPTION(dgdx.getLinearOp() == Teuchos::null &&
612 dgdx.isEmpty() ==
false,
614 "Error! Stokhos::SGModelEvaluator_Interlaced::evalModel: " <<
615 "Operator form of dg/dxdot is required!");
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)
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);
636 sg_basis, overlapped_stoch_row_map,
637 me->get_p_map(
j), product_map, sg_comm,
638 View, *(dgdp.getMultiVector())));
640 me_outargs.set_DgDp_sg(ii,
j,
641 SGDerivative(dgdp_sg,
642 dgdp.getMultiVectorOrientation()));
644 TEUCHOS_TEST_FOR_EXCEPTION(dgdp.getLinearOp() != Teuchos::null,
646 "Error! Stokhos::SGModelEvaluator_Interlaced::evalModel " <<
647 "cannot handle operator form of dg/dp!");
654 me->evalModel(me_inargs, me_outargs);
657 if ((W_out != Teuchos::null || (eval_W_with_f && f_out != Teuchos::null)) ) {
659 Teuchos::RCP<Epetra_Operator> W;
660 if (W_out != Teuchos::null)
664 Teuchos::RCP<Stokhos::SGOperator> W_sg =
665 Teuchos::rcp_dynamic_cast<Stokhos::SGOperator>(W,
true);
666 W_sg->setupOperator(W_sg_blocks);
670 if (f_out!=Teuchos::null){
672 for (
int i=0; i<sg_basis->size(); i++)
673 (*f_sg_blocks)[i].Scale(sg_basis->norm_squared(i));
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);
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);
699 *sg_x_init = x_sg_in;
702Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly>
712 *sg_p_init[i] = p_sg_in;
715Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly>
724 return sg_p_index_map;
730 return sg_g_index_map;
733Teuchos::Array< Teuchos::RCP<const Epetra_Map> >
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);
742Teuchos::RCP<const Epetra_BlockMap>
745 return overlapped_stoch_row_map;
748Teuchos::RCP<const Epetra_BlockMap>
751 return interlace_overlapped_x_map;
754Teuchos::RCP<const Epetra_Import>
757 return interlace_overlapped_x_importer;
760Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
764 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_x;
767 sg_basis, stoch_row_map, x_map, get_x_map(), sg_comm));
770 sg_basis, stoch_row_map, x_map, get_x_map(), sg_comm,
775Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
779 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_x;
782 sg_basis, overlapped_stoch_row_map, x_map,
783 get_x_sg_overlap_map(), sg_comm));
786 sg_basis, overlapped_stoch_row_map, x_map,
787 get_x_sg_overlap_map(), sg_comm, CV, *v));
791Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly>
795 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> sg_x;
798 sg_basis, stoch_row_map, x_map, get_x_map(), sg_comm,
802 sg_basis, stoch_row_map, x_map, get_x_map(), sg_comm,
807Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly>
813 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> sg_x;
816 sg_basis, overlapped_stoch_row_map, x_map,
817 get_x_sg_overlap_map(), sg_comm, num_vecs));
820 sg_basis, overlapped_stoch_row_map, x_map,
821 get_x_sg_overlap_map(), sg_comm, CV, *v));
825Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
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(),
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();
838 sg_basis, overlapped_stoch_p_map, me->get_p_map(l),
839 sg_p_map[ll], sg_comm));
842 sg_basis, overlapped_stoch_p_map, me->get_p_map(l),
843 sg_p_map[ll], sg_comm, CV, *v));
847Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
851 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_f;
854 sg_basis, stoch_row_map, f_map, interlace_f_map, sg_comm));
857 sg_basis, stoch_row_map, f_map, interlace_f_map, sg_comm,
862Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
866 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_f;
869 sg_basis, overlapped_stoch_row_map, f_map,
870 interlace_overlapped_f_map, sg_comm));
873 sg_basis, overlapped_stoch_row_map, f_map,
874 interlace_overlapped_f_map, sg_comm, CV, *v));
878Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly>
884 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> sg_f;
887 sg_basis, stoch_row_map, f_map, interlace_f_map, sg_comm,
891 sg_basis, stoch_row_map, f_map, interlace_f_map, sg_comm,
896Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly>
902 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> sg_f;
905 sg_basis, overlapped_stoch_row_map, f_map,
906 interlace_overlapped_f_map, sg_comm, num_vecs));
909 sg_basis, overlapped_stoch_row_map, f_map,
910 interlace_overlapped_f_map, sg_comm, CV, *v));
914Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
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(),
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();
927 sg_basis, overlapped_stoch_row_map,
929 sg_g_map[ll], sg_comm));
932 sg_basis, overlapped_stoch_row_map,
934 sg_g_map[ll], sg_comm, CV, *v));
938Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly>
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(),
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();
952 sg_basis, overlapped_stoch_row_map,
954 sg_g_map[ll], sg_comm, num_vecs));
957 sg_basis, overlapped_stoch_row_map,
959 sg_g_map[ll], sg_comm, CV, *v));
963Teuchos::RCP<Epetra_Map>
973 std::vector<int> interlacedUnks(stochaUnks*determUnks);
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;
979 return Teuchos::rcp(
new Epetra_Map(-1,interlacedUnks.size(),&interlacedUnks[0],0,determ_map.
Comm()));
985 std::size_t numBlocks = x_sg.
size();
986 Teuchos::RCP<const EpetraExt::BlockVector> bv_x = x_sg.
getBlockVector();
989 for(std::size_t blk=0;blk<numBlocks;blk++) {
992 for(
int dof=0;dof<v.
MyLength();dof++)
993 x[dof*numBlocks+blk] = v[dof];
1000 std::size_t numBlocks = x_sg.
size();
1001 Teuchos::RCP<EpetraExt::BlockVector> bv_x = x_sg.
getBlockVector();
1004 for(std::size_t blk=0;blk<numBlocks;blk++) {
1007 for(
int dof=0;dof<v.
MyLength();dof++)
1008 v[dof] = x[dof*numBlocks+blk];
int NumGlobalElements() const
const Epetra_Comm & Comm() const
int NumMyElements() 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.
InArgs createInArgs() const
Create InArgs.
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 > ¶ms, 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.
OutArgs createOutArgs() const
Create OutArgs.
Teuchos::RCP< const Epetra_Map > get_x_map() const
Return solution vector map.
ScalarType g(const Teuchos::Array< ScalarType > &x, const ScalarType &y)