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>
34 const Scalar& t_final,
35 const bool is_pseudotransient,
38 adjoint_residual_model_(adjoint_residual_model),
39 adjoint_solve_model_(adjoint_solve_model),
42 is_pseudotransient_(is_pseudotransient),
43 mass_matrix_is_computed_(false),
44 jacobian_matrix_is_computed_(false),
45 response_gradient_is_computed_(false),
46 t_interp_(Teuchos::ScalarTraits<Scalar>::rmax())
53 if (pList != Teuchos::null)
58 p_index_ = pl->
get<
int>(
"Sensitivity Parameter Index", 0);
59 g_index_ = pl->
get<
int>(
"Response Function Index", 0);
65 "AdjointSensitivityModelEvaluator currently does not support " <<
66 "non-constant mass matrix df/dx_dot!");
73 Thyra::multiVectorProductVectorSpace(
model_->get_p_space(p_index_),
77 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 if (me_inArgs.supports(MEB::IN_ARG_alpha))
88 inArgs.setSupports(MEB::IN_ARG_alpha);
89 if (me_inArgs.supports(MEB::IN_ARG_beta))
90 inArgs.setSupports(MEB::IN_ARG_beta);
93 inArgs.set_Np(me_inArgs.Np());
96 MEB::OutArgs<Scalar> me_outArgs =
model_->createOutArgs();
99 MEB::OutArgsSetup<Scalar> outArgs;
100 outArgs.setModelEvalDescription(this->
description());
101 outArgs.set_Np_Ng(me_inArgs.Np(),2);
102 outArgs.setSupports(MEB::OUT_ARG_f);
103 if (adj_mes_outArgs.supports(MEB::OUT_ARG_W_op))
104 outArgs.setSupports(MEB::OUT_ARG_W_op);
111 if (me_inArgs.supports(MEB::IN_ARG_x_dot)) {
115 MEB::DerivativeSupport dgdx_support =
116 me_outArgs.supports(MEB::OUT_ARG_DgDx, g_index_);
117 MEB::DerivativeSupport dgdp_support =
118 me_outArgs.supports(MEB::OUT_ARG_DgDp, g_index_, p_index_);
119 TEUCHOS_ASSERT(dgdx_support.supports(MEB::DERIV_MV_GRADIENT_FORM));
120 TEUCHOS_ASSERT(dgdp_support.supports(MEB::DERIV_MV_GRADIENT_FORM));
123 template <
typename Scalar>
131 template <
typename Scalar>
138 if (is_pseudotransient_)
139 forward_state_ = sh_->getCurrentState();
142 forward_state_ = Teuchos::null;
146 mass_matrix_is_computed_ =
false;
147 jacobian_matrix_is_computed_ =
false;
148 response_gradient_is_computed_ =
false;
151 template <
typename Scalar>
157 return model_->get_p_space(p);
160 template <
typename Scalar>
166 return model_->get_p_names(p);
169 template <
typename Scalar>
174 return adjoint_space_;
177 template <
typename Scalar>
182 return residual_space_;
185 template <
typename Scalar>
192 return response_space_;
193 return model_->get_g_space(g_index_);
196 template <
typename Scalar>
202 adjoint_solve_model_->create_W_op();
203 if (adjoint_op == Teuchos::null)
204 return Teuchos::null;
206 return Thyra::nonconstMultiVectorLinearOp(adjoint_op, num_adjoint_);
209 template <
typename Scalar>
215 using Teuchos::rcp_dynamic_cast;
218 RCP<const LOWSFB> alowsfb = adjoint_solve_model_->get_W_factory();
219 if (alowsfb == Teuchos::null)
220 return Teuchos::null;
222 return Thyra::multiVectorLinearOpWithSolveFactory(
223 alowsfb, residual_space_, adjoint_space_);
226 template <
typename Scalar>
231 return prototypeInArgs_;
234 template <
typename Scalar>
241 using Teuchos::rcp_dynamic_cast;
243 MEB::InArgs<Scalar> me_nominal = model_->getNominalValues();
244 MEB::InArgs<Scalar> nominal = this->createInArgs();
249 RCP< Thyra::VectorBase<Scalar> > x = Thyra::createMember(*adjoint_space_);
250 Thyra::assign(x.ptr(), zero);
253 if (me_nominal.supports(MEB::IN_ARG_x_dot)) {
254 RCP< Thyra::VectorBase<Scalar> > x_dot =
255 Thyra::createMember(*adjoint_space_);
256 Thyra::assign(x_dot.ptr(), zero);
257 nominal.set_x_dot(x_dot);
260 const int np = model_->Np();
261 for (
int i=0; i<np; ++i)
262 nominal.set_p(i, me_nominal.get_p(i));
267 template <
typename Scalar>
272 return prototypeOutArgs_;
275 template <
typename Scalar>
283 using Teuchos::rcp_dynamic_cast;
285 TEMPUS_FUNC_TIME_MONITOR_DIFF(
"Tempus::AdjointSensitivityModelEvaluator::evalModel()", TEMPUS_EVAL);
295 const Scalar t = inArgs.
get_t();
297 if (is_pseudotransient_)
298 forward_t = forward_state_->getTime();
300 forward_t = t_final_ + t_init_ - t;
301 if (forward_state_ == Teuchos::null || t_interp_ != t) {
302 if (forward_state_ == Teuchos::null)
303 forward_state_ = sh_->interpolateState(forward_t);
305 sh_->interpolateState(forward_t, forward_state_.get());
311 MEB::InArgs<Scalar> me_inArgs = model_->getNominalValues();
312 me_inArgs.set_x(forward_state_->getX());
313 if (me_inArgs.supports(MEB::IN_ARG_x_dot)) {
315 me_inArgs.set_x_dot(forward_state_->getXDot());
317 if (is_pseudotransient_) {
322 if (my_x_dot_ == Teuchos::null) {
323 my_x_dot_ = Thyra::createMember(model_->get_x_space());
324 Thyra::assign(my_x_dot_.ptr(), Scalar(0.0));
326 me_inArgs.set_x_dot(my_x_dot_);
331 me_inArgs.set_x_dot(Teuchos::null);
335 if (me_inArgs.supports(MEB::IN_ARG_t))
336 me_inArgs.set_t(forward_t);
337 const int np = me_inArgs.Np();
338 for (
int i=0; i<np; ++i)
339 me_inArgs.set_p(i, inArgs.
get_p(i));
345 RCP<Thyra::LinearOpBase<Scalar> > op;
346 if (outArgs.
supports(MEB::OUT_ARG_W_op))
348 if (op != Teuchos::null) {
349 TEMPUS_FUNC_TIME_MONITOR_DIFF(
"Tempus::AdjointSensitivityModelEvaluator::evalModel::W", TEMPUS_EVAL_W);
350 if (me_inArgs.supports(MEB::IN_ARG_alpha))
352 if (me_inArgs.supports(MEB::IN_ARG_beta)) {
353 if (me_inArgs.supports(MEB::IN_ARG_x_dot))
354 me_inArgs.set_beta(inArgs.
get_beta());
356 me_inArgs.set_beta(-inArgs.
get_beta());
359 RCP<Thyra::MultiVectorLinearOp<Scalar> > mv_adjoint_op =
361 RCP<Thyra::LinearOpBase<Scalar> > adjoint_op =
363 MEB::OutArgs<Scalar> adj_me_outArgs = adjoint_solve_model_->createOutArgs();
364 adj_me_outArgs.set_W_op(adjoint_op);
365 adjoint_solve_model_->evalModel(me_inArgs, adj_me_outArgs);
368 RCP<Thyra::VectorBase<Scalar> > adjoint_f = outArgs.
get_f();
369 RCP<Thyra::VectorBase<Scalar> > adjoint_g = outArgs.
get_g(0);
370 RCP<Thyra::VectorBase<Scalar> > g = outArgs.
get_g(1);
371 RCP<const Thyra::MultiVectorBase<Scalar> > adjoint_x_mv;
372 RCP<const Thyra::VectorBase<Scalar> > adjoint_x;
373 if (adjoint_f != Teuchos::null || adjoint_g != Teuchos::null) {
375 inArgs.
get_x().assert_not_null();
377 rcp_dynamic_cast<
const DMVPV>(adjoint_x,
true)->getMultiVector();
385 if (adjoint_f != Teuchos::null) {
386 TEMPUS_FUNC_TIME_MONITOR_DIFF(
"Tempus::AdjointSensitivityModelEvaluator::evalModel::f", TEMPUS_EVAL_F);
388 RCP<Thyra::MultiVectorBase<Scalar> > adjoint_f_mv =
389 rcp_dynamic_cast<
DMVPV>(adjoint_f,
true)->getNonconstMultiVector();
391 MEB::OutArgs<Scalar> me_outArgs = model_->createOutArgs();
392 MEB::OutArgs<Scalar> adj_me_outArgs = adjoint_residual_model_->createOutArgs();
397 TEMPUS_FUNC_TIME_MONITOR_DIFF(
"Tempus::AdjointSensitivityModelEvaluator::evalModel::dg/dx", TEMPUS_EVAL_DGDX);
398 if (my_dgdx_mv_ == Teuchos::null)
400 Thyra::createMembers(model_->get_x_space(),
401 model_->get_g_space(g_index_)->dim());
402 if (!response_gradient_is_computed_) {
403 me_outArgs.set_DgDx(g_index_,
404 MEB::Derivative<Scalar>(my_dgdx_mv_,
405 MEB::DERIV_MV_GRADIENT_FORM));
406 model_->evalModel(me_inArgs, me_outArgs);
407 me_outArgs.set_DgDx(g_index_, MEB::Derivative<Scalar>());
408 if (is_pseudotransient_)
409 response_gradient_is_computed_ =
true;
411 Thyra::assign(adjoint_f_mv.ptr(), *my_dgdx_mv_);
417 TEMPUS_FUNC_TIME_MONITOR_DIFF(
"Tempus::AdjointSensitivityModelEvaluator::evalModel::df/dx", TEMPUS_EVAL_DFDX);
418 if (my_dfdx_ == Teuchos::null)
419 my_dfdx_ = adjoint_residual_model_->create_W_op();
420 if (!jacobian_matrix_is_computed_) {
421 adj_me_outArgs.set_W_op(my_dfdx_);
422 if (me_inArgs.supports(MEB::IN_ARG_alpha))
423 me_inArgs.set_alpha(0.0);
424 if (me_inArgs.supports(MEB::IN_ARG_beta))
425 me_inArgs.set_beta(1.0);
426 adjoint_residual_model_->evalModel(me_inArgs, adj_me_outArgs);
427 if (is_pseudotransient_)
428 jacobian_matrix_is_computed_ =
true;
430 my_dfdx_->apply(
Thyra::NOTRANS, *adjoint_x_mv, adjoint_f_mv.ptr(),
431 Scalar(-1.0), Scalar(1.0));
438 if (me_inArgs.supports(MEB::IN_ARG_x_dot)) {
439 TEMPUS_FUNC_TIME_MONITOR_DIFF(
"Tempus::AdjointSensitivityModelEvaluator::evalModel::df/dx_dot", TEMPUS_EVAL_DFDXDOT);
440 RCP<const Thyra::VectorBase<Scalar> > adjoint_x_dot = inArgs.
get_x_dot();
441 if (adjoint_x_dot != Teuchos::null) {
442 RCP<const Thyra::MultiVectorBase<Scalar> > adjoint_x_dot_mv =
443 rcp_dynamic_cast<
const DMVPV>(adjoint_x_dot,
true)->getMultiVector();
444 if (mass_matrix_is_identity_) {
446 Thyra::V_StVpV(adjoint_f_mv.ptr(), Scalar(-1.0), *adjoint_f_mv,
450 if (my_dfdxdot_ == Teuchos::null)
451 my_dfdxdot_ = adjoint_residual_model_->create_W_op();
452 if (!mass_matrix_is_computed_) {
453 adj_me_outArgs.set_W_op(my_dfdxdot_);
454 me_inArgs.set_alpha(1.0);
455 me_inArgs.set_beta(0.0);
456 adjoint_residual_model_->evalModel(me_inArgs, adj_me_outArgs);
457 if (is_pseudotransient_ || mass_matrix_is_constant_)
458 mass_matrix_is_computed_ =
true;
461 adjoint_f_mv.ptr(), Scalar(1.0), Scalar(-1.0));
471 if (adjoint_g != Teuchos::null) {
472 TEMPUS_FUNC_TIME_MONITOR_DIFF(
"Tempus::AdjointSensitivityModelEvaluator::evalModel::g", TEMPUS_EVAL_G);
473 RCP<Thyra::MultiVectorBase<Scalar> > adjoint_g_mv =
474 rcp_dynamic_cast<
DMVPV>(adjoint_g,
true)->getNonconstMultiVector();
476 MEB::OutArgs<Scalar> me_outArgs = model_->createOutArgs();
479 MEB::DerivativeSupport dgdp_support =
480 me_outArgs.supports(MEB::OUT_ARG_DgDp, g_index_, p_index_);
481 if (dgdp_support.supports(MEB::DERIV_MV_GRADIENT_FORM)) {
482 me_outArgs.set_DgDp(g_index_, p_index_,
483 MEB::Derivative<Scalar>(adjoint_g_mv,
484 MEB::DERIV_MV_GRADIENT_FORM));
485 model_->evalModel(me_inArgs, me_outArgs);
487 else if (dgdp_support.supports(MEB::DERIV_MV_JACOBIAN_FORM)) {
488 const int num_g = model_->get_g_space(g_index_)->dim();
489 const int num_p = model_->get_p_space(p_index_)->dim();
490 RCP<Thyra::MultiVectorBase<Scalar> > dgdp_trans =
491 createMembers(model_->get_g_space(g_index_), num_p);
492 me_outArgs.set_DgDp(g_index_, p_index_,
493 MEB::Derivative<Scalar>(dgdp_trans,
494 MEB::DERIV_MV_JACOBIAN_FORM));
495 model_->evalModel(me_inArgs, me_outArgs);
498 for (
int i=0; i<num_p; ++i)
499 for (
int j=0; j<num_g; ++j)
500 dgdp_view(i,j) = dgdp_trans_view(j,i);
504 "Invalid dg/dp support");
505 me_outArgs.set_DgDp(g_index_, p_index_, MEB::Derivative<Scalar>());
508 MEB::DerivativeSupport dfdp_support =
509 me_outArgs.supports(MEB::OUT_ARG_DfDp, p_index_);
511 if (dfdp_support.supports(MEB::DERIV_LINEAR_OP)) {
512 if (my_dfdp_op_ == Teuchos::null)
513 my_dfdp_op_ = model_->create_DfDp_op(p_index_);
514 me_outArgs.set_DfDp(p_index_, MEB::Derivative<Scalar>(my_dfdp_op_));
517 else if (dfdp_support.supports(MEB::DERIV_MV_JACOBIAN_FORM)) {
518 if (my_dfdp_mv_ == Teuchos::null)
519 my_dfdp_mv_ = createMembers(model_->get_f_space(),
520 model_->get_p_space(p_index_)->dim());
521 me_outArgs.set_DfDp(p_index_,
522 MEB::Derivative<Scalar>(my_dfdp_mv_,
523 MEB::DERIV_MV_JACOBIAN_FORM));
524 my_dfdp_op_ = my_dfdp_mv_;
527 else if (dfdp_support.supports(MEB::DERIV_MV_GRADIENT_FORM)) {
528 if (my_dfdp_mv_ == Teuchos::null)
529 my_dfdp_mv_ = createMembers(model_->get_p_space(p_index_),
530 model_->get_f_space()->dim());
531 me_outArgs.set_DfDp(p_index_,
532 MEB::Derivative<Scalar>(my_dfdp_mv_,
533 MEB::DERIV_MV_GRADIENT_FORM));
534 my_dfdp_op_ = my_dfdp_mv_;
539 true, std::logic_error,
"Invalid df/dp support");
540 model_->evalModel(me_inArgs, me_outArgs);
541 my_dfdp_op_->apply(trans, *adjoint_x_mv, adjoint_g_mv.ptr(),
542 Scalar(-1.0), Scalar(1.0));
545 if (g != Teuchos::null) {
546 MEB::OutArgs<Scalar> me_outArgs = model_->createOutArgs();
547 me_outArgs.set_g(g_index_, g);
548 model_->evalModel(me_inArgs, me_outArgs);
552 template<
class Scalar>
558 pl->
set<
int>(
"Sensitivity Parameter Index", 0);
559 pl->
set<
int>(
"Response Function Index", 0);
560 pl->
set<
bool>(
"Mass Matrix Is Constant",
true);
561 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