9#ifndef Tempus_ModelEvaluatorPairPartIMEX_StaggeredFSA_impl_hpp
10#define Tempus_ModelEvaluatorPairPartIMEX_StaggeredFSA_impl_hpp
12#include "Thyra_VectorStdOps.hpp"
13#include "Thyra_MultiVectorStdOps.hpp"
17template <
typename Scalar>
21 const bool is_pseudotransient,
22 const Teuchos::RCP<const Teuchos::ParameterList>& pList) :
23 forwardModel_(forwardModel),
24 use_dfdp_as_tangent_(false),
29 using Thyra::multiVectorProductVectorSpace;
31 RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
32 if (pList != Teuchos::null)
37 pl->remove(
"Sensitivity Y Tangent Index");
45 const int sens_param_index = pl->get<
int>(
"Sensitivity Parameter Index");
46 const int num_sens_param =
48 RCP<const Thyra::VectorSpaceBase<Scalar> > explicit_y_space =
50 RCP<const Thyra::VectorSpaceBase<Scalar> > implicit_x_space =
53 multiVectorProductVectorSpace(explicit_y_space, num_sens_param);
55 multiVectorProductVectorSpace(implicit_x_space, num_sens_param);
62template <
typename Scalar>
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");
109template <
typename Scalar>
110Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
114 if (this->useImplicitModel_) {
115 if (i == this->parameterIndex_)
116 return explicit_dydp_prod_space_;
118 return appImplicitModel_->get_p_space(i);
121 return appExplicitModel_->get_p_space(i);
124template <
typename Scalar>
125Teuchos::RCP<Thyra::VectorBase<Scalar> >
130 using Teuchos::rcp_dynamic_cast;
133 using Thyra::multiVectorProductVector;
140 if(full == Teuchos::null)
141 return Teuchos::null;
143 if (this->numExplicitOnlyBlocks_==0)
146 const RCP<DMVPV> full_dmvpv = rcp_dynamic_cast<DMVPV>(full,
true);
147 const RCP<MultiVectorBase<Scalar> > full_mv =
148 full_dmvpv->getNonconstMultiVector();
149 const RCP<PMVB> blk_full_mv = rcp_dynamic_cast<PMVB>(full_mv,
true);
152 const int numBlocks = blk_full_mv->productSpace()->numBlocks();
153 const int numExplicitBlocks = this->numExplicitOnlyBlocks_;
154 if (numBlocks == numExplicitBlocks+1) {
155 const RCP<MultiVectorBase<Scalar> > imex_mv =
156 blk_full_mv->getNonconstMultiVectorBlock(numExplicitBlocks);
157 return multiVectorProductVector(imex_dxdp_prod_space_, imex_mv);
161 TEUCHOS_ASSERT(
false);
162 return Teuchos::null;
165template <
typename Scalar>
166Teuchos::RCP<const Thyra::VectorBase<Scalar> >
171 using Teuchos::rcp_dynamic_cast;
174 using Thyra::multiVectorProductVector;
181 if(full == Teuchos::null)
182 return Teuchos::null;
184 if (this->numExplicitOnlyBlocks_==0)
187 const RCP<const DMVPV> full_dmvpv = rcp_dynamic_cast<const DMVPV>(full,
true);
188 const RCP<const MultiVectorBase<Scalar> > full_mv =
189 full_dmvpv->getMultiVector();
190 const RCP<const PMVB> blk_full_mv =
191 rcp_dynamic_cast<const PMVB>(full_mv,
true);
194 const int numBlocks = blk_full_mv->productSpace()->numBlocks();
195 const int numExplicitBlocks = this->numExplicitOnlyBlocks_;
196 if (numBlocks == numExplicitBlocks+1) {
197 const RCP<const MultiVectorBase<Scalar> > imex_mv =
198 blk_full_mv->getMultiVectorBlock(numExplicitBlocks);
199 return multiVectorProductVector(imex_dxdp_prod_space_, imex_mv);
203 TEUCHOS_ASSERT(
false);
204 return Teuchos::null;
207template <
typename Scalar>
208Teuchos::RCP<Thyra::VectorBase<Scalar> >
214 using Teuchos::rcp_dynamic_cast;
217 using Thyra::multiVectorProductVectorSpace;
218 using Thyra::multiVectorProductVector;
225 if(full == Teuchos::null)
226 return Teuchos::null;
228 if (this->numExplicitOnlyBlocks_==0)
231 const RCP<DMVPV> full_dmvpv = rcp_dynamic_cast<DMVPV>(full,
true);
232 const RCP<MultiVectorBase<Scalar> > full_mv =
233 full_dmvpv->getNonconstMultiVector();
234 const RCP<PMVB> blk_full_mv = rcp_dynamic_cast<PMVB>(full_mv,
true);
237 const int numExplicitBlocks = this->numExplicitOnlyBlocks_;
238 if (numExplicitBlocks == 1) {
239 const RCP<MultiVectorBase<Scalar> > explicit_mv =
240 blk_full_mv->getNonconstMultiVectorBlock(0);
241 return multiVectorProductVector(explicit_dydp_prod_space_, explicit_mv);
245 TEUCHOS_ASSERT(
false);
246 return Teuchos::null;
249template <
typename Scalar>
250Teuchos::RCP<const Thyra::VectorBase<Scalar> >
256 using Teuchos::rcp_dynamic_cast;
259 using Thyra::multiVectorProductVectorSpace;
260 using Thyra::multiVectorProductVector;
267 if(full == Teuchos::null)
268 return Teuchos::null;
270 if (this->numExplicitOnlyBlocks_==0)
273 const RCP<const DMVPV> full_dmvpv = rcp_dynamic_cast<const DMVPV>(full,
true);
274 const RCP<const MultiVectorBase<Scalar> > full_mv =
275 full_dmvpv->getMultiVector();
276 const RCP<const PMVB> blk_full_mv =
277 rcp_dynamic_cast<const PMVB>(full_mv,
true);
280 const int numExplicitBlocks = this->numExplicitOnlyBlocks_;
281 if (numExplicitBlocks == 1) {
282 const RCP<const MultiVectorBase<Scalar> > explicit_mv =
283 blk_full_mv->getMultiVectorBlock(0);
284 return multiVectorProductVector(explicit_dydp_prod_space_, explicit_mv);
288 TEUCHOS_ASSERT(
false);
289 return Teuchos::null;
292template <
typename Scalar>
293Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >
297 return forwardModel_;
300template <
typename Scalar>
307 t_interp_ = Teuchos::ScalarTraits<Scalar>::rmax();
308 fsaExplicitModel_->setForwardSolutionHistory(sh);
312template <
typename Scalar>
321 fsaExplicitModel_->setForwardSolutionState(s);
322 fsaImplicitModel_->setForwardSolutionState(implicit_x_state_);
325template <
typename Scalar>
328setSolver(
const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
329 const bool force_W_update)
333 bool tf = forwardModel_->getUseImplicitModel();
334 Teuchos::RCP<WrapperModelEvaluatorPairPartIMEX_Basic<Scalar> > nc_forwardModel = Teuchos::rcp_const_cast<WrapperModelEvaluatorPairPartIMEX_Basic<Scalar> >(forwardModel_);
335 nc_forwardModel->setUseImplicitModel(
true);
336 fsaImplicitModel_->setSolver(solver, force_W_update);
337 nc_forwardModel->setUseImplicitModel(tf);
340template <
typename Scalar>
341Thyra::ModelEvaluatorBase::InArgs<Scalar>
346 using Teuchos::rcp_dynamic_cast;
347 using Thyra::createMember;
349 Thyra::ModelEvaluatorBase::InArgs<Scalar> inArgs = Base::createInArgs();
352 if (this->useImplicitModel_ ==
true) {
353 if (inArgs.get_p(this->parameterIndex_)!= Teuchos::null) {
355 rcp_dynamic_cast<DMVPV>(createMember(*explicit_dydp_prod_space_),
true);
356 Thyra::assign(dydp->getNonconstMultiVector().ptr(), Scalar(0.0));
357 inArgs.set_p(this->parameterIndex_, dydp);
363template <
typename Scalar>
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));
508template <
typename Scalar>
509Teuchos::RCP<const Teuchos::ParameterList>
513 Teuchos::RCP<const Teuchos::ParameterList> fsa_pl =
515 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList(*fsa_pl);
516 pl->set<
int>(
"Sensitivity Y Tangent Index", 3);
520template <
typename Scalar>
528 forward_state_->getMetaData(),
529 forwardModel_->getExplicitOnlyVector(forward_state_->getX()),
530 forwardModel_->getExplicitOnlyVector(forward_state_->getXDot()),
531 forwardModel_->getExplicitOnlyVector(forward_state_->getXDotDot()),
532 forward_state_->getStepperState(),
537 forward_state_->getMetaData(),
538 forwardModel_->getIMEXVector(forward_state_->getX()),
539 forwardModel_->getIMEXVector(forward_state_->getXDot()),
540 forwardModel_->getIMEXVector(forward_state_->getXDotDot()),
541 forward_state_->getStepperState(),
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Solution state for integrators and steppers. SolutionState contains the metadata for solutions and th...
static Teuchos::RCP< const Teuchos::ParameterList > getValidParameters()
ModelEvaluator pair for implicit and explicit (IMEX) evaulations.
void setup(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &explicitModel, const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &implicitModel, int numExplicitOnlyBlocks=0, int parameterIndex=-1)
Setup ME when using default constructor – for derived classes.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Teuchos::RCP< const DMVPVS > imex_dxdp_prod_space_
virtual void setForwardSolutionState(const Teuchos::RCP< const Tempus::SolutionState< Scalar > > &s)
Set solution state from forward state evaluation (for frozen state)
Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > appExplicitModel_
Teuchos::RCP< FSAME > fsaExplicitModel_
virtual Teuchos::RCP< Thyra::VectorBase< Scalar > > getExplicitOnlyVector(const Teuchos::RCP< Thyra::VectorBase< Scalar > > &full) const
Extract explicit-only vector from a full solution vector.
Teuchos::RCP< const WrapperModelEvaluatorPairPartIMEX_Basic< Scalar > > forwardModel_
virtual void setForwardSolutionHistory(const Teuchos::RCP< const Tempus::SolutionHistory< Scalar > > &sh)
Set solution history from forward state evaluation (for interpolation)
virtual void initialize()
Initialize after setting member data.
Teuchos::RCP< FSAME > fsaImplicitModel_
virtual void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
virtual void setSolver(const Teuchos::RCP< Thyra::NonlinearSolverBase< Scalar > > &solver, const bool force_W_update)
Set the solver of the underlying model if you want to reuse it.
virtual Teuchos::RCP< Thyra::VectorBase< Scalar > > getIMEXVector(const Teuchos::RCP< Thyra::VectorBase< Scalar > > &full) const
Extract IMEX vector from a full solution vector.
Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > appImplicitModel_
virtual Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > getForwardModel() const
Get the underlying forward model.
virtual Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
virtual Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_p_space(int i) const
Get the p space.
StaggeredForwardSensitivityModelEvaluator< Scalar > FSAME
Teuchos::RCP< const DMVPVS > explicit_dydp_prod_space_
void buildIMEXStates() const
Build implicit x and end explicit y states from forward_state_.
bool use_dfdp_as_tangent_
WrapperModelEvaluatorPairPartIMEX_StaggeredFSA(const Teuchos::RCP< const WrapperModelEvaluatorPairPartIMEX_Basic< Scalar > > &forwardModel, const bool is_pseudotransient, const Teuchos::RCP< const Teuchos::ParameterList > &pList=Teuchos::null)
Constructor.