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" 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)
103 Teuchos::rcp(
new Epetra_LocalMap(
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 =
156 Teuchos::rcp(
new Epetra_LocalMap(
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")
197 Teuchos::rcp(
new Epetra_LocalMap(
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(
271 Teuchos::RCP<const Epetra_Map>
274 return interlace_x_map;
277 Teuchos::RCP<const Epetra_Map>
280 return interlace_f_map;
283 Teuchos::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;
296 Teuchos::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);
304 Teuchos::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;
317 Teuchos::RCP<const Epetra_Vector>
320 return sg_x_init->getBlockVector();
323 Teuchos::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;
336 Teuchos::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;
356 EpetraExt::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));
378 EpetraExt::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 =
558 dgdx_dot.getLinearOp(),
true);
559 Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > sg_blocks =
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 =
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 =
591 dgdx.getLinearOp(),
true);
592 Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > sg_blocks =
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 =
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 =
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;
702 Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly>
712 *sg_p_init[i] = p_sg_in;
715 Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly>
724 return sg_p_index_map;
730 return sg_g_index_map;
733 Teuchos::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);
742 Teuchos::RCP<const Epetra_BlockMap>
745 return overlapped_stoch_row_map;
748 Teuchos::RCP<const Epetra_BlockMap>
751 return interlace_overlapped_x_map;
754 Teuchos::RCP<const Epetra_Import>
757 return interlace_overlapped_x_importer;
760 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
762 const Epetra_Vector* v)
const 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,
775 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
777 const Epetra_Vector* v)
const 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));
791 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly>
793 const Epetra_MultiVector* v)
const 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,
807 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly>
810 Epetra_DataAccess CV,
811 const Epetra_MultiVector* v)
const 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));
825 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
827 const Epetra_Vector* v)
const 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));
847 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
849 const Epetra_Vector* v)
const 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,
862 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
864 const Epetra_Vector* v)
const 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));
878 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly>
881 Epetra_DataAccess CV,
882 const Epetra_MultiVector* v)
const 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,
896 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly>
899 Epetra_DataAccess CV,
900 const Epetra_MultiVector* v)
const 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));
914 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
916 const Epetra_Vector* v)
const 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));
938 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly>
940 Epetra_DataAccess CV,
941 const Epetra_MultiVector* v)
const 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));
963 Teuchos::RCP<Epetra_Map>
966 int stochaUnks = stocha_map.NumMyElements();
967 int determUnks = determ_map.NumMyElements();
970 TEUCHOS_ASSERT(stocha_map.NumGlobalElements()==stochaUnks);
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++) {
990 const Epetra_Vector & v = *bv_x->GetBlock(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++) {
1005 Epetra_Vector & v = *bv_x->GetBlock(blk);
1007 for(
int dof=0;dof<v.MyLength();dof++)
1008 v[dof] =
x[dof*numBlocks+blk];
ordinal_type size() const
Return size.
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 Epetra_Vector > get_p_init(int l) const
Return initial parameters.
Teuchos::Array< int > get_g_sg_map_indices() const
Get indices of SG responses.
Teuchos::RCP< const Stokhos::EpetraVectorOrthogPoly > get_x_sg_init() const
Return initial SG x.
virtual void setupOperator(const Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > &poly)=0
Setup operator.
void evalModel(const InArgs &inArgs, const OutArgs &outArgs) const
Evaluate model on InArgs.
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.
static Teuchos::RCP< Epetra_Map > buildInterlaceMap(const Epetra_BlockMap &determ_map, const Epetra_BlockMap &stocha_map)
Teuchos::RCP< const Epetra_Map > get_f_map() const
Return residual vector map.
Teuchos::RCP< const Epetra_Map > x_map
Underlying unknown map.
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.
bool supports_x
Whether we support x (and thus f and W)
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.
Teuchos::RCP< const Stokhos::EpetraSparse3Tensor > serialCijk
Serial Epetra Cijk for dgdx*.
Teuchos::RCP< const Epetra_Map > interlace_x_map
Block SG unknown map.
unsigned int num_sg_blocks
Number of stochastic blocks.
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > x_dot_sg_blocks
x_dot stochastic Galerkin components
Teuchos::RCP< const Epetra_Vector > get_x_init() const
Return initial solution.
Teuchos::RCP< EpetraExt::BlockVector > getBlockVector()
Get block vector.
int num_g_sg
Number of stochastic response vectors.
Teuchos::RCP< Epetra_Import > interlace_overlapped_x_importer
Importer from SG to SG-overlapped maps.
int num_g
Number of response vectors of underlying model evaluator.
OutArgs createOutArgs() const
Create OutArgs.
Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > W_sg_blocks
W stochastic Galerkin components.
Teuchos::RCP< const Epetra_Map > interlace_overlapped_f_map
Block SG overlapped residual map.
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.
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< const Stokhos::EpetraSparse3Tensor > epetraCijk
Epetra Cijk.
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.
Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > sg_basis
Stochastic Galerkin basis.
Teuchos::Array< Teuchos::RCP< Stokhos::EpetraMultiVectorOrthogPoly > > dgdx_sg_blocks
dg/dx stochastic Galerkin components
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.
Abstract base class for orthogonal polynomial-based expansions.
Teuchos::RCP< const Epetra_Import > get_x_sg_importer() const
Return x sg importer.
static void copyToPolyOrthogVector(const Epetra_Vector &x, Stokhos::EpetraVectorOrthogPoly &x_sg)
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
InArgs createInArgs() const
Create InArgs.
void set_x_sg_init(const Stokhos::EpetraVectorOrthogPoly &x_sg_in)
Set initial solution polynomial.
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 Epetra_BlockMap > overlapped_stoch_row_map
Overlapped map for stochastic blocks (local map)
int num_p
Number of parameter vectors of underlying model evaluator.
Teuchos::RCP< const Epetra_BlockMap > get_overlap_stochastic_map() const
Return overlap stochastic map.
Abstract base class for quadrature methods.
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< const Epetra_BlockMap > get_x_sg_overlap_map() const
Return x sg overlap map.
Teuchos::Array< Teuchos::RCP< Teuchos::Array< std::string > > > sg_p_names
SG coefficient parameter names.
Teuchos::RCP< const EpetraExt::MultiComm > sg_comm
Parallel SG communicator.
virtual Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > getSGPolynomial()=0
Get SG polynomial.
Teuchos::RCP< const Stokhos::ParallelData > sg_parallel_data
Parallel SG data.
Teuchos::Array< Teuchos::RCP< const Epetra_Map > > get_g_sg_base_maps() const
Get base maps of SG responses.
Teuchos::Array< Teuchos::RCP< const Epetra_Map > > sg_g_map
Block SG response map.
Teuchos::Array< Teuchos::RCP< Stokhos::EpetraMultiVectorOrthogPoly > > dgdx_dot_sg_blocks
dg/dxdot 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.
Teuchos::RCP< const Epetra_BlockMap > stoch_row_map
Map for stochastic blocks.
Teuchos::RCP< Epetra_Vector > my_x
x pointer for evaluating preconditioner
Teuchos::RCP< Epetra_Operator > create_W() const
Create W = alpha*M + beta*J matrix.
Teuchos::RCP< const Epetra_Map > interlace_overlapped_x_map
Block SG overlapped unknown map.
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > sg_x_init
SG initial x.
unsigned int num_p_blocks
Number of p stochastic blocks (may be smaller than num_sg_blocks)
Teuchos::RCP< const Stokhos::EpetraVectorOrthogPoly > get_p_sg_init(int l) const
Return initial SG parameters.
void set_p_sg_init(int i, const Stokhos::EpetraVectorOrthogPoly &p_sg_in)
Set initial parameter polynomial.
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< Epetra_Export > interlace_overlapped_f_exporter
Exporter from SG-overlapped to SG maps.
Teuchos::RCP< const Epetra_Map > interlace_f_map
Block SG residual map.
An Epetra operator representing the block stochastic Galerkin operator generated by fully assembling ...
Teuchos::RCP< const Epetra_Map > get_p_map(int l) const
Return parameter vector map.
bool eval_W_with_f
Whether to always evaluate W with f.
Teuchos::RCP< const Epetra_BlockMap > overlapped_stoch_p_map
Overlapped map for p stochastic blocks (local 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.
Teuchos::RCP< const Epetra_Map > get_g_map(int l) const
Return response map.
unsigned int num_W_blocks
Number of W stochastic blocks (may be smaller than num_sg_blocks)
Teuchos::RCP< Teuchos::ParameterList > params
Algorithmic parameters.
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::Array< Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > > sg_p_init
SG initial p.
Teuchos::Array< int > get_p_sg_map_indices() const
Get indices of SG parameters.
Teuchos::RCP< const Epetra_Map > f_map
Underlying residual map.
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< const Teuchos::Array< std::string > > get_p_names(int l) const
Return array of parameter names.
Teuchos::RCP< EpetraExt::ModelEvaluator > me
Underlying model evaluator.
Teuchos::RCP< const Epetra_Map > get_x_map() const
Return solution vector map.
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > f_sg_blocks
f stochastic Galerkin components
int num_p_sg
Number of stochastic parameter vectors.
ScalarType g(const Teuchos::Array< ScalarType > &x, const ScalarType &y)
Teuchos::Array< Teuchos::RCP< const Epetra_Map > > sg_p_map
Block SG parameter map.
Teuchos::Array< int > sg_p_index_map
Index map between block-p and p_sg maps.
A container class storing an orthogonal polynomial whose coefficients are vectors, operators, or in general any type that would have an expensive copy constructor.
static void copyToInterlacedVector(const Stokhos::EpetraVectorOrthogPoly &x_sg, Epetra_Vector &x)
Teuchos::Array< int > sg_g_index_map
Index map between block-g and g_sg maps.