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 const Teuchos::RCP<const Teuchos::ParameterList>& pList) :
22 forwardModel_(forwardModel),
23 use_dfdp_as_tangent_(false),
28 using Thyra::multiVectorProductVectorSpace;
30 RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
31 if (pList != Teuchos::null)
36 pl->remove(
"Sensitivity Y Tangent Index");
44 const int sens_param_index = pl->get<
int>(
"Sensitivity Parameter Index");
45 const int num_sens_param =
47 RCP<const Thyra::VectorSpaceBase<Scalar> > explicit_y_space =
49 RCP<const Thyra::VectorSpaceBase<Scalar> > implicit_x_space =
52 multiVectorProductVectorSpace(explicit_y_space, num_sens_param);
54 multiVectorProductVectorSpace(implicit_x_space, num_sens_param);
61 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();
78 TEUCHOS_TEST_FOR_EXCEPTION( !(getIMEXVector(z)->space()->isCompatible(
79 *(this->implicitModel_->get_x_space()))),
81 "Error - WrapperModelEvaluatorPairIMEX_StaggeredFSA::initialize()\n"
82 " Explicit and Implicit vector x spaces are incompatible!\n"
83 " Explicit vector x space = " << *(getIMEXVector(z)->space()) <<
"\n"
84 " Implicit vector x space = " << *(this->implicitModel_->get_x_space()) <<
88 const RCP<const DMVPV> z_dmvpv = rcp_dynamic_cast<
const DMVPV>(z,
true);
89 const RCP<const Thyra::MultiVectorBase<Scalar> > z_mv =
90 z_dmvpv->getMultiVector();
91 RCP<const PMVB> zPVector = rcp_dynamic_cast<
const PMVB>(z_mv);
92 TEUCHOS_TEST_FOR_EXCEPTION( zPVector == Teuchos::null, std::logic_error,
93 "Error - WrapperModelEvaluatorPairPartIMEX_StaggeredFSA::initialize()\n"
94 " was given a VectorBase that could not be cast to a\n"
95 " ProductVectorBase!\n");
97 int numBlocks = zPVector->productSpace()->numBlocks();
99 TEUCHOS_TEST_FOR_EXCEPTION( !(0 <= this->numExplicitOnlyBlocks_ and
100 this->numExplicitOnlyBlocks_ < numBlocks),
102 "Error - WrapperModelEvaluatorPairPartIMEX_StaggeredFSA::initialize()\n"
103 "Invalid number of explicit-only blocks = " <<
104 this->numExplicitOnlyBlocks_ <<
"\n"
105 "Should be in the interval [0, numBlocks) = [0, " << numBlocks <<
")\n");
108 template <
typename Scalar>
109 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
113 if (this->useImplicitModel_) {
114 if (i == this->parameterIndex_)
115 return explicit_dydp_prod_space_;
117 return appImplicitModel_->get_p_space(i);
120 return appExplicitModel_->get_p_space(i);
123 template <
typename Scalar>
124 Teuchos::RCP<Thyra::VectorBase<Scalar> >
129 using Teuchos::rcp_dynamic_cast;
130 using Thyra::MultiVectorBase;
131 using Thyra::VectorBase;
132 using Thyra::multiVectorProductVector;
139 if(full == Teuchos::null)
140 return Teuchos::null;
142 if (this->numExplicitOnlyBlocks_==0)
145 const RCP<DMVPV> full_dmvpv = rcp_dynamic_cast<
DMVPV>(full,
true);
146 const RCP<MultiVectorBase<Scalar> > full_mv =
147 full_dmvpv->getNonconstMultiVector();
148 const RCP<PMVB> blk_full_mv = rcp_dynamic_cast<
PMVB>(full_mv,
true);
151 const int numBlocks = blk_full_mv->productSpace()->numBlocks();
152 const int numExplicitBlocks = this->numExplicitOnlyBlocks_;
153 if (numBlocks == numExplicitBlocks+1) {
154 const RCP<MultiVectorBase<Scalar> > imex_mv =
155 blk_full_mv->getNonconstMultiVectorBlock(numExplicitBlocks);
156 return multiVectorProductVector(imex_dxdp_prod_space_, imex_mv);
160 TEUCHOS_ASSERT(
false);
161 return Teuchos::null;
164 template <
typename Scalar>
165 Teuchos::RCP<const Thyra::VectorBase<Scalar> >
167 getIMEXVector(
const Teuchos::RCP<
const Thyra::VectorBase<Scalar> > & full)
const
170 using Teuchos::rcp_dynamic_cast;
171 using Thyra::MultiVectorBase;
172 using Thyra::VectorBase;
173 using Thyra::multiVectorProductVector;
180 if(full == Teuchos::null)
181 return Teuchos::null;
183 if (this->numExplicitOnlyBlocks_==0)
186 const RCP<const DMVPV> full_dmvpv = rcp_dynamic_cast<
const DMVPV>(full,
true);
187 const RCP<const MultiVectorBase<Scalar> > full_mv =
188 full_dmvpv->getMultiVector();
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);
202 TEUCHOS_ASSERT(
false);
203 return Teuchos::null;
206 template <
typename Scalar>
207 Teuchos::RCP<Thyra::VectorBase<Scalar> >
210 const Teuchos::RCP<Thyra::VectorBase<Scalar> > & full)
const
213 using Teuchos::rcp_dynamic_cast;
214 using Thyra::MultiVectorBase;
215 using Thyra::VectorBase;
216 using Thyra::multiVectorProductVectorSpace;
217 using Thyra::multiVectorProductVector;
224 if(full == Teuchos::null)
225 return Teuchos::null;
227 if (this->numExplicitOnlyBlocks_==0)
230 const RCP<DMVPV> full_dmvpv = rcp_dynamic_cast<
DMVPV>(full,
true);
231 const RCP<MultiVectorBase<Scalar> > full_mv =
232 full_dmvpv->getNonconstMultiVector();
233 const RCP<PMVB> blk_full_mv = rcp_dynamic_cast<
PMVB>(full_mv,
true);
236 const int numExplicitBlocks = this->numExplicitOnlyBlocks_;
237 if (numExplicitBlocks == 1) {
238 const RCP<MultiVectorBase<Scalar> > explicit_mv =
239 blk_full_mv->getNonconstMultiVectorBlock(0);
240 return multiVectorProductVector(explicit_dydp_prod_space_, explicit_mv);
244 TEUCHOS_ASSERT(
false);
245 return Teuchos::null;
248 template <
typename Scalar>
249 Teuchos::RCP<const Thyra::VectorBase<Scalar> >
252 const Teuchos::RCP<
const Thyra::VectorBase<Scalar> > & full)
const
255 using Teuchos::rcp_dynamic_cast;
256 using Thyra::MultiVectorBase;
257 using Thyra::VectorBase;
258 using Thyra::multiVectorProductVectorSpace;
259 using Thyra::multiVectorProductVector;
266 if(full == Teuchos::null)
267 return Teuchos::null;
269 if (this->numExplicitOnlyBlocks_==0)
272 const RCP<const DMVPV> full_dmvpv = rcp_dynamic_cast<
const DMVPV>(full,
true);
273 const RCP<const MultiVectorBase<Scalar> > full_mv =
274 full_dmvpv->getMultiVector();
275 const RCP<const PMVB> blk_full_mv =
276 rcp_dynamic_cast<
const PMVB>(full_mv,
true);
279 const int numExplicitBlocks = this->numExplicitOnlyBlocks_;
280 if (numExplicitBlocks == 1) {
281 const RCP<const MultiVectorBase<Scalar> > explicit_mv =
282 blk_full_mv->getMultiVectorBlock(0);
283 return multiVectorProductVector(explicit_dydp_prod_space_, explicit_mv);
287 TEUCHOS_ASSERT(
false);
288 return Teuchos::null;
291 template <
typename Scalar>
292 Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >
296 return forwardModel_;
299 template <
typename Scalar>
306 t_interp_ = Teuchos::ScalarTraits<Scalar>::rmax();
307 fsaExplicitModel_->setForwardSolutionHistory(sh);
311 template <
typename Scalar>
320 fsaExplicitModel_->setForwardSolutionState(s);
321 fsaImplicitModel_->setForwardSolutionState(implicit_x_state_);
324 template <
typename Scalar>
327 setSolver(
const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
328 const bool force_W_update)
332 bool tf = forwardModel_->getUseImplicitModel();
335 fsaImplicitModel_->setSolver(solver, force_W_update);
336 nc_forwardModel->setUseImplicitModel(tf);
339 template <
typename Scalar>
340 Thyra::ModelEvaluatorBase::InArgs<Scalar>
345 using Teuchos::rcp_dynamic_cast;
346 using Thyra::createMember;
348 Thyra::ModelEvaluatorBase::InArgs<Scalar> inArgs = Base::createInArgs();
351 if (this->useImplicitModel_ ==
true) {
352 if (inArgs.get_p(this->parameterIndex_)!= Teuchos::null) {
354 rcp_dynamic_cast<
DMVPV>(createMember(*explicit_dydp_prod_space_),
true);
355 Thyra::assign(dydp->getNonconstMultiVector().ptr(), Scalar(0.0));
356 inArgs.set_p(this->parameterIndex_, dydp);
362 template <
typename Scalar>
366 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> & outArgs)
const
368 typedef Thyra::ModelEvaluatorBase MEB;
370 using Teuchos::rcp_dynamic_cast;
371 using Teuchos::Range1D;
376 if (sh_ != Teuchos::null) {
377 forward_t = inArgs.get_t();
378 if (t_interp_ != forward_t) {
379 if (nc_forward_state_ == Teuchos::null)
380 nc_forward_state_ = sh_->interpolateState(forward_t);
382 sh_->interpolateState(forward_t, nc_forward_state_.get());
383 forward_state_ = nc_forward_state_;
384 t_interp_ = forward_t;
386 fsaImplicitModel_->setForwardSolutionState(implicit_x_state_);
390 TEUCHOS_ASSERT(forward_state_ != Teuchos::null);
391 forward_t = forward_state_->getTime();
394 const int p_index = this->parameterIndex_;
399 RCP<const Thyra::VectorBase<Scalar> > x = inArgs.get_x();
400 RCP<Thyra::VectorBase<Scalar> > x_dot =
401 Thyra::createMember(fsaImplicitModel_->get_x_space());
402 this->timeDer_->compute(x, x_dot);
404 MEB::InArgs<Scalar> fsaImplicitInArgs (this->wrapperImplicitInArgs_);
405 MEB::OutArgs<Scalar> fsaImplicitOutArgs(this->wrapperImplicitOutArgs_);
406 if (fsaImplicitInArgs.supports(MEB::IN_ARG_t))
407 fsaImplicitInArgs.set_t(inArgs.get_t());
408 fsaImplicitInArgs.set_x(x);
409 fsaImplicitInArgs.set_x_dot(x_dot);
410 for (
int i=0; i<fsaImplicitModel_->Np(); ++i) {
412 if ((inArgs.get_p(i) != Teuchos::null) and (i != p_index))
413 fsaImplicitInArgs.set_p(i, inArgs.get_p(i));
420 TEUCHOS_ASSERT(explicit_y_state_ != Teuchos::null);
421 RCP<const Thyra::VectorBase<Scalar> > y =
422 explicit_y_state_->getX();
423 RCP<const Thyra::MultiVectorBase<Scalar> > dydp;
424 if (fsaImplicitInArgs.get_p(p_index) != Teuchos::null) {
425 RCP<const Thyra::VectorBase<Scalar> > p =
426 fsaImplicitInArgs.get_p(p_index);
427 dydp = rcp_dynamic_cast<
const DMVPV>(p,
true)->getMultiVector();
428 fsaImplicitInArgs.set_p(p_index, y);
430 if (use_dfdp_as_tangent_) {
431 RCP< const Thyra::VectorBase<Scalar> > dydp_vec =
432 Thyra::multiVectorProductVector(explicit_dydp_prod_space_, dydp);
433 fsaImplicitInArgs.set_p(y_tangent_index_, dydp_vec);
436 fsaImplicitOutArgs.set_f(outArgs.get_f());
437 fsaImplicitOutArgs.set_W_op(outArgs.get_W_op());
439 fsaImplicitModel_->evalModel(fsaImplicitInArgs,fsaImplicitOutArgs);
443 if (!use_dfdp_as_tangent_ && outArgs.get_f() != Teuchos::null) {
444 MEB::InArgs<Scalar> appImplicitInArgs =
445 appImplicitModel_->getNominalValues();
446 TEUCHOS_ASSERT(forward_state_ != Teuchos::null);
447 RCP< const Thyra::VectorBase<Scalar> > app_x =
448 implicit_x_state_->getX();
449 RCP< const Thyra::VectorBase<Scalar> > app_x_dot =
450 implicit_x_state_->getXDot();
451 appImplicitInArgs.set_x(app_x);
452 appImplicitInArgs.set_x_dot(app_x_dot);
453 for (
int i=0; i<appImplicitModel_->Np(); ++i) {
455 appImplicitInArgs.set_p(i, inArgs.get_p(i));
457 appImplicitInArgs.set_p(p_index, y);
458 if (appImplicitInArgs.supports(MEB::IN_ARG_t))
459 appImplicitInArgs.set_t(forward_t);
460 MEB::OutArgs<Scalar> appImplicitOutArgs =
461 appImplicitModel_->createOutArgs();
462 MEB::DerivativeSupport dfdp_support =
463 appImplicitOutArgs.supports(MEB::OUT_ARG_DfDp, p_index);
464 Thyra::EOpTransp trans = Thyra::NOTRANS;
465 if (dfdp_support.supports(MEB::DERIV_LINEAR_OP)) {
466 if (my_dfdp_op_ == Teuchos::null)
467 my_dfdp_op_ = appImplicitModel_->create_DfDp_op(p_index);
468 appImplicitOutArgs.set_DfDp(p_index,
469 MEB::Derivative<Scalar>(my_dfdp_op_));
470 trans = Thyra::NOTRANS;
472 else if (dfdp_support.supports(MEB::DERIV_MV_JACOBIAN_FORM)) {
473 if (my_dfdp_mv_ == Teuchos::null)
474 my_dfdp_mv_ = Thyra::createMembers(
475 appImplicitModel_->get_f_space(),
476 appImplicitModel_->get_p_space(p_index)->dim());
477 appImplicitOutArgs.set_DfDp(
478 p_index, MEB::Derivative<Scalar>(my_dfdp_mv_,
479 MEB::DERIV_MV_JACOBIAN_FORM));
480 my_dfdp_op_ = my_dfdp_mv_;
481 trans = Thyra::NOTRANS;
483 else if (dfdp_support.supports(MEB::DERIV_MV_GRADIENT_FORM)) {
484 if (my_dfdp_mv_ == Teuchos::null)
485 my_dfdp_mv_ = Thyra::createMembers(
486 appImplicitModel_->get_p_space(p_index),
487 appImplicitModel_->get_f_space()->dim());
488 appImplicitOutArgs.set_DfDp(
489 p_index, MEB::Derivative<Scalar>(my_dfdp_mv_,
490 MEB::DERIV_MV_GRADIENT_FORM));
491 my_dfdp_op_ = my_dfdp_mv_;
492 trans = Thyra::TRANS;
495 TEUCHOS_TEST_FOR_EXCEPTION(
496 true, std::logic_error,
"Invalid df/dp support");
498 appImplicitModel_->evalModel(appImplicitInArgs, appImplicitOutArgs);
501 RCP<Thyra::MultiVectorBase<Scalar> > dfdp =
502 rcp_dynamic_cast<
DMVPV>(outArgs.get_f(),
true)->getNonconstMultiVector();
503 my_dfdp_op_->apply(trans, *dydp, dfdp.ptr(), Scalar(1.0), Scalar(1.0));
507 template <
typename Scalar>
508 Teuchos::RCP<const Teuchos::ParameterList>
512 Teuchos::RCP<const Teuchos::ParameterList> fsa_pl =
514 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList(*fsa_pl);
515 pl->set<
int>(
"Sensitivity Y Tangent Index", 3);
519 template <
typename Scalar>
527 forward_state_->getMetaData(),
528 forwardModel_->getExplicitOnlyVector(forward_state_->getX()),
529 forwardModel_->getExplicitOnlyVector(forward_state_->getXDot()),
530 forwardModel_->getExplicitOnlyVector(forward_state_->getXDotDot()),
531 forward_state_->getStepperState(),
536 forward_state_->getMetaData(),
537 forwardModel_->getIMEXVector(forward_state_->getX()),
538 forwardModel_->getIMEXVector(forward_state_->getXDot()),
539 forwardModel_->getIMEXVector(forward_state_->getXDotDot()),
540 forward_state_->getStepperState(),
546 #endif // Tempus_ModelEvaluatorPairPartIMEX_StaggeredFSA_impl_hpp
WrapperModelEvaluatorPairPartIMEX_StaggeredFSA(const Teuchos::RCP< const WrapperModelEvaluatorPairPartIMEX_Basic< Scalar > > &forwardModel, const Teuchos::RCP< const Teuchos::ParameterList > &pList=Teuchos::null)
Constructor.
virtual Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > appImplicitModel_
virtual void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
static Teuchos::RCP< const Teuchos::ParameterList > getValidParameters()
Thyra::ProductMultiVectorBase< Scalar > PMVB
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.
ModelEvaluator pair for implicit and explicit (IMEX) evaulations.
Teuchos::RCP< FSAME > fsaExplicitModel_
StaggeredForwardSensitivityModelEvaluator< Scalar > FSAME
virtual void setUseImplicitModel(bool tf)
Set parameter to switch wrapperME base functions between explicit and implicit functions.
virtual void setForwardSolutionHistory(const Teuchos::RCP< const Tempus::SolutionHistory< Scalar > > &sh)
Set solution history from forward state evaluation (for interpolation)
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.
Teuchos::RCP< const DMVPVS > imex_dxdp_prod_space_
bool use_dfdp_as_tangent_
Thyra::DefaultMultiVectorProductVector< Scalar > DMVPV
Solution state for integrators and steppers. SolutionState contains the metadata for solutions and th...
Teuchos::RCP< FSAME > fsaImplicitModel_