9 #ifndef Tempus_AdjointAuxSensitivityModelEvaluator_impl_hpp
10 #define Tempus_AdjointAuxSensitivityModelEvaluator_impl_hpp
16 #include "Thyra_DefaultBlockedLinearOp.hpp"
18 #include "Thyra_VectorStdOps.hpp"
19 #include "Thyra_MultiVectorStdOps.hpp"
23 template <
typename Scalar>
28 const Scalar& t_init,
const Scalar& t_final,
31 adjoint_model_(adjoint_model),
34 mass_matrix_is_computed_(false),
35 t_interp_(Teuchos::ScalarTraits<Scalar>::rmax())
45 if (pList != Teuchos::null) *pl = *pList;
49 p_index_ = pl->
get<
int>(
"Sensitivity Parameter Index", 0);
50 g_index_ = pl->
get<
int>(
"Response Function Index", 0);
56 "AdjointAuxSensitivityModelEvaluator currently does not support "
57 <<
"non-constant mass matrix df/dx_dot!");
65 Array<RCP<const VectorSpaceBase<Scalar> > > x_spaces(2);
66 Array<RCP<const VectorSpaceBase<Scalar> > > f_spaces(2);
75 MEB::InArgs<Scalar> me_inArgs =
model_->createInArgs();
78 MEB::InArgsSetup<Scalar> inArgs;
79 inArgs.setModelEvalDescription(this->
description());
80 inArgs.setSupports(MEB::IN_ARG_x);
81 inArgs.setSupports(MEB::IN_ARG_t);
82 if (me_inArgs.supports(MEB::IN_ARG_x_dot))
83 inArgs.setSupports(MEB::IN_ARG_x_dot);
84 inArgs.setSupports(MEB::IN_ARG_alpha);
85 inArgs.setSupports(MEB::IN_ARG_beta);
88 inArgs.set_Np(me_inArgs.Np());
91 MEB::OutArgs<Scalar> me_outArgs =
model_->createOutArgs();
92 MEB::OutArgs<Scalar> adj_me_outArgs =
adjoint_model_->createOutArgs();
93 MEB::OutArgsSetup<Scalar> outArgs;
94 outArgs.setModelEvalDescription(this->
description());
95 outArgs.set_Np_Ng(me_inArgs.Np(), 0);
96 outArgs.setSupports(MEB::OUT_ARG_f);
97 outArgs.setSupports(MEB::OUT_ARG_W_op);
104 if (me_inArgs.supports(MEB::IN_ARG_x_dot)) {
110 template <
typename Scalar>
112 const Scalar t_final)
117 template <
typename Scalar>
123 forward_state_ = Teuchos::null;
126 template <
typename Scalar>
131 return model_->get_p_space(p);
134 template <
typename Scalar>
139 return model_->get_p_names(p);
142 template <
typename Scalar>
146 return x_prod_space_;
149 template <
typename Scalar>
153 return f_prod_space_;
156 template <
typename Scalar>
163 RCP<LinearOpBase<Scalar> > adjoint_op = adjoint_model_->create_W_op();
164 RCP<LinearOpBase<Scalar> > mv_adjoint_op =
165 Thyra::nonconstMultiVectorLinearOp(adjoint_op, num_adjoint_);
166 RCP<const Thyra::VectorSpaceBase<Scalar> > g_space = response_space_;
167 RCP<LinearOpBase<Scalar> > g_op = Thyra::scaledIdentity(g_space, Scalar(1.0));
168 RCP<LinearOpBase<Scalar> > null_op;
169 return nonconstBlock2x2(mv_adjoint_op, null_op, null_op, g_op);
172 template <
typename Scalar>
177 using Teuchos::rcp_dynamic_cast;
180 RCP<const LOWSFB> alowsfb = adjoint_model_->get_W_factory();
181 if (alowsfb == Teuchos::null)
182 return Teuchos::null;
184 RCP<const LOWSFB> mv_alowsfb = Thyra::multiVectorLinearOpWithSolveFactory(
185 alowsfb, residual_space_, adjoint_space_);
187 RCP<const Thyra::VectorSpaceBase<Scalar> > g_space = response_space_;
188 RCP<const LOWSFB> g_lowsfb =
189 Thyra::scaledIdentitySolveFactory(g_space, Scalar(1.0));
192 lowsfbs[0] = mv_alowsfb;
193 lowsfbs[1] = g_lowsfb;
194 return Thyra::blockedTriangularLinearOpWithSolveFactory(lowsfbs);
197 template <
typename Scalar>
201 return prototypeInArgs_;
204 template <
typename Scalar>
210 using Teuchos::rcp_dynamic_cast;
212 MEB::InArgs<Scalar> me_nominal = model_->getNominalValues();
213 MEB::InArgs<Scalar> nominal = this->createInArgs();
218 RCP<Thyra::VectorBase<Scalar> > x = Thyra::createMember(*x_prod_space_);
219 Thyra::assign(x.ptr(), zero);
222 if (me_nominal.supports(MEB::IN_ARG_x_dot)) {
223 RCP<Thyra::VectorBase<Scalar> > x_dot = Thyra::createMember(*x_prod_space_);
224 Thyra::assign(x_dot.ptr(), zero);
225 nominal.set_x_dot(x_dot);
228 const int np = model_->Np();
229 for (
int i = 0; i < np; ++i) nominal.set_p(i, me_nominal.get_p(i));
234 template <
typename Scalar>
238 return prototypeOutArgs_;
241 template <
typename Scalar>
248 using Teuchos::rcp_dynamic_cast;
258 const Scalar t = inArgs.
get_t();
259 const Scalar forward_t = t_final_ + t_init_ - t;
260 if (forward_state_ == Teuchos::null || t_interp_ != t) {
261 if (forward_state_ == Teuchos::null)
262 forward_state_ = sh_->interpolateState(forward_t);
264 sh_->interpolateState(forward_t, forward_state_.get());
269 MEB::InArgs<Scalar> me_inArgs = model_->getNominalValues();
270 me_inArgs.set_x(forward_state_->getX());
271 if (me_inArgs.supports(MEB::IN_ARG_x_dot) &&
273 me_inArgs.set_x_dot(forward_state_->getXDot());
274 if (me_inArgs.supports(MEB::IN_ARG_t)) me_inArgs.set_t(forward_t);
275 const int np = me_inArgs.Np();
276 for (
int i = 0; i < np; ++i) me_inArgs.set_p(i, inArgs.
get_p(i));
279 RCP<Thyra::LinearOpBase<Scalar> > op = outArgs.
get_W_op();
280 if (op != Teuchos::null) {
281 if (me_inArgs.supports(MEB::IN_ARG_alpha))
283 if (me_inArgs.supports(MEB::IN_ARG_beta))
284 me_inArgs.set_beta(inArgs.
get_beta());
287 RCP<Thyra::DefaultBlockedLinearOp<Scalar> > block_op =
289 RCP<Thyra::MultiVectorLinearOp<Scalar> > mv_adjoint_op =
291 block_op->getNonconstBlock(0, 0),
true);
292 RCP<Thyra::LinearOpBase<Scalar> > adjoint_op =
294 MEB::OutArgs<Scalar> adj_me_outArgs = adjoint_model_->createOutArgs();
295 adj_me_outArgs.set_W_op(adjoint_op);
296 adjoint_model_->evalModel(me_inArgs, adj_me_outArgs);
299 RCP<Thyra::ScaledIdentityLinearOpWithSolve<Scalar> > si_op =
301 block_op->getNonconstBlock(1, 1),
true);
310 RCP<Thyra::VectorBase<Scalar> > f = outArgs.
get_f();
311 if (f != Teuchos::null) {
312 RCP<const Thyra::VectorBase<Scalar> > x = inArgs.
get_x().assert_not_null();
313 RCP<const DPV> prod_x = rcp_dynamic_cast<
const DPV>(x,
true);
314 RCP<const Thyra::VectorBase<Scalar> > adjoint_x = prod_x->
getVectorBlock(0);
315 RCP<const Thyra::MultiVectorBase<Scalar> > adjoint_x_mv =
316 rcp_dynamic_cast<
const DMVPV>(adjoint_x,
true)->getMultiVector();
318 RCP<DPV> prod_f = rcp_dynamic_cast<
DPV>(f,
true);
319 RCP<Thyra::VectorBase<Scalar> > adjoint_f =
321 RCP<Thyra::MultiVectorBase<Scalar> > adjoint_f_mv =
322 rcp_dynamic_cast<DMVPV>(adjoint_f,
true)->getNonconstMultiVector();
324 MEB::OutArgs<Scalar> adj_me_outArgs = adjoint_model_->createOutArgs();
326 if (my_dfdx_ == Teuchos::null) my_dfdx_ = adjoint_model_->create_W_op();
327 adj_me_outArgs.set_W_op(my_dfdx_);
328 if (me_inArgs.supports(MEB::IN_ARG_alpha)) me_inArgs.set_alpha(0.0);
329 if (me_inArgs.supports(MEB::IN_ARG_beta)) me_inArgs.set_beta(1.0);
330 adjoint_model_->evalModel(me_inArgs, adj_me_outArgs);
333 my_dfdx_->apply(
Thyra::NOTRANS, *adjoint_x_mv, adjoint_f_mv.ptr(),
334 Scalar(-1.0), Scalar(0.0));
338 RCP<const DPV> prod_x_dot;
339 if (me_inArgs.supports(MEB::IN_ARG_x_dot)) {
340 RCP<const Thyra::VectorBase<Scalar> > x_dot = inArgs.
get_x_dot();
341 if (x_dot != Teuchos::null) {
342 prod_x_dot = rcp_dynamic_cast<
const DPV>(x_dot,
true);
343 RCP<const Thyra::VectorBase<Scalar> > adjoint_x_dot =
345 RCP<const Thyra::MultiVectorBase<Scalar> > adjoint_x_dot_mv =
346 rcp_dynamic_cast<
const DMVPV>(adjoint_x_dot,
true)
348 if (mass_matrix_is_identity_) {
350 Thyra::V_StVpV(adjoint_f_mv.ptr(), Scalar(-1.0), *adjoint_f_mv,
354 if (my_dfdxdot_ == Teuchos::null)
355 my_dfdxdot_ = adjoint_model_->create_W_op();
356 if (!mass_matrix_is_constant_ || !mass_matrix_is_computed_) {
357 adj_me_outArgs.set_W_op(my_dfdxdot_);
358 me_inArgs.set_alpha(1.0);
359 me_inArgs.set_beta(0.0);
360 adjoint_model_->evalModel(me_inArgs, adj_me_outArgs);
362 mass_matrix_is_computed_ =
true;
365 adjoint_f_mv.ptr(), Scalar(1.0), Scalar(-1.0));
372 RCP<Thyra::VectorBase<Scalar> > adjoint_g =
373 prod_f->getNonconstVectorBlock(1);
374 RCP<Thyra::MultiVectorBase<Scalar> > adjoint_g_mv =
375 rcp_dynamic_cast<DMVPV>(adjoint_g,
true)->getNonconstMultiVector();
377 MEB::OutArgs<Scalar> me_outArgs2 = model_->createOutArgs();
378 MEB::DerivativeSupport dfdp_support =
379 me_outArgs2.supports(MEB::OUT_ARG_DfDp, p_index_);
381 if (dfdp_support.supports(MEB::DERIV_LINEAR_OP)) {
382 if (my_dfdp_op_ == Teuchos::null)
383 my_dfdp_op_ = model_->create_DfDp_op(p_index_);
384 me_outArgs2.set_DfDp(p_index_, MEB::Derivative<Scalar>(my_dfdp_op_));
387 else if (dfdp_support.supports(MEB::DERIV_MV_JACOBIAN_FORM)) {
388 if (my_dfdp_mv_ == Teuchos::null)
389 my_dfdp_mv_ = createMembers(model_->get_f_space(),
390 model_->get_p_space(p_index_)->dim());
391 me_outArgs2.set_DfDp(
393 MEB::Derivative<Scalar>(my_dfdp_mv_, MEB::DERIV_MV_JACOBIAN_FORM));
394 my_dfdp_op_ = my_dfdp_mv_;
397 else if (dfdp_support.supports(MEB::DERIV_MV_GRADIENT_FORM)) {
398 if (my_dfdp_mv_ == Teuchos::null)
399 my_dfdp_mv_ = createMembers(model_->get_p_space(p_index_),
400 model_->get_f_space()->dim());
401 me_outArgs2.set_DfDp(
403 MEB::Derivative<Scalar>(my_dfdp_mv_, MEB::DERIV_MV_GRADIENT_FORM));
404 my_dfdp_op_ = my_dfdp_mv_;
409 "Invalid df/dp support");
410 model_->evalModel(me_inArgs, me_outArgs2);
411 my_dfdp_op_->apply(trans, *adjoint_x_mv, adjoint_g_mv.ptr(), Scalar(1.0),
414 if (prod_x_dot != Teuchos::null) {
415 RCP<const Thyra::VectorBase<Scalar> > z_dot =
416 prod_x_dot->getVectorBlock(1);
417 RCP<const Thyra::MultiVectorBase<Scalar> > z_dot_mv =
418 rcp_dynamic_cast<
const DMVPV>(z_dot,
true)->getMultiVector();
419 Thyra::V_VmV(adjoint_g_mv.ptr(), *z_dot_mv, *adjoint_g_mv);
424 template <
class Scalar>
429 pl->
set<
int>(
"Sensitivity Parameter Index", 0);
430 pl->
set<
int>(
"Response Function Index", 0);
431 pl->
set<
bool>(
"Mass Matrix Is Constant",
true);
432 pl->
set<
bool>(
"Mass Matrix Is Identity",
false);
Teuchos::RCP< const Thyra::LinearOpWithSolveFactoryBase< Scalar > > get_W_factory() const
Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > adjoint_model_
Implicit concrete LinearOpBase subclass that takes a flattended out multi-vector and performs a multi...
void setScale(const Scalar &s)
RCP< const VectorBase< Scalar > > get_x_dot() const
void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
static Teuchos::RCP< const Teuchos::ParameterList > getValidParameters()
T & get(const std::string &name, T def_value)
RCP< const VectorBase< Scalar > > getVectorBlock(const int k) const
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< Thyra::LinearOpBase< Scalar > > create_W_op() const
bool mass_matrix_is_constant_
RCP< LinearOpBase< Scalar > > getNonconstLinearOp()
Evaluation< VectorBase< Scalar > > get_f() const
Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
static magnitudeType rmax()
bool mass_matrix_is_identity_
Teuchos::RCP< const DPVS > f_prod_space_
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int p) const
Thyra::ModelEvaluatorBase::InArgs< Scalar > getNominalValues() const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Thyra::ModelEvaluatorBase::OutArgs< Scalar > createOutArgsImpl() const
virtual std::string description() const
void validateParametersAndSetDefaults(ParameterList const &validParamList, int const depth=1000)
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_f_space() const
Thyra::ModelEvaluatorBase::InArgs< Scalar > prototypeInArgs_
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
void setForwardSolutionHistory(const Teuchos::RCP< const Tempus::SolutionHistory< Scalar > > &sh)
Set solution history from forward evaluation.
RCP< VectorBase< Scalar > > getNonconstVectorBlock(const int k)
Implicit concrete LinearOpBase subclass that takes a flattended out multi-vector and performs a multi...
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_p_space(int p) const
Teuchos::RCP< const DPVS > x_prod_space_
Teuchos::RCP< const DMVPVS > residual_space_
AdjointAuxSensitivityModelEvaluator(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &adjoint_model, const Scalar &t_init, const Scalar &t_final, const Teuchos::RCP< const Teuchos::ParameterList > &pList=Teuchos::null)
Constructor.
Teuchos::RCP< const DMVPVS > response_space_
Thyra::ModelEvaluatorBase::OutArgs< Scalar > prototypeOutArgs_
Teuchos::RCP< const DMVPVS > adjoint_space_
#define TEUCHOS_ASSERT(assertion_test)
Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > model_
RCP< LinearOpBase< Scalar > > get_W_op() const
RCP< const VectorBase< Scalar > > get_x() const
void setFinalTime(const Scalar t_final)
Set the final time from the forward evaluation.
RCP< const VectorBase< Scalar > > get_p(int l) const