9 #ifndef Tempus_ModelEvaluatorPairPartIMEX_StaggeredFSA_impl_hpp
10 #define Tempus_ModelEvaluatorPairPartIMEX_StaggeredFSA_impl_hpp
12 #include "Thyra_VectorStdOps.hpp"
13 #include "Thyra_MultiVectorStdOps.hpp"
17 template <
typename Scalar>
21 Scalar> >& forwardModel,
22 const bool is_pseudotransient,
24 : forwardModel_(forwardModel),
25 use_dfdp_as_tangent_(false),
30 using Thyra::multiVectorProductVectorSpace;
32 RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
33 if (pList != Teuchos::null) *pl = *pList;
37 pl->remove(
"Sensitivity Y Tangent Index");
47 const int sens_param_index = pl->get<
int>(
"Sensitivity Parameter Index");
48 const int num_sens_param =
50 RCP<const Thyra::VectorSpaceBase<Scalar> > explicit_y_space =
52 RCP<const Thyra::VectorSpaceBase<Scalar> > implicit_x_space =
55 multiVectorProductVectorSpace(explicit_y_space, num_sens_param);
57 multiVectorProductVectorSpace(implicit_x_space, num_sens_param);
63 template <
typename Scalar>
67 using Teuchos::rcp_dynamic_cast;
69 this->useImplicitModel_ =
true;
70 this->wrapperImplicitInArgs_ = this->createInArgs();
71 this->wrapperImplicitOutArgs_ = this->createOutArgs();
72 this->useImplicitModel_ =
false;
74 RCP<const Thyra::VectorBase<Scalar> > z =
75 this->explicitModel_->getNominalValues().get_x();
79 !(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 = "
85 << *(getIMEXVector(z)->space())
86 <<
"\n Implicit vector x space = "
87 << *(this->implicitModel_->get_x_space()) <<
"\n");
90 const RCP<const DMVPV> z_dmvpv = rcp_dynamic_cast<
const DMVPV>(z,
true);
91 const RCP<const Thyra::MultiVectorBase<Scalar> > z_mv =
93 RCP<const PMVB> zPVector = rcp_dynamic_cast<
const PMVB>(z_mv);
95 zPVector == Teuchos::null, std::logic_error,
96 "Error - WrapperModelEvaluatorPairPartIMEX_StaggeredFSA::initialize()\n"
97 " was given a VectorBase that could not be cast to a\n"
98 " ProductVectorBase!\n");
100 int numBlocks = zPVector->productSpace()->numBlocks();
102 TEUCHOS_TEST_FOR_EXCEPTION(
103 !(0 <= this->numExplicitOnlyBlocks_ &&
104 this->numExplicitOnlyBlocks_ < numBlocks),
106 "Error - WrapperModelEvaluatorPairPartIMEX_StaggeredFSA::initialize()\n"
107 <<
"Invalid number of explicit-only blocks = "
108 << this->numExplicitOnlyBlocks_
109 <<
"\nShould be in the interval [0, numBlocks) = [0, "
110 << numBlocks <<
")\n");
113 template <
typename Scalar>
117 if (this->useImplicitModel_) {
118 if (i == this->parameterIndex_)
119 return explicit_dydp_prod_space_;
121 return appImplicitModel_->get_p_space(i);
124 return appExplicitModel_->get_p_space(i);
127 template <
typename Scalar>
133 using Teuchos::rcp_dynamic_cast;
135 using Thyra::multiVectorProductVector;
143 if (full == Teuchos::null)
return Teuchos::null;
145 if (this->numExplicitOnlyBlocks_ == 0)
return full;
147 const RCP<DMVPV> full_dmvpv = rcp_dynamic_cast<
DMVPV>(full,
true);
148 const RCP<MultiVectorBase<Scalar> > full_mv =
150 const RCP<PMVB> blk_full_mv = rcp_dynamic_cast<
PMVB>(full_mv,
true);
153 const int numBlocks = blk_full_mv->
productSpace()->numBlocks();
154 const int numExplicitBlocks = this->numExplicitOnlyBlocks_;
155 if (numBlocks == numExplicitBlocks + 1) {
156 const RCP<MultiVectorBase<Scalar> > imex_mv =
157 blk_full_mv->getNonconstMultiVectorBlock(numExplicitBlocks);
158 return multiVectorProductVector(imex_dxdp_prod_space_, imex_mv);
163 return Teuchos::null;
166 template <
typename Scalar>
172 using Teuchos::rcp_dynamic_cast;
174 using Thyra::multiVectorProductVector;
182 if (full == Teuchos::null)
return Teuchos::null;
184 if (this->numExplicitOnlyBlocks_ == 0)
return full;
186 const RCP<const DMVPV> full_dmvpv = rcp_dynamic_cast<
const DMVPV>(full,
true);
187 const RCP<const MultiVectorBase<Scalar> > full_mv =
189 const RCP<const PMVB> blk_full_mv =
190 rcp_dynamic_cast<
const PMVB>(full_mv,
true);
193 const int numBlocks = blk_full_mv->
productSpace()->numBlocks();
194 const int numExplicitBlocks = this->numExplicitOnlyBlocks_;
195 if (numBlocks == numExplicitBlocks + 1) {
196 const RCP<const MultiVectorBase<Scalar> > imex_mv =
197 blk_full_mv->getMultiVectorBlock(numExplicitBlocks);
198 return multiVectorProductVector(imex_dxdp_prod_space_, imex_mv);
203 return Teuchos::null;
206 template <
typename Scalar>
212 using Teuchos::rcp_dynamic_cast;
214 using Thyra::multiVectorProductVector;
215 using Thyra::multiVectorProductVectorSpace;
223 if (full == Teuchos::null)
return Teuchos::null;
225 if (this->numExplicitOnlyBlocks_ == 0)
return full;
227 const RCP<DMVPV> full_dmvpv = rcp_dynamic_cast<
DMVPV>(full,
true);
228 const RCP<MultiVectorBase<Scalar> > full_mv =
230 const RCP<PMVB> blk_full_mv = rcp_dynamic_cast<
PMVB>(full_mv,
true);
233 const int numExplicitBlocks = this->numExplicitOnlyBlocks_;
234 if (numExplicitBlocks == 1) {
235 const RCP<MultiVectorBase<Scalar> > explicit_mv =
237 return multiVectorProductVector(explicit_dydp_prod_space_, explicit_mv);
242 return Teuchos::null;
245 template <
typename Scalar>
251 using Teuchos::rcp_dynamic_cast;
253 using Thyra::multiVectorProductVector;
254 using Thyra::multiVectorProductVectorSpace;
262 if (full == Teuchos::null)
return Teuchos::null;
264 if (this->numExplicitOnlyBlocks_ == 0)
return full;
266 const RCP<const DMVPV> full_dmvpv = rcp_dynamic_cast<
const DMVPV>(full,
true);
267 const RCP<const MultiVectorBase<Scalar> > full_mv =
269 const RCP<const PMVB> blk_full_mv =
270 rcp_dynamic_cast<
const PMVB>(full_mv,
true);
273 const int numExplicitBlocks = this->numExplicitOnlyBlocks_;
274 if (numExplicitBlocks == 1) {
275 const RCP<const MultiVectorBase<Scalar> > explicit_mv =
277 return multiVectorProductVector(explicit_dydp_prod_space_, explicit_mv);
282 return Teuchos::null;
285 template <
typename Scalar>
289 return forwardModel_;
292 template <
typename Scalar>
299 fsaExplicitModel_->setForwardSolutionHistory(sh);
303 template <
typename Scalar>
311 fsaExplicitModel_->setForwardSolutionState(s);
312 fsaImplicitModel_->setForwardSolutionState(implicit_x_state_);
315 template <
typename Scalar>
318 const bool force_W_update)
322 bool tf = forwardModel_->getUseImplicitModel();
324 nc_forwardModel = Teuchos::rcp_const_cast<
326 nc_forwardModel->setUseImplicitModel(
true);
327 fsaImplicitModel_->setSolver(solver, force_W_update);
328 nc_forwardModel->setUseImplicitModel(tf);
331 template <
typename Scalar>
336 using Teuchos::rcp_dynamic_cast;
337 using Thyra::createMember;
342 if (this->useImplicitModel_ ==
true) {
343 if (inArgs.
get_p(this->parameterIndex_) != Teuchos::null) {
344 RCP<DMVPV> dydp = rcp_dynamic_cast<
DMVPV>(
345 createMember(*explicit_dydp_prod_space_),
true);
346 Thyra::assign(dydp->getNonconstMultiVector().ptr(), Scalar(0.0));
347 inArgs.
set_p(this->parameterIndex_, dydp);
353 template <
typename Scalar>
361 using Teuchos::rcp_dynamic_cast;
366 if (sh_ != Teuchos::null) {
367 forward_t = inArgs.
get_t();
368 if (t_interp_ != forward_t) {
369 if (nc_forward_state_ == Teuchos::null)
370 nc_forward_state_ = sh_->interpolateState(forward_t);
372 sh_->interpolateState(forward_t, nc_forward_state_.get());
373 forward_state_ = nc_forward_state_;
374 t_interp_ = forward_t;
376 fsaImplicitModel_->setForwardSolutionState(implicit_x_state_);
381 forward_t = forward_state_->getTime();
384 const int p_index = this->parameterIndex_;
389 RCP<const Thyra::VectorBase<Scalar> > x = inArgs.
get_x();
390 RCP<Thyra::VectorBase<Scalar> > x_dot =
391 Thyra::createMember(fsaImplicitModel_->get_x_space());
392 this->timeDer_->compute(x, x_dot);
394 MEB::InArgs<Scalar> fsaImplicitInArgs(this->wrapperImplicitInArgs_);
395 MEB::OutArgs<Scalar> fsaImplicitOutArgs(this->wrapperImplicitOutArgs_);
396 if (fsaImplicitInArgs.supports(MEB::IN_ARG_t))
397 fsaImplicitInArgs.set_t(inArgs.
get_t());
398 fsaImplicitInArgs.set_x(x);
399 fsaImplicitInArgs.set_x_dot(x_dot);
400 for (
int i = 0; i < fsaImplicitModel_->Np(); ++i) {
402 if ((inArgs.
get_p(i) != Teuchos::null) && (i != p_index))
403 fsaImplicitInArgs.set_p(i, inArgs.
get_p(i));
411 RCP<const Thyra::VectorBase<Scalar> > y = explicit_y_state_->getX();
412 RCP<const Thyra::MultiVectorBase<Scalar> > dydp;
413 if (fsaImplicitInArgs.get_p(p_index) != Teuchos::null) {
414 RCP<const Thyra::VectorBase<Scalar> > p = fsaImplicitInArgs.get_p(p_index);
415 dydp = rcp_dynamic_cast<
const DMVPV>(p,
true)->getMultiVector();
416 fsaImplicitInArgs.set_p(p_index, y);
418 if (use_dfdp_as_tangent_) {
419 RCP<const Thyra::VectorBase<Scalar> > dydp_vec =
420 Thyra::multiVectorProductVector(explicit_dydp_prod_space_, dydp);
421 fsaImplicitInArgs.set_p(y_tangent_index_, dydp_vec);
424 fsaImplicitOutArgs.set_f(outArgs.
get_f());
425 fsaImplicitOutArgs.set_W_op(outArgs.
get_W_op());
427 fsaImplicitModel_->evalModel(fsaImplicitInArgs, fsaImplicitOutArgs);
431 if (!use_dfdp_as_tangent_ && outArgs.
get_f() != Teuchos::null) {
432 MEB::InArgs<Scalar> appImplicitInArgs =
433 appImplicitModel_->getNominalValues();
435 RCP<const Thyra::VectorBase<Scalar> > app_x = implicit_x_state_->getX();
436 RCP<const Thyra::VectorBase<Scalar> > app_x_dot =
437 implicit_x_state_->getXDot();
438 appImplicitInArgs.set_x(app_x);
439 appImplicitInArgs.set_x_dot(app_x_dot);
440 for (
int i = 0; i < appImplicitModel_->Np(); ++i) {
441 if (i != p_index) appImplicitInArgs.set_p(i, inArgs.
get_p(i));
443 appImplicitInArgs.set_p(p_index, y);
444 if (appImplicitInArgs.supports(MEB::IN_ARG_t))
445 appImplicitInArgs.set_t(forward_t);
446 MEB::OutArgs<Scalar> appImplicitOutArgs =
447 appImplicitModel_->createOutArgs();
448 MEB::DerivativeSupport dfdp_support =
449 appImplicitOutArgs.supports(MEB::OUT_ARG_DfDp, p_index);
451 if (dfdp_support.supports(MEB::DERIV_LINEAR_OP)) {
452 if (my_dfdp_op_ == Teuchos::null)
453 my_dfdp_op_ = appImplicitModel_->create_DfDp_op(p_index);
454 appImplicitOutArgs.set_DfDp(p_index,
455 MEB::Derivative<Scalar>(my_dfdp_op_));
458 else if (dfdp_support.supports(MEB::DERIV_MV_JACOBIAN_FORM)) {
459 if (my_dfdp_mv_ == Teuchos::null)
460 my_dfdp_mv_ = Thyra::createMembers(
461 appImplicitModel_->get_f_space(),
462 appImplicitModel_->get_p_space(p_index)->dim());
463 appImplicitOutArgs.set_DfDp(
465 MEB::Derivative<Scalar>(my_dfdp_mv_, MEB::DERIV_MV_JACOBIAN_FORM));
466 my_dfdp_op_ = my_dfdp_mv_;
469 else if (dfdp_support.supports(MEB::DERIV_MV_GRADIENT_FORM)) {
470 if (my_dfdp_mv_ == Teuchos::null)
472 Thyra::createMembers(appImplicitModel_->get_p_space(p_index),
473 appImplicitModel_->get_f_space()->dim());
474 appImplicitOutArgs.set_DfDp(
476 MEB::Derivative<Scalar>(my_dfdp_mv_, MEB::DERIV_MV_GRADIENT_FORM));
477 my_dfdp_op_ = my_dfdp_mv_;
482 "Invalid df/dp support");
484 appImplicitModel_->evalModel(appImplicitInArgs, appImplicitOutArgs);
487 RCP<Thyra::MultiVectorBase<Scalar> > dfdp =
488 rcp_dynamic_cast<
DMVPV>(outArgs.
get_f(),
true)
489 ->getNonconstMultiVector();
490 my_dfdp_op_->apply(trans, *dydp, dfdp.ptr(), Scalar(1.0), Scalar(1.0));
494 template <
typename Scalar>
502 pl->
set<
int>(
"Sensitivity Y Tangent Index", 3);
506 template <
typename Scalar>
511 forward_state_->getMetaData(),
512 forwardModel_->getExplicitOnlyVector(forward_state_->getX()),
513 forwardModel_->getExplicitOnlyVector(forward_state_->getXDot()),
514 forwardModel_->getExplicitOnlyVector(forward_state_->getXDotDot()),
515 forward_state_->getStepperState(), Teuchos::null));
517 forward_state_->getMetaData(),
518 forwardModel_->getIMEXVector(forward_state_->getX()),
519 forwardModel_->getIMEXVector(forward_state_->getXDot()),
520 forwardModel_->getIMEXVector(forward_state_->getXDotDot()),
521 forward_state_->getStepperState(), Teuchos::null));
526 #endif // Tempus_ModelEvaluatorPairPartIMEX_StaggeredFSA_impl_hpp
virtual Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
virtual Teuchos::RCP< const ProductVectorSpaceBase< Scalar > > productSpace() const =0
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > appImplicitModel_
virtual void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
virtual Teuchos::RCP< MultiVectorBase< Scalar > > getNonconstMultiVectorBlock(const int k)=0
Evaluation< VectorBase< Scalar > > get_f() const
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
virtual Teuchos::RCP< const MultiVectorBase< Scalar > > getMultiVectorBlock(const int k) const =0
static Teuchos::RCP< const Teuchos::ParameterList > getValidParameters()
static magnitudeType rmax()
virtual Teuchos::RCP< Thyra::VectorBase< Scalar > > getExplicitOnlyVector(const Teuchos::RCP< Thyra::VectorBase< Scalar > > &full) const
Extract explicit-only vector from a full solution vector.
virtual Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > getForwardModel() const
Get the underlying forward model.
Teuchos::RCP< const DMVPVS > explicit_dydp_prod_space_
virtual Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_p_space(int i) const
Get the p space.
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.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
ModelEvaluator pair for implicit and explicit (IMEX) evaulations.
Teuchos::RCP< FSAME > fsaExplicitModel_
StaggeredForwardSensitivityModelEvaluator< Scalar > FSAME
WrapperModelEvaluatorPairPartIMEX_StaggeredFSA(const Teuchos::RCP< const WrapperModelEvaluatorPairPartIMEX_Basic< Scalar > > &forwardModel, const bool is_pseudotransient, const Teuchos::RCP< const Teuchos::ParameterList > &pList=Teuchos::null)
Constructor.
virtual void setForwardSolutionHistory(const Teuchos::RCP< const Tempus::SolutionHistory< Scalar > > &sh)
Set solution history from forward state evaluation (for interpolation)
RCP< const MultiVectorBase< Scalar > > getMultiVector() const
void validateParametersAndSetDefaults(ParameterList const &validParamList, int const depth=1000)
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.
virtual void setForwardSolutionState(const Teuchos::RCP< const Tempus::SolutionState< Scalar > > &s)
Set solution state from forward state evaluation (for frozen state)
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
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 > > appExplicitModel_
void buildIMEXStates() const
Build implicit x and end explicit y states from forward_state_.
Teuchos::RCP< const WrapperModelEvaluatorPairPartIMEX_Basic< Scalar > > forwardModel_
virtual void initialize()
Initialize after setting member data.
void set_p(int l, const RCP< const VectorBase< Scalar > > &p_l)
Teuchos::RCP< const DMVPVS > imex_dxdp_prod_space_
RCP< MultiVectorBase< Scalar > > getNonconstMultiVector()
bool use_dfdp_as_tangent_
#define TEUCHOS_ASSERT(assertion_test)
RCP< LinearOpBase< Scalar > > get_W_op() const
RCP< const VectorBase< Scalar > > get_x() const
Solution state for integrators and steppers.
RCP< const VectorBase< Scalar > > get_p(int l) const
Teuchos::RCP< FSAME > fsaImplicitModel_