9#ifndef Tempus_AdjointSensitivityModelEvaluator_impl_hpp
10#define Tempus_AdjointSensitivityModelEvaluator_impl_hpp
12#include "Teuchos_TimeMonitor.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"
27template <
typename Scalar>
34 const Scalar& t_final,
35 const bool is_pseudotransient,
36 const Teuchos::RCP<const Teuchos::ParameterList>& pList) :
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())
48 typedef Thyra::ModelEvaluatorBase MEB;
51 Teuchos::RCP<Teuchos::ParameterList> pl =
52 Teuchos::rcp(
new Teuchos::ParameterList);
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);
63 TEUCHOS_TEST_FOR_EXCEPTION(
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);
109 TEUCHOS_ASSERT(me_inArgs.supports(MEB::IN_ARG_x));
110 TEUCHOS_ASSERT(adj_mer_outArgs.supports(MEB::OUT_ARG_W_op));
111 if (me_inArgs.supports(MEB::IN_ARG_x_dot)) {
112 TEUCHOS_ASSERT(me_inArgs.supports(MEB::IN_ARG_alpha));
113 TEUCHOS_ASSERT(me_inArgs.supports(MEB::IN_ARG_beta));
115 MEB::DerivativeSupport dgdx_support =
116 me_outArgs.supports(MEB::OUT_ARG_DgDx,
g_index_);
117 MEB::DerivativeSupport dgdp_support =
119 TEUCHOS_ASSERT(dgdx_support.supports(MEB::DERIV_MV_GRADIENT_FORM));
120 TEUCHOS_ASSERT(dgdp_support.supports(MEB::DERIV_MV_GRADIENT_FORM));
123template <
typename Scalar>
131template <
typename Scalar>
138 if (is_pseudotransient_)
139 forward_state_ = sh_->getCurrentState();
141 t_interp_ = Teuchos::ScalarTraits<Scalar>::rmax();
142 forward_state_ = Teuchos::null;
146 mass_matrix_is_computed_ =
false;
147 jacobian_matrix_is_computed_ =
false;
148 response_gradient_is_computed_ =
false;
151template <
typename Scalar>
152Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
156 TEUCHOS_ASSERT(p < model_->Np());
157 return model_->get_p_space(p);
160template <
typename Scalar>
161Teuchos::RCP<const Teuchos::Array<std::string> >
165 TEUCHOS_ASSERT(p < model_->Np());
166 return model_->get_p_names(p);
169template <
typename Scalar>
170Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
174 return adjoint_space_;
177template <
typename Scalar>
178Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
182 return residual_space_;
185template <
typename Scalar>
186Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
190 TEUCHOS_ASSERT(j == 0 || j == 1);
192 return response_space_;
193 return model_->get_g_space(g_index_);
196template <
typename Scalar>
197Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
201 Teuchos::RCP<Thyra::LinearOpBase<Scalar> > adjoint_op =
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_);
209template <
typename Scalar>
210Teuchos::RCP<const Thyra::LinearOpWithSolveFactoryBase<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_);
226template <
typename Scalar>
227Thyra::ModelEvaluatorBase::InArgs<Scalar>
231 return prototypeInArgs_;
234template <
typename Scalar>
235Thyra::ModelEvaluatorBase::InArgs<Scalar>
239 typedef Thyra::ModelEvaluatorBase MEB;
241 using Teuchos::rcp_dynamic_cast;
243 MEB::InArgs<Scalar> me_nominal = model_->getNominalValues();
244 MEB::InArgs<Scalar> nominal = this->createInArgs();
246 const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
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));
267template <
typename Scalar>
268Thyra::ModelEvaluatorBase::OutArgs<Scalar>
272 return prototypeOutArgs_;
275template <
typename Scalar>
278evalModelImpl(
const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
279 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs)
const
281 typedef Thyra::ModelEvaluatorBase MEB;
283 using Teuchos::rcp_dynamic_cast;
285 TEMPUS_FUNC_TIME_MONITOR_DIFF(
"Tempus::AdjointSensitivityModelEvaluator::evalModel()", TEMPUS_EVAL);
294 TEUCHOS_ASSERT(sh_ != Teuchos::null);
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)) {
314 if (inArgs.get_x_dot() != Teuchos::null)
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))
347 op = outArgs.get_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))
351 me_inArgs.set_alpha(inArgs.get_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 =
360 rcp_dynamic_cast<Thyra::MultiVectorLinearOp<Scalar> >(op,
true);
361 RCP<Thyra::LinearOpBase<Scalar> > adjoint_op =
362 mv_adjoint_op->getNonconstLinearOp();
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;
460 my_dfdxdot_->apply(Thyra::NOTRANS, *adjoint_x_dot_mv,
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);
496 Thyra::DetachedMultiVectorView<Scalar> dgdp_view(*adjoint_g_mv);
497 Thyra::DetachedMultiVectorView<Scalar> dgdp_trans_view(*dgdp_trans);
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);
503 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
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_);
510 Thyra::EOpTransp trans = Thyra::CONJTRANS;
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_));
515 trans = Thyra::CONJTRANS;
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_;
525 trans = Thyra::CONJTRANS;
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_;
538 TEUCHOS_TEST_FOR_EXCEPTION(
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);
552template<
class Scalar>
553Teuchos::RCP<const Teuchos::ParameterList>
557 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
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);
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_p_space(int p) const
Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > adjoint_residual_model_
Thyra::ModelEvaluatorBase::InArgs< Scalar > getNominalValues() const
bool mass_matrix_is_identity_
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
void setForwardSolutionHistory(const Teuchos::RCP< const Tempus::SolutionHistory< Scalar > > &sh)
Set solution history from forward evaluation.
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_f_space() const
Thyra::ModelEvaluatorBase::OutArgs< Scalar > createOutArgsImpl() const
static Teuchos::RCP< const Teuchos::ParameterList > getValidParameters()
Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > model_
Teuchos::RCP< const DMVPVS > residual_space_
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.
Thyra::ModelEvaluatorBase::InArgs< Scalar > prototypeInArgs_
Teuchos::RCP< const Thyra::LinearOpWithSolveFactoryBase< Scalar > > get_W_factory() const
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int p) const
Teuchos::RCP< const DMVPVS > response_space_
Thyra::ModelEvaluatorBase::OutArgs< Scalar > prototypeOutArgs_
Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > adjoint_solve_model_
void setFinalTime(const Scalar t_final)
Set the final time from the forward evaluation.
bool mass_matrix_is_constant_
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_g_space(int j) const
Teuchos::RCP< const DMVPVS > adjoint_space_
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_W_op() const
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...