9 #ifndef Tempus_AdjointSensitivityModelEvaluator_impl_hpp
10 #define Tempus_AdjointSensitivityModelEvaluator_impl_hpp
13 #include "Tempus_config.hpp"
15 #include "Thyra_DefaultMultiVectorProductVectorSpace.hpp"
16 #include "Thyra_DefaultMultiVectorProductVector.hpp"
19 #include "Thyra_DefaultScaledAdjointLinearOp.hpp"
20 #include "Thyra_DefaultAdjointLinearOpWithSolve.hpp"
22 #include "Thyra_VectorStdOps.hpp"
23 #include "Thyra_MultiVectorStdOps.hpp"
27 template <
typename Scalar>
31 adjoint_residual_model,
34 const Scalar& t_init,
const Scalar& t_final,
const bool is_pseudotransient,
37 adjoint_residual_model_(adjoint_residual_model),
38 adjoint_solve_model_(adjoint_solve_model),
41 is_pseudotransient_(is_pseudotransient),
42 mass_matrix_is_computed_(false),
43 jacobian_matrix_is_computed_(false),
44 response_gradient_is_computed_(false),
45 t_interp_(Teuchos::ScalarTraits<Scalar>::rmax())
52 if (pList != Teuchos::null) *pl = *pList;
56 p_index_ = pl->
get<
int>(
"Sensitivity Parameter Index", 0);
57 g_index_ = pl->
get<
int>(
"Response Function Index", 0);
63 "AdjointSensitivityModelEvaluator currently does not support "
64 <<
"non-constant mass matrix df/dx_dot!");
74 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 if (me_inArgs.supports(MEB::IN_ARG_alpha))
85 inArgs.setSupports(MEB::IN_ARG_alpha);
86 if (me_inArgs.supports(MEB::IN_ARG_beta))
87 inArgs.setSupports(MEB::IN_ARG_beta);
90 inArgs.set_Np(me_inArgs.Np());
93 MEB::OutArgs<Scalar> me_outArgs =
model_->createOutArgs();
94 MEB::OutArgs<Scalar> adj_mer_outArgs =
97 MEB::OutArgsSetup<Scalar> outArgs;
98 outArgs.setModelEvalDescription(this->
description());
99 outArgs.set_Np_Ng(me_inArgs.Np(), 2);
100 outArgs.setSupports(MEB::OUT_ARG_f);
101 if (adj_mes_outArgs.supports(MEB::OUT_ARG_W_op))
102 outArgs.setSupports(MEB::OUT_ARG_W_op);
109 if (me_inArgs.supports(MEB::IN_ARG_x_dot)) {
113 MEB::DerivativeSupport dgdx_support =
114 me_outArgs.supports(MEB::OUT_ARG_DgDx, g_index_);
115 MEB::DerivativeSupport dgdp_support =
116 me_outArgs.supports(MEB::OUT_ARG_DgDp, g_index_, p_index_);
117 TEUCHOS_ASSERT(dgdx_support.supports(MEB::DERIV_MV_GRADIENT_FORM));
118 TEUCHOS_ASSERT(dgdp_support.supports(MEB::DERIV_MV_GRADIENT_FORM));
121 template <
typename Scalar>
123 const Scalar t_final)
128 template <
typename Scalar>
133 if (is_pseudotransient_)
134 forward_state_ = sh_->getCurrentState();
137 forward_state_ = Teuchos::null;
141 mass_matrix_is_computed_ =
false;
142 jacobian_matrix_is_computed_ =
false;
143 response_gradient_is_computed_ =
false;
146 template <
typename Scalar>
151 return model_->get_p_space(p);
154 template <
typename Scalar>
159 return model_->get_p_names(p);
162 template <
typename Scalar>
166 return adjoint_space_;
169 template <
typename Scalar>
173 return residual_space_;
176 template <
typename Scalar>
181 if (j == 0)
return response_space_;
182 return model_->get_g_space(g_index_);
185 template <
typename Scalar>
190 adjoint_solve_model_->create_W_op();
191 if (adjoint_op == Teuchos::null)
return Teuchos::null;
193 return Thyra::nonconstMultiVectorLinearOp(adjoint_op, num_adjoint_);
196 template <
typename Scalar>
201 using Teuchos::rcp_dynamic_cast;
204 RCP<const LOWSFB> alowsfb = adjoint_solve_model_->get_W_factory();
205 if (alowsfb == Teuchos::null)
206 return Teuchos::null;
208 return Thyra::multiVectorLinearOpWithSolveFactory(alowsfb, residual_space_,
212 template <
typename Scalar>
216 return prototypeInArgs_;
219 template <
typename Scalar>
225 using Teuchos::rcp_dynamic_cast;
227 MEB::InArgs<Scalar> me_nominal = model_->getNominalValues();
228 MEB::InArgs<Scalar> nominal = this->createInArgs();
233 RCP<Thyra::VectorBase<Scalar> > x = Thyra::createMember(*adjoint_space_);
234 Thyra::assign(x.ptr(), zero);
237 if (me_nominal.supports(MEB::IN_ARG_x_dot)) {
238 RCP<Thyra::VectorBase<Scalar> > x_dot =
239 Thyra::createMember(*adjoint_space_);
240 Thyra::assign(x_dot.ptr(), zero);
241 nominal.set_x_dot(x_dot);
244 const int np = model_->Np();
245 for (
int i = 0; i < np; ++i) nominal.set_p(i, me_nominal.get_p(i));
250 template <
typename Scalar>
254 return prototypeOutArgs_;
257 template <
typename Scalar>
264 using Teuchos::rcp_dynamic_cast;
266 TEMPUS_FUNC_TIME_MONITOR_DIFF(
267 "Tempus::AdjointSensitivityModelEvaluator::evalModel()", TEMPUS_EVAL);
277 const Scalar t = inArgs.
get_t();
279 if (is_pseudotransient_)
280 forward_t = forward_state_->getTime();
282 forward_t = t_final_ + t_init_ - t;
283 if (forward_state_ == Teuchos::null || t_interp_ != t) {
284 if (forward_state_ == Teuchos::null)
285 forward_state_ = sh_->interpolateState(forward_t);
287 sh_->interpolateState(forward_t, forward_state_.get());
293 MEB::InArgs<Scalar> me_inArgs = model_->getNominalValues();
294 me_inArgs.set_x(forward_state_->getX());
295 if (me_inArgs.supports(MEB::IN_ARG_x_dot)) {
297 me_inArgs.set_x_dot(forward_state_->getXDot());
299 if (is_pseudotransient_) {
304 if (my_x_dot_ == Teuchos::null) {
305 my_x_dot_ = Thyra::createMember(model_->get_x_space());
306 Thyra::assign(my_x_dot_.ptr(), Scalar(0.0));
308 me_inArgs.set_x_dot(my_x_dot_);
313 me_inArgs.set_x_dot(Teuchos::null);
317 if (me_inArgs.supports(MEB::IN_ARG_t)) me_inArgs.set_t(forward_t);
318 const int np = me_inArgs.Np();
319 for (
int i = 0; i < np; ++i) me_inArgs.set_p(i, inArgs.
get_p(i));
325 RCP<Thyra::LinearOpBase<Scalar> > op;
327 if (op != Teuchos::null) {
328 TEMPUS_FUNC_TIME_MONITOR_DIFF(
329 "Tempus::AdjointSensitivityModelEvaluator::evalModel::W",
331 if (me_inArgs.supports(MEB::IN_ARG_alpha))
333 if (me_inArgs.supports(MEB::IN_ARG_beta)) {
334 if (me_inArgs.supports(MEB::IN_ARG_x_dot))
335 me_inArgs.set_beta(inArgs.
get_beta());
337 me_inArgs.set_beta(-inArgs.
get_beta());
340 RCP<Thyra::MultiVectorLinearOp<Scalar> > mv_adjoint_op =
342 RCP<Thyra::LinearOpBase<Scalar> > adjoint_op =
344 MEB::OutArgs<Scalar> adj_me_outArgs = adjoint_solve_model_->createOutArgs();
345 adj_me_outArgs.set_W_op(adjoint_op);
346 adjoint_solve_model_->evalModel(me_inArgs, adj_me_outArgs);
349 RCP<Thyra::VectorBase<Scalar> > adjoint_f = outArgs.
get_f();
350 RCP<Thyra::VectorBase<Scalar> > adjoint_g = outArgs.
get_g(0);
351 RCP<Thyra::VectorBase<Scalar> > g = outArgs.
get_g(1);
352 RCP<const Thyra::MultiVectorBase<Scalar> > adjoint_x_mv;
353 RCP<const Thyra::VectorBase<Scalar> > adjoint_x;
354 if (adjoint_f != Teuchos::null || adjoint_g != Teuchos::null) {
355 adjoint_x = inArgs.
get_x().assert_not_null();
357 rcp_dynamic_cast<
const DMVPV>(adjoint_x,
true)->getMultiVector();
365 if (adjoint_f != Teuchos::null) {
366 TEMPUS_FUNC_TIME_MONITOR_DIFF(
367 "Tempus::AdjointSensitivityModelEvaluator::evalModel::f",
370 RCP<Thyra::MultiVectorBase<Scalar> > adjoint_f_mv =
371 rcp_dynamic_cast<
DMVPV>(adjoint_f,
true)->getNonconstMultiVector();
373 MEB::OutArgs<Scalar> me_outArgs = model_->createOutArgs();
374 MEB::OutArgs<Scalar> adj_me_outArgs =
375 adjoint_residual_model_->createOutArgs();
380 TEMPUS_FUNC_TIME_MONITOR_DIFF(
381 "Tempus::AdjointSensitivityModelEvaluator::evalModel::dg/dx",
383 if (my_dgdx_mv_ == Teuchos::null)
384 my_dgdx_mv_ = Thyra::createMembers(
385 model_->get_x_space(), model_->get_g_space(g_index_)->dim());
386 if (!response_gradient_is_computed_) {
389 MEB::Derivative<Scalar>(my_dgdx_mv_, MEB::DERIV_MV_GRADIENT_FORM));
390 model_->evalModel(me_inArgs, me_outArgs);
391 me_outArgs.set_DgDx(g_index_, MEB::Derivative<Scalar>());
392 if (is_pseudotransient_) response_gradient_is_computed_ =
true;
394 Thyra::assign(adjoint_f_mv.ptr(), *my_dgdx_mv_);
400 TEMPUS_FUNC_TIME_MONITOR_DIFF(
401 "Tempus::AdjointSensitivityModelEvaluator::evalModel::df/dx",
403 if (my_dfdx_ == Teuchos::null)
404 my_dfdx_ = adjoint_residual_model_->create_W_op();
405 if (!jacobian_matrix_is_computed_) {
406 adj_me_outArgs.set_W_op(my_dfdx_);
407 if (me_inArgs.supports(MEB::IN_ARG_alpha)) me_inArgs.set_alpha(0.0);
408 if (me_inArgs.supports(MEB::IN_ARG_beta)) me_inArgs.set_beta(1.0);
409 adjoint_residual_model_->evalModel(me_inArgs, adj_me_outArgs);
410 if (is_pseudotransient_) jacobian_matrix_is_computed_ =
true;
412 my_dfdx_->apply(
Thyra::NOTRANS, *adjoint_x_mv, adjoint_f_mv.ptr(),
413 Scalar(-1.0), Scalar(1.0));
420 if (me_inArgs.supports(MEB::IN_ARG_x_dot)) {
421 TEMPUS_FUNC_TIME_MONITOR_DIFF(
422 "Tempus::AdjointSensitivityModelEvaluator::evalModel::df/dx_dot",
423 TEMPUS_EVAL_DFDXDOT);
424 RCP<const Thyra::VectorBase<Scalar> > adjoint_x_dot = inArgs.
get_x_dot();
425 if (adjoint_x_dot != Teuchos::null) {
426 RCP<const Thyra::MultiVectorBase<Scalar> > adjoint_x_dot_mv =
427 rcp_dynamic_cast<
const DMVPV>(adjoint_x_dot,
true)
429 if (mass_matrix_is_identity_) {
431 Thyra::V_StVpV(adjoint_f_mv.ptr(), Scalar(-1.0), *adjoint_f_mv,
435 if (my_dfdxdot_ == Teuchos::null)
436 my_dfdxdot_ = adjoint_residual_model_->create_W_op();
437 if (!mass_matrix_is_computed_) {
438 adj_me_outArgs.set_W_op(my_dfdxdot_);
439 me_inArgs.set_alpha(1.0);
440 me_inArgs.set_beta(0.0);
441 adjoint_residual_model_->evalModel(me_inArgs, adj_me_outArgs);
442 if (is_pseudotransient_ || mass_matrix_is_constant_)
443 mass_matrix_is_computed_ =
true;
446 adjoint_f_mv.ptr(), Scalar(1.0), Scalar(-1.0));
456 if (adjoint_g != Teuchos::null) {
457 TEMPUS_FUNC_TIME_MONITOR_DIFF(
458 "Tempus::AdjointSensitivityModelEvaluator::evalModel::g",
460 RCP<Thyra::MultiVectorBase<Scalar> > adjoint_g_mv =
461 rcp_dynamic_cast<
DMVPV>(adjoint_g,
true)->getNonconstMultiVector();
463 MEB::OutArgs<Scalar> me_outArgs = model_->createOutArgs();
466 MEB::DerivativeSupport dgdp_support =
467 me_outArgs.supports(MEB::OUT_ARG_DgDp, g_index_, p_index_);
468 if (dgdp_support.supports(MEB::DERIV_MV_GRADIENT_FORM)) {
471 MEB::Derivative<Scalar>(adjoint_g_mv, MEB::DERIV_MV_GRADIENT_FORM));
472 model_->evalModel(me_inArgs, me_outArgs);
474 else if (dgdp_support.supports(MEB::DERIV_MV_JACOBIAN_FORM)) {
475 const int num_g = model_->get_g_space(g_index_)->dim();
476 const int num_p = model_->get_p_space(p_index_)->dim();
477 RCP<Thyra::MultiVectorBase<Scalar> > dgdp_trans =
478 createMembers(model_->get_g_space(g_index_), num_p);
481 MEB::Derivative<Scalar>(dgdp_trans, MEB::DERIV_MV_JACOBIAN_FORM));
482 model_->evalModel(me_inArgs, me_outArgs);
485 for (
int i = 0; i < num_p; ++i)
486 for (
int j = 0; j < num_g; ++j) dgdp_view(i, j) = dgdp_trans_view(j, i);
490 "Invalid dg/dp support");
491 me_outArgs.set_DgDp(g_index_, p_index_, MEB::Derivative<Scalar>());
494 MEB::DerivativeSupport dfdp_support =
495 me_outArgs.supports(MEB::OUT_ARG_DfDp, p_index_);
497 if (dfdp_support.supports(MEB::DERIV_LINEAR_OP)) {
498 if (my_dfdp_op_ == Teuchos::null)
499 my_dfdp_op_ = model_->create_DfDp_op(p_index_);
500 me_outArgs.set_DfDp(p_index_, MEB::Derivative<Scalar>(my_dfdp_op_));
503 else if (dfdp_support.supports(MEB::DERIV_MV_JACOBIAN_FORM)) {
504 if (my_dfdp_mv_ == Teuchos::null)
505 my_dfdp_mv_ = createMembers(model_->get_f_space(),
506 model_->get_p_space(p_index_)->dim());
509 MEB::Derivative<Scalar>(my_dfdp_mv_, MEB::DERIV_MV_JACOBIAN_FORM));
510 my_dfdp_op_ = my_dfdp_mv_;
513 else if (dfdp_support.supports(MEB::DERIV_MV_GRADIENT_FORM)) {
514 if (my_dfdp_mv_ == Teuchos::null)
515 my_dfdp_mv_ = createMembers(model_->get_p_space(p_index_),
516 model_->get_f_space()->dim());
519 MEB::Derivative<Scalar>(my_dfdp_mv_, MEB::DERIV_MV_GRADIENT_FORM));
520 my_dfdp_op_ = my_dfdp_mv_;
525 "Invalid df/dp support");
526 model_->evalModel(me_inArgs, me_outArgs);
527 my_dfdp_op_->apply(trans, *adjoint_x_mv, adjoint_g_mv.ptr(), Scalar(-1.0),
531 if (g != Teuchos::null) {
532 MEB::OutArgs<Scalar> me_outArgs = model_->createOutArgs();
533 me_outArgs.set_g(g_index_, g);
534 model_->evalModel(me_inArgs, me_outArgs);
538 template <
class Scalar>
543 pl->
set<
int>(
"Sensitivity Parameter Index", 0);
544 pl->
set<
int>(
"Response Function Index", 0);
545 pl->
set<
bool>(
"Mass Matrix Is Constant",
true);
546 pl->
set<
bool>(
"Mass Matrix Is Identity",
false);
static Teuchos::RCP< const Teuchos::ParameterList > getValidParameters()
Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > adjoint_solve_model_
Teuchos::RCP< const DMVPVS > residual_space_
Teuchos::RCP< const DMVPVS > adjoint_space_
RCP< const VectorBase< Scalar > > get_x_dot() const
bool mass_matrix_is_identity_
T & get(const std::string &name, T def_value)
Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() 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< const Thyra::ModelEvaluator< Scalar > > adjoint_residual_model_
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
RCP< LinearOpBase< Scalar > > getNonconstLinearOp()
Evaluation< VectorBase< Scalar > > get_g(int j) const
Evaluation< VectorBase< Scalar > > get_f() const
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int p) const
Teuchos::RCP< const DMVPVS > response_space_
bool supports(EOutArgsMembers arg) const
Teuchos::RCP< const Thyra::LinearOpWithSolveFactoryBase< Scalar > > get_W_factory() const
static magnitudeType rmax()
bool mass_matrix_is_constant_
void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
Thyra::ModelEvaluatorBase::OutArgs< Scalar > createOutArgsImpl() const
Thyra::ModelEvaluatorBase::InArgs< Scalar > prototypeInArgs_
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Thyra::ModelEvaluatorBase::OutArgs< Scalar > prototypeOutArgs_
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_f_space() const
virtual std::string description() const
void validateParametersAndSetDefaults(ParameterList const &validParamList, int const depth=1000)
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_g_space(int j) const
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > model_
AdjointSensitivityModelEvaluator(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &adjoint_residual_model, const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &adjoint_solve_model, const Scalar &t_init, const Scalar &t_final, const bool is_pseudotransient, const Teuchos::RCP< const Teuchos::ParameterList > &pList=Teuchos::null)
Constructor.
void setForwardSolutionHistory(const Teuchos::RCP< const Tempus::SolutionHistory< Scalar > > &sh)
Set solution history from forward evaluation.
Implicit concrete LinearOpBase subclass that takes a flattended out multi-vector and performs a multi...
Thyra::ModelEvaluatorBase::InArgs< Scalar > getNominalValues() const
void setFinalTime(const Scalar t_final)
Set the final time from the forward evaluation.
#define TEUCHOS_ASSERT(assertion_test)
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_W_op() const
RCP< LinearOpBase< Scalar > > get_W_op() const
RCP< const VectorBase< Scalar > > get_x() const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_p_space(int p) const
RCP< const VectorBase< Scalar > > get_p(int l) const