68 using Teuchos::rcp_dynamic_cast;
70 this->useImplicitModel_ =
true;
71 this->wrapperImplicitInArgs_ = this->createInArgs();
72 this->wrapperImplicitOutArgs_ = this->createOutArgs();
73 this->useImplicitModel_ =
false;
75 RCP<const Thyra::VectorBase<Scalar> > z =
76 this->explicitModel_->getNominalValues().get_x();
79 TEUCHOS_TEST_FOR_EXCEPTION( !(getIMEXVector(z)->space()->isCompatible(
80 *(this->implicitModel_->get_x_space()))),
82 "Error - WrapperModelEvaluatorPairIMEX_StaggeredFSA::initialize()\n"
83 " Explicit and Implicit vector x spaces are incompatible!\n"
84 " Explicit vector x space = " << *(getIMEXVector(z)->space()) <<
"\n"
85 " Implicit vector x space = " << *(this->implicitModel_->get_x_space()) <<
89 const RCP<const DMVPV> z_dmvpv = rcp_dynamic_cast<const DMVPV>(z,
true);
90 const RCP<const Thyra::MultiVectorBase<Scalar> > z_mv =
91 z_dmvpv->getMultiVector();
92 RCP<const PMVB> zPVector = rcp_dynamic_cast<const PMVB>(z_mv);
93 TEUCHOS_TEST_FOR_EXCEPTION( zPVector == Teuchos::null, std::logic_error,
94 "Error - WrapperModelEvaluatorPairPartIMEX_StaggeredFSA::initialize()\n"
95 " was given a VectorBase that could not be cast to a\n"
96 " ProductVectorBase!\n");
98 int numBlocks = zPVector->productSpace()->numBlocks();
100 TEUCHOS_TEST_FOR_EXCEPTION( !(0 <= this->numExplicitOnlyBlocks_ &&
101 this->numExplicitOnlyBlocks_ < numBlocks),
103 "Error - WrapperModelEvaluatorPairPartIMEX_StaggeredFSA::initialize()\n"
104 "Invalid number of explicit-only blocks = " <<
105 this->numExplicitOnlyBlocks_ <<
"\n"
106 "Should be in the interval [0, numBlocks) = [0, " << numBlocks <<
")\n");
366evalModelImpl(
const Thyra::ModelEvaluatorBase::InArgs<Scalar> & inArgs,
367 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> & outArgs)
const
369 typedef Thyra::ModelEvaluatorBase MEB;
371 using Teuchos::rcp_dynamic_cast;
372 using Teuchos::Range1D;
377 if (sh_ != Teuchos::null) {
378 forward_t = inArgs.get_t();
379 if (t_interp_ != forward_t) {
380 if (nc_forward_state_ == Teuchos::null)
381 nc_forward_state_ = sh_->interpolateState(forward_t);
383 sh_->interpolateState(forward_t, nc_forward_state_.get());
384 forward_state_ = nc_forward_state_;
385 t_interp_ = forward_t;
387 fsaImplicitModel_->setForwardSolutionState(implicit_x_state_);
391 TEUCHOS_ASSERT(forward_state_ != Teuchos::null);
392 forward_t = forward_state_->getTime();
395 const int p_index = this->parameterIndex_;
400 RCP<const Thyra::VectorBase<Scalar> > x = inArgs.get_x();
401 RCP<Thyra::VectorBase<Scalar> > x_dot =
402 Thyra::createMember(fsaImplicitModel_->get_x_space());
403 this->timeDer_->compute(x, x_dot);
405 MEB::InArgs<Scalar> fsaImplicitInArgs (this->wrapperImplicitInArgs_);
406 MEB::OutArgs<Scalar> fsaImplicitOutArgs(this->wrapperImplicitOutArgs_);
407 if (fsaImplicitInArgs.supports(MEB::IN_ARG_t))
408 fsaImplicitInArgs.set_t(inArgs.get_t());
409 fsaImplicitInArgs.set_x(x);
410 fsaImplicitInArgs.set_x_dot(x_dot);
411 for (
int i=0; i<fsaImplicitModel_->Np(); ++i) {
413 if ((inArgs.get_p(i) != Teuchos::null) && (i != p_index))
414 fsaImplicitInArgs.set_p(i, inArgs.get_p(i));
421 TEUCHOS_ASSERT(explicit_y_state_ != Teuchos::null);
422 RCP<const Thyra::VectorBase<Scalar> > y =
423 explicit_y_state_->getX();
424 RCP<const Thyra::MultiVectorBase<Scalar> > dydp;
425 if (fsaImplicitInArgs.get_p(p_index) != Teuchos::null) {
426 RCP<const Thyra::VectorBase<Scalar> > p =
427 fsaImplicitInArgs.get_p(p_index);
428 dydp = rcp_dynamic_cast<const DMVPV>(p,
true)->getMultiVector();
429 fsaImplicitInArgs.set_p(p_index, y);
431 if (use_dfdp_as_tangent_) {
432 RCP< const Thyra::VectorBase<Scalar> > dydp_vec =
433 Thyra::multiVectorProductVector(explicit_dydp_prod_space_, dydp);
434 fsaImplicitInArgs.set_p(y_tangent_index_, dydp_vec);
437 fsaImplicitOutArgs.set_f(outArgs.get_f());
438 fsaImplicitOutArgs.set_W_op(outArgs.get_W_op());
440 fsaImplicitModel_->evalModel(fsaImplicitInArgs,fsaImplicitOutArgs);
444 if (!use_dfdp_as_tangent_ && outArgs.get_f() != Teuchos::null) {
445 MEB::InArgs<Scalar> appImplicitInArgs =
446 appImplicitModel_->getNominalValues();
447 TEUCHOS_ASSERT(forward_state_ != Teuchos::null);
448 RCP< const Thyra::VectorBase<Scalar> > app_x =
449 implicit_x_state_->getX();
450 RCP< const Thyra::VectorBase<Scalar> > app_x_dot =
451 implicit_x_state_->getXDot();
452 appImplicitInArgs.set_x(app_x);
453 appImplicitInArgs.set_x_dot(app_x_dot);
454 for (
int i=0; i<appImplicitModel_->Np(); ++i) {
456 appImplicitInArgs.set_p(i, inArgs.get_p(i));
458 appImplicitInArgs.set_p(p_index, y);
459 if (appImplicitInArgs.supports(MEB::IN_ARG_t))
460 appImplicitInArgs.set_t(forward_t);
461 MEB::OutArgs<Scalar> appImplicitOutArgs =
462 appImplicitModel_->createOutArgs();
463 MEB::DerivativeSupport dfdp_support =
464 appImplicitOutArgs.supports(MEB::OUT_ARG_DfDp, p_index);
465 Thyra::EOpTransp trans = Thyra::NOTRANS;
466 if (dfdp_support.supports(MEB::DERIV_LINEAR_OP)) {
467 if (my_dfdp_op_ == Teuchos::null)
468 my_dfdp_op_ = appImplicitModel_->create_DfDp_op(p_index);
469 appImplicitOutArgs.set_DfDp(p_index,
470 MEB::Derivative<Scalar>(my_dfdp_op_));
471 trans = Thyra::NOTRANS;
473 else if (dfdp_support.supports(MEB::DERIV_MV_JACOBIAN_FORM)) {
474 if (my_dfdp_mv_ == Teuchos::null)
475 my_dfdp_mv_ = Thyra::createMembers(
476 appImplicitModel_->get_f_space(),
477 appImplicitModel_->get_p_space(p_index)->dim());
478 appImplicitOutArgs.set_DfDp(
479 p_index, MEB::Derivative<Scalar>(my_dfdp_mv_,
480 MEB::DERIV_MV_JACOBIAN_FORM));
481 my_dfdp_op_ = my_dfdp_mv_;
482 trans = Thyra::NOTRANS;
484 else if (dfdp_support.supports(MEB::DERIV_MV_GRADIENT_FORM)) {
485 if (my_dfdp_mv_ == Teuchos::null)
486 my_dfdp_mv_ = Thyra::createMembers(
487 appImplicitModel_->get_p_space(p_index),
488 appImplicitModel_->get_f_space()->dim());
489 appImplicitOutArgs.set_DfDp(
490 p_index, MEB::Derivative<Scalar>(my_dfdp_mv_,
491 MEB::DERIV_MV_GRADIENT_FORM));
492 my_dfdp_op_ = my_dfdp_mv_;
493 trans = Thyra::TRANS;
496 TEUCHOS_TEST_FOR_EXCEPTION(
497 true, std::logic_error,
"Invalid df/dp support");
499 appImplicitModel_->evalModel(appImplicitInArgs, appImplicitOutArgs);
502 RCP<Thyra::MultiVectorBase<Scalar> > dfdp =
503 rcp_dynamic_cast<DMVPV>(outArgs.get_f(),
true)->getNonconstMultiVector();
504 my_dfdp_op_->apply(trans, *dydp, dfdp.ptr(), Scalar(1.0), Scalar(1.0));