Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tempus_AdjointAuxSensitivityModelEvaluator_impl.hpp
Go to the documentation of this file.
1 //@HEADER
2 // *****************************************************************************
3 // Tempus: Time Integration and Sensitivity Analysis Package
4 //
5 // Copyright 2017 NTESS and the Tempus contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 //@HEADER
9 
10 #ifndef Tempus_AdjointAuxSensitivityModelEvaluator_impl_hpp
11 #define Tempus_AdjointAuxSensitivityModelEvaluator_impl_hpp
12 
17 #include "Thyra_DefaultBlockedLinearOp.hpp"
19 #include "Thyra_VectorStdOps.hpp"
20 #include "Thyra_MultiVectorStdOps.hpp"
21 
22 namespace Tempus {
23 
24 template <typename Scalar>
27  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
28  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& adjoint_model,
29  const Scalar& t_init, const Scalar& t_final,
31  : model_(model),
32  adjoint_model_(adjoint_model),
33  t_init_(t_init),
34  t_final_(t_final),
35  mass_matrix_is_computed_(false),
36  t_interp_(Teuchos::ScalarTraits<Scalar>::rmax())
37 {
38  using Teuchos::Array;
39  using Teuchos::RCP;
41  typedef Thyra::ModelEvaluatorBase MEB;
42 
43  // Set parameters
46  if (pList != Teuchos::null) *pl = *pList;
48  mass_matrix_is_constant_ = pl->get<bool>("Mass Matrix Is Constant");
49  mass_matrix_is_identity_ = pl->get<bool>("Mass Matrix Is Identity");
50  p_index_ = pl->get<int>("Sensitivity Parameter Index", 0);
51  g_index_ = pl->get<int>("Response Function Index", 0);
52  num_adjoint_ = model_->get_g_space(g_index_)->dim();
53 
54  // We currently do not support a non-constant mass matrix
56  mass_matrix_is_constant_ == false, std::logic_error,
57  "AdjointAuxSensitivityModelEvaluator currently does not support "
58  << "non-constant mass matrix df/dx_dot!");
59 
61  Thyra::multiVectorProductVectorSpace(model_->get_f_space(), num_adjoint_);
63  Thyra::multiVectorProductVectorSpace(model_->get_x_space(), num_adjoint_);
64  response_space_ = Thyra::multiVectorProductVectorSpace(
65  model_->get_p_space(p_index_), num_adjoint_);
66  Array<RCP<const VectorSpaceBase<Scalar> > > x_spaces(2);
67  Array<RCP<const VectorSpaceBase<Scalar> > > f_spaces(2);
68  x_spaces[0] = adjoint_space_;
69  x_spaces[1] = response_space_;
70  f_spaces[0] = residual_space_;
71  f_spaces[1] = response_space_;
72  x_prod_space_ = Thyra::productVectorSpace(x_spaces());
73  f_prod_space_ = Thyra::productVectorSpace(f_spaces());
74 
75  // forward and adjoint models must support same InArgs
76  MEB::InArgs<Scalar> me_inArgs = model_->createInArgs();
77  me_inArgs.assertSameSupport(adjoint_model_->createInArgs());
78 
79  MEB::InArgsSetup<Scalar> inArgs;
80  inArgs.setModelEvalDescription(this->description());
81  inArgs.setSupports(MEB::IN_ARG_x);
82  inArgs.setSupports(MEB::IN_ARG_t);
83  if (me_inArgs.supports(MEB::IN_ARG_x_dot))
84  inArgs.setSupports(MEB::IN_ARG_x_dot);
85  inArgs.setSupports(MEB::IN_ARG_alpha);
86  inArgs.setSupports(MEB::IN_ARG_beta);
87 
88  // Support additional parameters for x and xdot
89  inArgs.set_Np(me_inArgs.Np());
90  prototypeInArgs_ = inArgs;
91 
92  MEB::OutArgs<Scalar> me_outArgs = model_->createOutArgs();
93  MEB::OutArgs<Scalar> adj_me_outArgs = adjoint_model_->createOutArgs();
94  MEB::OutArgsSetup<Scalar> outArgs;
95  outArgs.setModelEvalDescription(this->description());
96  outArgs.set_Np_Ng(me_inArgs.Np(), 0);
97  outArgs.setSupports(MEB::OUT_ARG_f);
98  outArgs.setSupports(MEB::OUT_ARG_W_op);
99  prototypeOutArgs_ = outArgs;
100 
101  // Adjoint ME must support W_op to define adjoint ODE/DAE.
102  // Must support alpha, beta if it suports x_dot
103  TEUCHOS_ASSERT(me_inArgs.supports(MEB::IN_ARG_x));
104  TEUCHOS_ASSERT(adj_me_outArgs.supports(MEB::OUT_ARG_W_op));
105  if (me_inArgs.supports(MEB::IN_ARG_x_dot)) {
106  TEUCHOS_ASSERT(me_inArgs.supports(MEB::IN_ARG_alpha));
107  TEUCHOS_ASSERT(me_inArgs.supports(MEB::IN_ARG_beta));
108  }
109 }
110 
111 template <typename Scalar>
113  const Scalar t_final)
114 {
115  t_final_ = t_final;
116 }
117 
118 template <typename Scalar>
121 {
122  sh_ = sh;
124  forward_state_ = Teuchos::null;
125 }
126 
127 template <typename Scalar>
130 {
131  TEUCHOS_ASSERT(p < model_->Np());
132  return model_->get_p_space(p);
133 }
134 
135 template <typename Scalar>
138 {
139  TEUCHOS_ASSERT(p < model_->Np());
140  return model_->get_p_names(p);
141 }
142 
143 template <typename Scalar>
146 {
147  return x_prod_space_;
148 }
149 
150 template <typename Scalar>
153 {
154  return f_prod_space_;
155 }
156 
157 template <typename Scalar>
160 {
161  using Teuchos::RCP;
162  using Thyra::LinearOpBase;
163 
164  RCP<LinearOpBase<Scalar> > adjoint_op = adjoint_model_->create_W_op();
165  RCP<LinearOpBase<Scalar> > mv_adjoint_op =
166  Thyra::nonconstMultiVectorLinearOp(adjoint_op, num_adjoint_);
167  RCP<const Thyra::VectorSpaceBase<Scalar> > g_space = response_space_;
168  RCP<LinearOpBase<Scalar> > g_op = Thyra::scaledIdentity(g_space, Scalar(1.0));
169  RCP<LinearOpBase<Scalar> > null_op;
170  return nonconstBlock2x2(mv_adjoint_op, null_op, null_op, g_op);
171 }
172 
173 template <typename Scalar>
176 {
177  using Teuchos::RCP;
178  using Teuchos::rcp_dynamic_cast;
180 
181  RCP<const LOWSFB> alowsfb = adjoint_model_->get_W_factory();
182  if (alowsfb == Teuchos::null)
183  return Teuchos::null; // model_ doesn't support W_factory
184 
185  RCP<const LOWSFB> mv_alowsfb = Thyra::multiVectorLinearOpWithSolveFactory(
186  alowsfb, residual_space_, adjoint_space_);
187 
188  RCP<const Thyra::VectorSpaceBase<Scalar> > g_space = response_space_;
189  RCP<const LOWSFB> g_lowsfb =
190  Thyra::scaledIdentitySolveFactory(g_space, Scalar(1.0));
191 
193  lowsfbs[0] = mv_alowsfb;
194  lowsfbs[1] = g_lowsfb;
195  return Thyra::blockedTriangularLinearOpWithSolveFactory(lowsfbs);
196 }
197 
198 template <typename Scalar>
201 {
202  return prototypeInArgs_;
203 }
204 
205 template <typename Scalar>
208 {
209  typedef Thyra::ModelEvaluatorBase MEB;
210  using Teuchos::RCP;
211  using Teuchos::rcp_dynamic_cast;
212 
213  MEB::InArgs<Scalar> me_nominal = model_->getNominalValues();
214  MEB::InArgs<Scalar> nominal = this->createInArgs();
215 
216  const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
217 
218  // Set initial x, x_dot
219  RCP<Thyra::VectorBase<Scalar> > x = Thyra::createMember(*x_prod_space_);
220  Thyra::assign(x.ptr(), zero);
221  nominal.set_x(x);
222 
223  if (me_nominal.supports(MEB::IN_ARG_x_dot)) {
224  RCP<Thyra::VectorBase<Scalar> > x_dot = Thyra::createMember(*x_prod_space_);
225  Thyra::assign(x_dot.ptr(), zero);
226  nominal.set_x_dot(x_dot);
227  }
228 
229  const int np = model_->Np();
230  for (int i = 0; i < np; ++i) nominal.set_p(i, me_nominal.get_p(i));
231 
232  return nominal;
233 }
234 
235 template <typename Scalar>
238 {
239  return prototypeOutArgs_;
240 }
241 
242 template <typename Scalar>
245  const Thyra::ModelEvaluatorBase::OutArgs<Scalar>& outArgs) const
246 {
247  typedef Thyra::ModelEvaluatorBase MEB;
248  using Teuchos::RCP;
249  using Teuchos::rcp_dynamic_cast;
250 
251  // Note: adjoint_model computes the transposed W (either explicitly or
252  // implicitly. Thus we need to always call adjoint_model->evalModel()
253  // whenever computing the adjoint operator, and subsequent calls to apply()
254  // do not transpose it.
255 
256  // Interpolate forward solution at supplied time, reusing previous
257  // interpolation if possible
258  TEUCHOS_ASSERT(sh_ != Teuchos::null);
259  const Scalar t = inArgs.get_t();
260  const Scalar forward_t = t_final_ + t_init_ - t;
261  if (forward_state_ == Teuchos::null || t_interp_ != t) {
262  if (forward_state_ == Teuchos::null)
263  forward_state_ = sh_->interpolateState(forward_t);
264  else
265  sh_->interpolateState(forward_t, forward_state_.get());
266  t_interp_ = t;
267  }
268 
269  // setup input arguments for model
270  MEB::InArgs<Scalar> me_inArgs = model_->getNominalValues();
271  me_inArgs.set_x(forward_state_->getX());
272  if (me_inArgs.supports(MEB::IN_ARG_x_dot) &&
273  inArgs.get_x_dot() != Teuchos::null)
274  me_inArgs.set_x_dot(forward_state_->getXDot());
275  if (me_inArgs.supports(MEB::IN_ARG_t)) me_inArgs.set_t(forward_t);
276  const int np = me_inArgs.Np();
277  for (int i = 0; i < np; ++i) me_inArgs.set_p(i, inArgs.get_p(i));
278 
279  // compute W
280  RCP<Thyra::LinearOpBase<Scalar> > op = outArgs.get_W_op();
281  if (op != Teuchos::null) {
282  if (me_inArgs.supports(MEB::IN_ARG_alpha))
283  me_inArgs.set_alpha(inArgs.get_alpha());
284  if (me_inArgs.supports(MEB::IN_ARG_beta))
285  me_inArgs.set_beta(inArgs.get_beta());
286 
287  // Adjoint W
288  RCP<Thyra::DefaultBlockedLinearOp<Scalar> > block_op =
289  rcp_dynamic_cast<Thyra::DefaultBlockedLinearOp<Scalar> >(op, true);
290  RCP<Thyra::MultiVectorLinearOp<Scalar> > mv_adjoint_op =
291  rcp_dynamic_cast<Thyra::MultiVectorLinearOp<Scalar> >(
292  block_op->getNonconstBlock(0, 0), true);
293  RCP<Thyra::LinearOpBase<Scalar> > adjoint_op =
294  mv_adjoint_op->getNonconstLinearOp();
295  MEB::OutArgs<Scalar> adj_me_outArgs = adjoint_model_->createOutArgs();
296  adj_me_outArgs.set_W_op(adjoint_op);
297  adjoint_model_->evalModel(me_inArgs, adj_me_outArgs);
298 
299  // g W
300  RCP<Thyra::ScaledIdentityLinearOpWithSolve<Scalar> > si_op =
302  block_op->getNonconstBlock(1, 1), true);
303  si_op->setScale(inArgs.get_alpha());
304  }
305 
306  // Compute adjoint residual F(y):
307  // * For implicit form, F(y) = d/dt( df/dx_dot^T*y ) + df/dx^T*y
308  // * For explict form, F(y) = -df/dx^T*y
309  // For implicit form, we assume df/dx_dot is constant w.r.t. x, x_dot, and t,
310  // so the residual becomes F(y) = df/dx_dot^T*y_dot + df/dx^T*y
311  RCP<Thyra::VectorBase<Scalar> > f = outArgs.get_f();
312  if (f != Teuchos::null) {
313  RCP<const Thyra::VectorBase<Scalar> > x = inArgs.get_x().assert_not_null();
314  RCP<const DPV> prod_x = rcp_dynamic_cast<const DPV>(x, true);
315  RCP<const Thyra::VectorBase<Scalar> > adjoint_x = prod_x->getVectorBlock(0);
316  RCP<const Thyra::MultiVectorBase<Scalar> > adjoint_x_mv =
317  rcp_dynamic_cast<const DMVPV>(adjoint_x, true)->getMultiVector();
318 
319  RCP<DPV> prod_f = rcp_dynamic_cast<DPV>(f, true);
320  RCP<Thyra::VectorBase<Scalar> > adjoint_f =
321  prod_f->getNonconstVectorBlock(0);
322  RCP<Thyra::MultiVectorBase<Scalar> > adjoint_f_mv =
323  rcp_dynamic_cast<DMVPV>(adjoint_f, true)->getNonconstMultiVector();
324 
325  MEB::OutArgs<Scalar> adj_me_outArgs = adjoint_model_->createOutArgs();
326 
327  if (my_dfdx_ == Teuchos::null) my_dfdx_ = adjoint_model_->create_W_op();
328  adj_me_outArgs.set_W_op(my_dfdx_);
329  if (me_inArgs.supports(MEB::IN_ARG_alpha)) me_inArgs.set_alpha(0.0);
330  if (me_inArgs.supports(MEB::IN_ARG_beta)) me_inArgs.set_beta(1.0);
331  adjoint_model_->evalModel(me_inArgs, adj_me_outArgs);
332 
333  // Explicit form residual F(y) = -df/dx^T*y
334  my_dfdx_->apply(Thyra::NOTRANS, *adjoint_x_mv, adjoint_f_mv.ptr(),
335  Scalar(-1.0), Scalar(0.0));
336 
337  // Implicit form residual df/dx_dot^T*y_dot + df/dx^T*y using the second
338  // scalar argument to apply() to change the explicit term above
339  RCP<const DPV> prod_x_dot;
340  if (me_inArgs.supports(MEB::IN_ARG_x_dot)) {
341  RCP<const Thyra::VectorBase<Scalar> > x_dot = inArgs.get_x_dot();
342  if (x_dot != Teuchos::null) {
343  prod_x_dot = rcp_dynamic_cast<const DPV>(x_dot, true);
344  RCP<const Thyra::VectorBase<Scalar> > adjoint_x_dot =
345  prod_x_dot->getVectorBlock(0);
346  RCP<const Thyra::MultiVectorBase<Scalar> > adjoint_x_dot_mv =
347  rcp_dynamic_cast<const DMVPV>(adjoint_x_dot, true)
348  ->getMultiVector();
349  if (mass_matrix_is_identity_) {
350  // F = -F + y_dot
351  Thyra::V_StVpV(adjoint_f_mv.ptr(), Scalar(-1.0), *adjoint_f_mv,
352  *adjoint_x_dot_mv);
353  }
354  else {
355  if (my_dfdxdot_ == Teuchos::null)
356  my_dfdxdot_ = adjoint_model_->create_W_op();
357  if (!mass_matrix_is_constant_ || !mass_matrix_is_computed_) {
358  adj_me_outArgs.set_W_op(my_dfdxdot_);
359  me_inArgs.set_alpha(1.0);
360  me_inArgs.set_beta(0.0);
361  adjoint_model_->evalModel(me_inArgs, adj_me_outArgs);
362 
363  mass_matrix_is_computed_ = true;
364  }
365  my_dfdxdot_->apply(Thyra::NOTRANS, *adjoint_x_dot_mv,
366  adjoint_f_mv.ptr(), Scalar(1.0), Scalar(-1.0));
367  }
368  }
369  }
370 
371  // Compute g = z_dot - df/dp^T*y for computing the model parameter term
372  // in the adjoint sensitivity formula
373  RCP<Thyra::VectorBase<Scalar> > adjoint_g =
374  prod_f->getNonconstVectorBlock(1);
375  RCP<Thyra::MultiVectorBase<Scalar> > adjoint_g_mv =
376  rcp_dynamic_cast<DMVPV>(adjoint_g, true)->getNonconstMultiVector();
377 
378  MEB::OutArgs<Scalar> me_outArgs2 = model_->createOutArgs();
379  MEB::DerivativeSupport dfdp_support =
380  me_outArgs2.supports(MEB::OUT_ARG_DfDp, p_index_);
382  if (dfdp_support.supports(MEB::DERIV_LINEAR_OP)) {
383  if (my_dfdp_op_ == Teuchos::null)
384  my_dfdp_op_ = model_->create_DfDp_op(p_index_);
385  me_outArgs2.set_DfDp(p_index_, MEB::Derivative<Scalar>(my_dfdp_op_));
386  trans = Thyra::CONJTRANS;
387  }
388  else if (dfdp_support.supports(MEB::DERIV_MV_JACOBIAN_FORM)) {
389  if (my_dfdp_mv_ == Teuchos::null)
390  my_dfdp_mv_ = createMembers(model_->get_f_space(),
391  model_->get_p_space(p_index_)->dim());
392  me_outArgs2.set_DfDp(
393  p_index_,
394  MEB::Derivative<Scalar>(my_dfdp_mv_, MEB::DERIV_MV_JACOBIAN_FORM));
395  my_dfdp_op_ = my_dfdp_mv_;
396  trans = Thyra::CONJTRANS;
397  }
398  else if (dfdp_support.supports(MEB::DERIV_MV_GRADIENT_FORM)) {
399  if (my_dfdp_mv_ == Teuchos::null)
400  my_dfdp_mv_ = createMembers(model_->get_p_space(p_index_),
401  model_->get_f_space()->dim());
402  me_outArgs2.set_DfDp(
403  p_index_,
404  MEB::Derivative<Scalar>(my_dfdp_mv_, MEB::DERIV_MV_GRADIENT_FORM));
405  my_dfdp_op_ = my_dfdp_mv_;
406  trans = Thyra::CONJ;
407  }
408  else
409  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
410  "Invalid df/dp support");
411  model_->evalModel(me_inArgs, me_outArgs2);
412  my_dfdp_op_->apply(trans, *adjoint_x_mv, adjoint_g_mv.ptr(), Scalar(1.0),
413  Scalar(0.0));
414 
415  if (prod_x_dot != Teuchos::null) {
416  RCP<const Thyra::VectorBase<Scalar> > z_dot =
417  prod_x_dot->getVectorBlock(1);
418  RCP<const Thyra::MultiVectorBase<Scalar> > z_dot_mv =
419  rcp_dynamic_cast<const DMVPV>(z_dot, true)->getMultiVector();
420  Thyra::V_VmV(adjoint_g_mv.ptr(), *z_dot_mv, *adjoint_g_mv);
421  }
422  }
423 }
424 
425 template <class Scalar>
428 {
429  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
430  pl->set<int>("Sensitivity Parameter Index", 0);
431  pl->set<int>("Response Function Index", 0);
432  pl->set<bool>("Mass Matrix Is Constant", true);
433  pl->set<bool>("Mass Matrix Is Identity", false);
434  return pl;
435 }
436 
437 } // namespace Tempus
438 
439 #endif
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...
EOpTransp
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
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()
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
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
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.
#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