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>
29 const Scalar& t_final,
32 adjoint_model_(adjoint_model),
35 mass_matrix_is_computed_(false),
36 t_interp_(Teuchos::ScalarTraits<Scalar>::rmax())
46 if (pList != Teuchos::null)
51 p_index_ = pl->
get<
int>(
"Sensitivity Parameter Index", 0);
52 g_index_ = pl->
get<
int>(
"Response Function Index", 0);
58 "AdjointAuxSensitivityModelEvaluator currently does not support " <<
59 "non-constant mass matrix df/dx_dot!");
66 Thyra::multiVectorProductVectorSpace(
model_->get_p_space(p_index_),
68 Array< RCP<const VectorSpaceBase<Scalar> > > x_spaces(2);
69 Array< RCP<const VectorSpaceBase<Scalar> > > f_spaces(2);
78 MEB::InArgs<Scalar> me_inArgs =
model_->createInArgs();
81 MEB::InArgsSetup<Scalar> inArgs;
82 inArgs.setModelEvalDescription(this->
description());
83 inArgs.setSupports(MEB::IN_ARG_x);
84 inArgs.setSupports(MEB::IN_ARG_t);
85 if (me_inArgs.supports(MEB::IN_ARG_x_dot))
86 inArgs.setSupports(MEB::IN_ARG_x_dot);
87 inArgs.setSupports(MEB::IN_ARG_alpha);
88 inArgs.setSupports(MEB::IN_ARG_beta);
91 inArgs.set_Np(me_inArgs.Np());
94 MEB::OutArgs<Scalar> me_outArgs =
model_->createOutArgs();
95 MEB::OutArgs<Scalar> adj_me_outArgs =
adjoint_model_->createOutArgs();
96 MEB::OutArgsSetup<Scalar> outArgs;
97 outArgs.setModelEvalDescription(this->
description());
98 outArgs.set_Np_Ng(me_inArgs.Np(),0);
99 outArgs.setSupports(MEB::OUT_ARG_f);
100 outArgs.setSupports(MEB::OUT_ARG_W_op);
107 if (me_inArgs.supports(MEB::IN_ARG_x_dot)) {
113 template <
typename Scalar>
121 template <
typename Scalar>
129 forward_state_ = Teuchos::null;
132 template <
typename Scalar>
138 return model_->get_p_space(p);
141 template <
typename Scalar>
147 return model_->get_p_names(p);
150 template <
typename Scalar>
155 return x_prod_space_;
158 template <
typename Scalar>
163 return f_prod_space_;
166 template <
typename Scalar>
174 RCP<LinearOpBase<Scalar> > adjoint_op = adjoint_model_->create_W_op();
175 RCP<LinearOpBase<Scalar> > mv_adjoint_op =
176 Thyra::nonconstMultiVectorLinearOp(adjoint_op, num_adjoint_);
177 RCP<const Thyra::VectorSpaceBase<Scalar> > g_space = response_space_;
178 RCP<LinearOpBase<Scalar> > g_op =
179 Thyra::scaledIdentity(g_space, Scalar(1.0));
180 RCP<LinearOpBase<Scalar> > null_op;
181 return nonconstBlock2x2(mv_adjoint_op, null_op, null_op, g_op);
184 template <
typename Scalar>
190 using Teuchos::rcp_dynamic_cast;
193 RCP<const LOWSFB > alowsfb = adjoint_model_->get_W_factory();
194 if (alowsfb == Teuchos::null)
195 return Teuchos::null;
197 RCP<const LOWSFB > mv_alowsfb =
198 Thyra::multiVectorLinearOpWithSolveFactory(alowsfb, residual_space_,
201 RCP<const Thyra::VectorSpaceBase<Scalar> > g_space = response_space_;
202 RCP<const LOWSFB > g_lowsfb =
203 Thyra::scaledIdentitySolveFactory(g_space, Scalar(1.0));
206 lowsfbs[0] = mv_alowsfb;
207 lowsfbs[1] = g_lowsfb;
208 return Thyra::blockedTriangularLinearOpWithSolveFactory(lowsfbs);
211 template <
typename Scalar>
216 return prototypeInArgs_;
219 template <
typename Scalar>
226 using Teuchos::rcp_dynamic_cast;
228 MEB::InArgs<Scalar> me_nominal = model_->getNominalValues();
229 MEB::InArgs<Scalar> nominal = this->createInArgs();
234 RCP< Thyra::VectorBase<Scalar> > x = Thyra::createMember(*x_prod_space_);
235 Thyra::assign(x.ptr(), zero);
238 if (me_nominal.supports(MEB::IN_ARG_x_dot)) {
239 RCP< Thyra::VectorBase<Scalar> > x_dot =
240 Thyra::createMember(*x_prod_space_);
241 Thyra::assign(x_dot.ptr(), zero);
242 nominal.set_x_dot(x_dot);
245 const int np = model_->Np();
246 for (
int i=0; i<np; ++i)
247 nominal.set_p(i, me_nominal.get_p(i));
252 template <
typename Scalar>
257 return prototypeOutArgs_;
260 template <
typename Scalar>
268 using Teuchos::rcp_dynamic_cast;
278 const Scalar t = inArgs.
get_t();
279 const Scalar forward_t = t_final_ + t_init_ - t;
280 if (forward_state_ == Teuchos::null || t_interp_ != t) {
281 if (forward_state_ == Teuchos::null)
282 forward_state_ = sh_->interpolateState(forward_t);
284 sh_->interpolateState(forward_t, forward_state_.get());
289 MEB::InArgs<Scalar> me_inArgs = model_->getNominalValues();
290 me_inArgs.set_x(forward_state_->getX());
291 if (me_inArgs.supports(MEB::IN_ARG_x_dot) &&
293 me_inArgs.set_x_dot(forward_state_->getXDot());
294 if (me_inArgs.supports(MEB::IN_ARG_t))
295 me_inArgs.set_t(forward_t);
296 const int np = me_inArgs.Np();
297 for (
int i=0; i<np; ++i)
298 me_inArgs.set_p(i, inArgs.
get_p(i));
301 RCP<Thyra::LinearOpBase<Scalar> > op = outArgs.
get_W_op();
302 if (op != Teuchos::null) {
303 if (me_inArgs.supports(MEB::IN_ARG_alpha))
305 if (me_inArgs.supports(MEB::IN_ARG_beta))
306 me_inArgs.set_beta(inArgs.
get_beta());
309 RCP<Thyra::DefaultBlockedLinearOp<Scalar> > block_op =
311 RCP<Thyra::MultiVectorLinearOp<Scalar> > mv_adjoint_op =
313 block_op->getNonconstBlock(0,0),
true);
314 RCP<Thyra::LinearOpBase<Scalar> > adjoint_op =
316 MEB::OutArgs<Scalar> adj_me_outArgs = adjoint_model_->createOutArgs();
317 adj_me_outArgs.set_W_op(adjoint_op);
318 adjoint_model_->evalModel(me_inArgs, adj_me_outArgs);
321 RCP<Thyra::ScaledIdentityLinearOpWithSolve<Scalar> > si_op =
323 block_op->getNonconstBlock(1,1),
true);
332 RCP<Thyra::VectorBase<Scalar> > f = outArgs.
get_f();
333 if (f != Teuchos::null) {
334 RCP<const Thyra::VectorBase<Scalar> > x = inArgs.
get_x().assert_not_null();
335 RCP<const DPV> prod_x = rcp_dynamic_cast<
const DPV>(x,
true);
336 RCP<const Thyra::VectorBase<Scalar> > adjoint_x = prod_x->
getVectorBlock(0);
337 RCP<const Thyra::MultiVectorBase<Scalar> >adjoint_x_mv =
338 rcp_dynamic_cast<
const DMVPV>(adjoint_x,
true)->getMultiVector();
340 RCP<DPV> prod_f = rcp_dynamic_cast<
DPV>(f,
true);
341 RCP<Thyra::VectorBase<Scalar> > adjoint_f =
343 RCP<Thyra::MultiVectorBase<Scalar> > adjoint_f_mv =
344 rcp_dynamic_cast<DMVPV>(adjoint_f,
true)->getNonconstMultiVector();
346 MEB::OutArgs<Scalar> adj_me_outArgs = adjoint_model_->createOutArgs();
348 if (my_dfdx_ == Teuchos::null)
349 my_dfdx_ = adjoint_model_->create_W_op();
350 adj_me_outArgs.set_W_op(my_dfdx_);
351 if (me_inArgs.supports(MEB::IN_ARG_alpha))
352 me_inArgs.set_alpha(0.0);
353 if (me_inArgs.supports(MEB::IN_ARG_beta))
354 me_inArgs.set_beta(1.0);
355 adjoint_model_->evalModel(me_inArgs, adj_me_outArgs);
358 my_dfdx_->apply(
Thyra::NOTRANS, *adjoint_x_mv, adjoint_f_mv.ptr(),
359 Scalar(-1.0), Scalar(0.0));
363 RCP<const DPV> prod_x_dot;
364 if (me_inArgs.supports(MEB::IN_ARG_x_dot)) {
365 RCP<const Thyra::VectorBase<Scalar> > x_dot = inArgs.
get_x_dot();
366 if (x_dot != Teuchos::null) {
367 prod_x_dot = rcp_dynamic_cast<
const DPV>(x_dot,
true);
368 RCP<const Thyra::VectorBase<Scalar> > adjoint_x_dot =
370 RCP<const Thyra::MultiVectorBase<Scalar> > adjoint_x_dot_mv =
371 rcp_dynamic_cast<
const DMVPV>(adjoint_x_dot,
true)->getMultiVector();
372 if (mass_matrix_is_identity_) {
374 Thyra::V_StVpV(adjoint_f_mv.ptr(), Scalar(-1.0), *adjoint_f_mv,
378 if (my_dfdxdot_ == Teuchos::null)
379 my_dfdxdot_ = adjoint_model_->create_W_op();
380 if (!mass_matrix_is_constant_ || !mass_matrix_is_computed_) {
381 adj_me_outArgs.set_W_op(my_dfdxdot_);
382 me_inArgs.set_alpha(1.0);
383 me_inArgs.set_beta(0.0);
384 adjoint_model_->evalModel(me_inArgs, adj_me_outArgs);
386 mass_matrix_is_computed_ =
true;
389 adjoint_f_mv.ptr(), Scalar(1.0), Scalar(-1.0));
396 RCP<Thyra::VectorBase<Scalar> > adjoint_g =
397 prod_f->getNonconstVectorBlock(1);
398 RCP<Thyra::MultiVectorBase<Scalar> > adjoint_g_mv =
399 rcp_dynamic_cast<DMVPV>(adjoint_g,
true)->getNonconstMultiVector();
401 MEB::OutArgs<Scalar> me_outArgs2 = model_->createOutArgs();
402 MEB::DerivativeSupport dfdp_support =
403 me_outArgs2.supports(MEB::OUT_ARG_DfDp, p_index_);
405 if (dfdp_support.supports(MEB::DERIV_LINEAR_OP)) {
406 if (my_dfdp_op_ == Teuchos::null)
407 my_dfdp_op_ = model_->create_DfDp_op(p_index_);
408 me_outArgs2.set_DfDp(p_index_, MEB::Derivative<Scalar>(my_dfdp_op_));
411 else if (dfdp_support.supports(MEB::DERIV_MV_JACOBIAN_FORM)) {
412 if (my_dfdp_mv_ == Teuchos::null)
413 my_dfdp_mv_ = createMembers(model_->get_f_space(),
414 model_->get_p_space(p_index_)->dim());
415 me_outArgs2.set_DfDp(p_index_,
416 MEB::Derivative<Scalar>(my_dfdp_mv_,
417 MEB::DERIV_MV_JACOBIAN_FORM));
418 my_dfdp_op_ = my_dfdp_mv_;
421 else if (dfdp_support.supports(MEB::DERIV_MV_GRADIENT_FORM)) {
422 if (my_dfdp_mv_ == Teuchos::null)
423 my_dfdp_mv_ = createMembers(model_->get_p_space(p_index_),
424 model_->get_f_space()->dim());
425 me_outArgs2.set_DfDp(p_index_,
426 MEB::Derivative<Scalar>(my_dfdp_mv_,
427 MEB::DERIV_MV_GRADIENT_FORM));
428 my_dfdp_op_ = my_dfdp_mv_;
433 true, std::logic_error,
"Invalid df/dp support");
434 model_->evalModel(me_inArgs, me_outArgs2);
435 my_dfdp_op_->apply(trans, *adjoint_x_mv, adjoint_g_mv.ptr(),
436 Scalar(1.0), Scalar(0.0));
438 if (prod_x_dot != Teuchos::null) {
439 RCP<const Thyra::VectorBase<Scalar> > z_dot =
440 prod_x_dot->getVectorBlock(1);
441 RCP<const Thyra::MultiVectorBase<Scalar> > z_dot_mv =
442 rcp_dynamic_cast<
const DMVPV>(z_dot,
true)->getMultiVector();
443 Thyra::V_VmV(adjoint_g_mv.ptr(), *z_dot_mv, *adjoint_g_mv);
448 template<
class Scalar>
454 pl->
set<
int>(
"Sensitivity Parameter Index", 0);
455 pl->
set<
int>(
"Response Function Index", 0);
456 pl->
set<
bool>(
"Mass Matrix Is Constant",
true);
457 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