9 #ifndef Tempus_ModelEvaluatorPairPartIMEX_CombinedFSA_impl_hpp
10 #define Tempus_ModelEvaluatorPairPartIMEX_CombinedFSA_impl_hpp
12 #include "Thyra_VectorStdOps.hpp"
13 #include "Thyra_MultiVectorStdOps.hpp"
17 template <
typename Scalar>
21 Scalar> >& forwardModel,
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) *pl = *pList;
36 pl->remove(
"Sensitivity Y Tangent Index");
46 const int sens_param_index = pl->get<
int>(
"Sensitivity Parameter Index");
47 const int num_sens_param =
49 RCP<const Thyra::VectorSpaceBase<Scalar> > explicit_y_space =
51 RCP<const Thyra::VectorSpaceBase<Scalar> > implicit_x_space =
54 multiVectorProductVectorSpace(explicit_y_space, num_sens_param + 1);
56 multiVectorProductVectorSpace(explicit_y_space, num_sens_param);
58 multiVectorProductVectorSpace(implicit_x_space, num_sens_param + 1);
64 template <
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();
80 !(getIMEXVector(z)->space()->isCompatible(
81 *(this->implicitModel_->get_x_space()))),
83 "Error - WrapperModelEvaluatorPairIMEX_CombinedFSA::initialize()\n"
84 <<
" Explicit and Implicit vector x spaces are incompatible!\n"
85 <<
" Explicit vector x space = "
86 << *(getIMEXVector(z)->space())
87 <<
"\n Implicit vector x space = "
88 << *(this->implicitModel_->get_x_space()) <<
"\n");
91 const RCP<const DMVPV> z_dmvpv = rcp_dynamic_cast<
const DMVPV>(z,
true);
92 const RCP<const Thyra::MultiVectorBase<Scalar> > z_mv =
94 RCP<const PMVB> zPVector = rcp_dynamic_cast<
const PMVB>(z_mv);
96 zPVector == Teuchos::null, std::logic_error,
97 "Error - WrapperModelEvaluatorPairPartIMEX_CombinedFSA::initialize()\n"
98 " was given a VectorBase that could not be cast to a\n"
99 " ProductVectorBase!\n");
101 int numBlocks = zPVector->productSpace()->numBlocks();
103 TEUCHOS_TEST_FOR_EXCEPTION(
104 !(0 <= this->numExplicitOnlyBlocks_ &&
105 this->numExplicitOnlyBlocks_ < numBlocks),
107 "Error - WrapperModelEvaluatorPairPartIMEX_CombinedFSA::initialize()\n"
108 "Invalid number of explicit-only blocks = "
109 << this->numExplicitOnlyBlocks_
110 <<
"\nShould be in the interval [0, numBlocks) = [0, "
111 << numBlocks <<
")\n");
114 template <
typename Scalar>
118 if (this->useImplicitModel_) {
119 if (i == this->parameterIndex_)
120 return explicit_y_dydp_prod_space_;
122 return appImplicitModel_->get_p_space(i);
125 return appExplicitModel_->get_p_space(i);
128 template <
typename Scalar>
134 using Teuchos::rcp_dynamic_cast;
136 using Thyra::multiVectorProductVector;
144 if (full == Teuchos::null)
return Teuchos::null;
146 if (this->numExplicitOnlyBlocks_ == 0)
return full;
148 const RCP<DMVPV> full_dmvpv = rcp_dynamic_cast<
DMVPV>(full,
true);
149 const RCP<MultiVectorBase<Scalar> > full_mv =
151 const RCP<PMVB> blk_full_mv = rcp_dynamic_cast<
PMVB>(full_mv,
true);
154 const int numBlocks = blk_full_mv->
productSpace()->numBlocks();
155 const int numExplicitBlocks = this->numExplicitOnlyBlocks_;
156 if (numBlocks == numExplicitBlocks + 1) {
157 const RCP<MultiVectorBase<Scalar> > imex_mv =
158 blk_full_mv->getNonconstMultiVectorBlock(numExplicitBlocks);
159 return multiVectorProductVector(imex_x_dxdp_prod_space_, imex_mv);
164 return Teuchos::null;
167 template <
typename Scalar>
173 using Teuchos::rcp_dynamic_cast;
175 using Thyra::multiVectorProductVector;
183 if (full == Teuchos::null)
return Teuchos::null;
185 if (this->numExplicitOnlyBlocks_ == 0)
return full;
187 const RCP<const DMVPV> full_dmvpv = rcp_dynamic_cast<
const DMVPV>(full,
true);
188 const RCP<const MultiVectorBase<Scalar> > full_mv =
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_x_dxdp_prod_space_, imex_mv);
204 return Teuchos::null;
207 template <
typename Scalar>
213 using Teuchos::rcp_dynamic_cast;
215 using Thyra::multiVectorProductVector;
216 using Thyra::multiVectorProductVectorSpace;
224 if (full == Teuchos::null)
return Teuchos::null;
226 if (this->numExplicitOnlyBlocks_ == 0)
return full;
228 const RCP<DMVPV> full_dmvpv = rcp_dynamic_cast<
DMVPV>(full,
true);
229 const RCP<MultiVectorBase<Scalar> > full_mv =
231 const RCP<PMVB> blk_full_mv = rcp_dynamic_cast<
PMVB>(full_mv,
true);
234 const int numExplicitBlocks = this->numExplicitOnlyBlocks_;
235 if (numExplicitBlocks == 1) {
236 const RCP<MultiVectorBase<Scalar> > explicit_mv =
238 return multiVectorProductVector(explicit_y_dydp_prod_space_, explicit_mv);
243 return Teuchos::null;
246 template <
typename Scalar>
252 using Teuchos::rcp_dynamic_cast;
254 using Thyra::multiVectorProductVector;
255 using Thyra::multiVectorProductVectorSpace;
263 if (full == Teuchos::null)
return Teuchos::null;
265 if (this->numExplicitOnlyBlocks_ == 0)
return full;
267 const RCP<const DMVPV> full_dmvpv = rcp_dynamic_cast<
const DMVPV>(full,
true);
268 const RCP<const MultiVectorBase<Scalar> > full_mv =
270 const RCP<const PMVB> blk_full_mv =
271 rcp_dynamic_cast<
const PMVB>(full_mv,
true);
274 const int numExplicitBlocks = this->numExplicitOnlyBlocks_;
275 if (numExplicitBlocks == 1) {
276 const RCP<const MultiVectorBase<Scalar> > explicit_mv =
278 return multiVectorProductVector(explicit_y_dydp_prod_space_, explicit_mv);
283 return Teuchos::null;
286 template <
typename Scalar>
290 return forwardModel_;
293 template <
typename Scalar>
298 using Teuchos::rcp_dynamic_cast;
299 using Thyra::createMember;
304 if (this->useImplicitModel_ ==
true) {
305 RCP<const Thyra::VectorBase<Scalar> > y =
306 inArgs.
get_p(this->parameterIndex_);
307 if (y != Teuchos::null) {
309 rcp_dynamic_cast<
DMVPV>(createMember(*explicit_y_dydp_prod_space_));
310 Thyra::assign(y_dydp->getNonconstMultiVector().ptr(), Scalar(0.0));
311 Thyra::assign(y_dydp->getNonconstMultiVector()->col(0).ptr(), *y);
312 inArgs.
set_p(this->parameterIndex_, y_dydp);
318 template <
typename Scalar>
326 using Teuchos::rcp_dynamic_cast;
328 const int p_index = this->parameterIndex_;
333 RCP<const Thyra::VectorBase<Scalar> > x = inArgs.
get_x();
334 RCP<Thyra::VectorBase<Scalar> > x_dot =
335 Thyra::createMember(fsaImplicitModel_->get_x_space());
336 this->timeDer_->compute(x, x_dot);
338 MEB::InArgs<Scalar> fsaImplicitInArgs(this->wrapperImplicitInArgs_);
339 MEB::OutArgs<Scalar> fsaImplicitOutArgs(this->wrapperImplicitOutArgs_);
340 fsaImplicitInArgs.set_x(x);
341 fsaImplicitInArgs.set_x_dot(x_dot);
342 for (
int i = 0; i < fsaImplicitModel_->Np(); ++i) {
344 if ((inArgs.
get_p(i) != Teuchos::null) && (i != p_index))
345 fsaImplicitInArgs.set_p(i, inArgs.
get_p(i));
352 RCP<const Thyra::VectorBase<Scalar> > y;
353 RCP<const Thyra::MultiVectorBase<Scalar> > dydp;
354 if (fsaImplicitInArgs.get_p(p_index) != Teuchos::null) {
355 RCP<const Thyra::VectorBase<Scalar> > p = fsaImplicitInArgs.get_p(p_index);
356 RCP<const Thyra::MultiVectorBase<Scalar> > p_mv =
357 rcp_dynamic_cast<
const DMVPV>(p,
true)->getMultiVector();
358 const int num_param = p_mv->domain()->dim() - 1;
360 dydp = p_mv->subView(
Range1D(1, num_param));
361 fsaImplicitInArgs.set_p(p_index, y);
363 if (use_dfdp_as_tangent_) {
364 RCP<const Thyra::VectorBase<Scalar> > dydp_vec =
365 Thyra::multiVectorProductVector(explicit_dydp_prod_space_, dydp);
366 fsaImplicitInArgs.set_p(y_tangent_index_, dydp_vec);
369 fsaImplicitOutArgs.set_f(outArgs.
get_f());
370 fsaImplicitOutArgs.set_W_op(outArgs.
get_W_op());
372 fsaImplicitModel_->evalModel(fsaImplicitInArgs, fsaImplicitOutArgs);
376 if (!use_dfdp_as_tangent_ && outArgs.
get_f() != Teuchos::null) {
377 MEB::InArgs<Scalar> appImplicitInArgs =
378 appImplicitModel_->getNominalValues();
379 RCP<const Thyra::VectorBase<Scalar> > app_x =
380 rcp_dynamic_cast<
const DMVPV>(x,
true)->getMultiVector()->
col(0);
381 RCP<const Thyra::VectorBase<Scalar> > app_x_dot =
382 rcp_dynamic_cast<
const DMVPV>(x_dot,
true)->getMultiVector()->
col(0);
383 appImplicitInArgs.set_x(app_x);
384 appImplicitInArgs.set_x_dot(app_x_dot);
385 for (
int i = 0; i < appImplicitModel_->Np(); ++i) {
386 if (i != p_index) appImplicitInArgs.set_p(i, inArgs.
get_p(i));
388 appImplicitInArgs.set_p(p_index, y);
389 if (appImplicitInArgs.supports(MEB::IN_ARG_t))
390 appImplicitInArgs.set_t(inArgs.
get_t());
391 MEB::OutArgs<Scalar> appImplicitOutArgs =
392 appImplicitModel_->createOutArgs();
393 MEB::DerivativeSupport dfdp_support =
394 appImplicitOutArgs.supports(MEB::OUT_ARG_DfDp, p_index);
396 if (dfdp_support.supports(MEB::DERIV_LINEAR_OP)) {
397 if (my_dfdp_op_ == Teuchos::null)
398 my_dfdp_op_ = appImplicitModel_->create_DfDp_op(p_index);
399 appImplicitOutArgs.set_DfDp(p_index,
400 MEB::Derivative<Scalar>(my_dfdp_op_));
403 else if (dfdp_support.supports(MEB::DERIV_MV_JACOBIAN_FORM)) {
404 if (my_dfdp_mv_ == Teuchos::null)
405 my_dfdp_mv_ = Thyra::createMembers(
406 appImplicitModel_->get_f_space(),
407 appImplicitModel_->get_p_space(p_index)->dim());
408 appImplicitOutArgs.set_DfDp(
410 MEB::Derivative<Scalar>(my_dfdp_mv_, MEB::DERIV_MV_JACOBIAN_FORM));
411 my_dfdp_op_ = my_dfdp_mv_;
414 else if (dfdp_support.supports(MEB::DERIV_MV_GRADIENT_FORM)) {
415 if (my_dfdp_mv_ == Teuchos::null)
417 Thyra::createMembers(appImplicitModel_->get_p_space(p_index),
418 appImplicitModel_->get_f_space()->dim());
419 appImplicitOutArgs.set_DfDp(
421 MEB::Derivative<Scalar>(my_dfdp_mv_, MEB::DERIV_MV_GRADIENT_FORM));
422 my_dfdp_op_ = my_dfdp_mv_;
427 "Invalid df/dp support");
429 appImplicitModel_->evalModel(appImplicitInArgs, appImplicitOutArgs);
432 RCP<DMVPV> f_dfdp = rcp_dynamic_cast<
DMVPV>(outArgs.
get_f(),
true);
434 RCP<Thyra::MultiVectorBase<Scalar> > dfdp =
435 f_dfdp->getNonconstMultiVector()->subView(
Range1D(1, num_param));
436 my_dfdp_op_->apply(trans, *dydp, dfdp.ptr(), Scalar(1.0), Scalar(1.0));
440 template <
typename Scalar>
448 pl->
set<
int>(
"Sensitivity Y Tangent Index", 3);
454 #endif // Tempus_ModelEvaluatorPairPartIMEX_CombinedFSA_impl_hpp
WrapperModelEvaluatorPairPartIMEX_CombinedFSA(const Teuchos::RCP< const WrapperModelEvaluatorPairPartIMEX_Basic< Scalar > > &forwardModel, const Teuchos::RCP< const Teuchos::ParameterList > &pList=Teuchos::null)
Constructor.
virtual void initialize()
Initialize after setting member data.
Teuchos::RCP< const DMVPVS > imex_x_dxdp_prod_space_
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)
virtual Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
static Teuchos::RCP< const Teuchos::ParameterList > getValidParameters()
virtual Teuchos::RCP< MultiVectorBase< Scalar > > getNonconstMultiVectorBlock(const int k)=0
Evaluation< VectorBase< Scalar > > get_f() const
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 MultiVectorBase< Scalar > > getMultiVectorBlock(const int k) const =0
virtual Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > getForwardModel() const
Get the underlying forward model.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
ModelEvaluator pair for implicit and explicit (IMEX) evaulations.
virtual Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_p_space(int i) const
Get the p space.
RCP< const VectorBase< Scalar > > col(Ordinal j) const
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.
CombinedForwardSensitivityModelEvaluator< Scalar > FSAME
Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > appImplicitModel_
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
void set_p(int l, const RCP< const VectorBase< Scalar > > &p_l)
Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > appExplicitModel_
Teuchos::RCP< const DMVPVS > explicit_dydp_prod_space_
Teuchos::RCP< const FSAME > fsaExplicitModel_
RCP< MultiVectorBase< Scalar > > getNonconstMultiVector()
Teuchos::RCP< const FSAME > fsaImplicitModel_
#define TEUCHOS_ASSERT(assertion_test)
bool use_dfdp_as_tangent_
Teuchos::RCP< const WrapperModelEvaluatorPairPartIMEX_Basic< Scalar > > forwardModel_
RCP< LinearOpBase< Scalar > > get_W_op() const
virtual void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
RCP< const VectorBase< Scalar > > get_x() const
Teuchos::RCP< const DMVPVS > explicit_y_dydp_prod_space_
RCP< const VectorBase< Scalar > > get_p(int l) const
virtual Teuchos::RCP< Thyra::VectorBase< Scalar > > getIMEXVector(const Teuchos::RCP< Thyra::VectorBase< Scalar > > &full) const
Extract IMEX vector from a full solution vector.