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: Copyright (2017) Sandia Corporation
4 //
5 // Distributed under BSD 3-clause license (See accompanying file Copyright.txt)
6 // ****************************************************************************
7 // @HEADER
8 
9 #ifndef Tempus_AdjointAuxSensitivityModelEvaluator_impl_hpp
10 #define Tempus_AdjointAuxSensitivityModelEvaluator_impl_hpp
11 
16 #include "Thyra_DefaultBlockedLinearOp.hpp"
18 #include "Thyra_VectorStdOps.hpp"
19 #include "Thyra_MultiVectorStdOps.hpp"
20 
21 namespace Tempus {
22 
23 template <typename Scalar>
26  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > & model,
27  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > & adjoint_model,
28  const Scalar& t_init,
29  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::RCP;
39  using Teuchos::Array;
41  typedef Thyra::ModelEvaluatorBase MEB;
42 
43  // Set parameters
46  if (pList != Teuchos::null)
47  *pl = *pList;
49  mass_matrix_is_constant_ = pl->get<bool>("Mass Matrix Is Constant");
50  mass_matrix_is_identity_ = pl->get<bool>("Mass Matrix Is Identity");
51  p_index_ = pl->get<int>("Sensitivity Parameter Index", 0);
52  g_index_ = pl->get<int>("Response Function Index", 0);
53  num_adjoint_ = model_->get_g_space(g_index_)->dim();
54 
55  // We currently do not support a non-constant mass matrix
57  mass_matrix_is_constant_ == false, std::logic_error,
58  "AdjointAuxSensitivityModelEvaluator currently does not support " <<
59  "non-constant mass matrix df/dx_dot!");
60 
62  Thyra::multiVectorProductVectorSpace(model_->get_f_space(), num_adjoint_);
64  Thyra::multiVectorProductVectorSpace(model_->get_x_space(), num_adjoint_);
66  Thyra::multiVectorProductVectorSpace(model_->get_p_space(p_index_),
67  num_adjoint_);
68  Array< RCP<const VectorSpaceBase<Scalar> > > x_spaces(2);
69  Array< RCP<const VectorSpaceBase<Scalar> > > f_spaces(2);
70  x_spaces[0] = adjoint_space_;
71  x_spaces[1] = response_space_;
72  f_spaces[0] = residual_space_;
73  f_spaces[1] = response_space_;
74  x_prod_space_ = Thyra::productVectorSpace(x_spaces());
75  f_prod_space_ = Thyra::productVectorSpace(f_spaces());
76 
77  // forward and adjoint models must support same InArgs
78  MEB::InArgs<Scalar> me_inArgs = model_->createInArgs();
79  me_inArgs.assertSameSupport(adjoint_model_->createInArgs());
80 
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);
89 
90  // Support additional parameters for x and xdot
91  inArgs.set_Np(me_inArgs.Np());
92  prototypeInArgs_ = inArgs;
93 
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);
101  prototypeOutArgs_ = outArgs;
102 
103  // Adjoint ME must support W_op to define adjoint ODE/DAE.
104  // Must support alpha, beta if it suports x_dot
105  TEUCHOS_ASSERT(me_inArgs.supports(MEB::IN_ARG_x));
106  TEUCHOS_ASSERT(adj_me_outArgs.supports(MEB::OUT_ARG_W_op));
107  if (me_inArgs.supports(MEB::IN_ARG_x_dot)) {
108  TEUCHOS_ASSERT(me_inArgs.supports(MEB::IN_ARG_alpha));
109  TEUCHOS_ASSERT(me_inArgs.supports(MEB::IN_ARG_beta));
110  }
111 }
112 
113 template <typename Scalar>
114 void
116 setFinalTime(const Scalar t_final)
117 {
118  t_final_ = t_final;
119 }
120 
121 template <typename Scalar>
122 void
126 {
127  sh_ = sh;
129  forward_state_ = Teuchos::null;
130 }
131 
132 template <typename Scalar>
135 get_p_space(int p) const
136 {
137  TEUCHOS_ASSERT(p < model_->Np());
138  return model_->get_p_space(p);
139 }
140 
141 template <typename Scalar>
144 get_p_names(int p) const
145 {
146  TEUCHOS_ASSERT(p < model_->Np());
147  return model_->get_p_names(p);
148 }
149 
150 template <typename Scalar>
153 get_x_space() const
154 {
155  return x_prod_space_;
156 }
157 
158 template <typename Scalar>
161 get_f_space() const
162 {
163  return f_prod_space_;
164 }
165 
166 template <typename Scalar>
169 create_W_op() const
170 {
171  using Teuchos::RCP;
172  using Thyra::LinearOpBase;
173 
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);
182 }
183 
184 template <typename Scalar>
188 {
189  using Teuchos::RCP;
190  using Teuchos::rcp_dynamic_cast;
192 
193  RCP<const LOWSFB > alowsfb = adjoint_model_->get_W_factory();
194  if (alowsfb == Teuchos::null)
195  return Teuchos::null; // model_ doesn't support W_factory
196 
197  RCP<const LOWSFB > mv_alowsfb =
198  Thyra::multiVectorLinearOpWithSolveFactory(alowsfb, residual_space_,
199  adjoint_space_);
200 
201  RCP<const Thyra::VectorSpaceBase<Scalar> > g_space = response_space_;
202  RCP<const LOWSFB > g_lowsfb =
203  Thyra::scaledIdentitySolveFactory(g_space, Scalar(1.0));
204 
206  lowsfbs[0] = mv_alowsfb;
207  lowsfbs[1] = g_lowsfb;
208  return Thyra::blockedTriangularLinearOpWithSolveFactory(lowsfbs);
209 }
210 
211 template <typename Scalar>
215 {
216  return prototypeInArgs_;
217 }
218 
219 template <typename Scalar>
223 {
224  typedef Thyra::ModelEvaluatorBase MEB;
225  using Teuchos::RCP;
226  using Teuchos::rcp_dynamic_cast;
227 
228  MEB::InArgs<Scalar> me_nominal = model_->getNominalValues();
229  MEB::InArgs<Scalar> nominal = this->createInArgs();
230 
231  const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
232 
233  // Set initial x, x_dot
234  RCP< Thyra::VectorBase<Scalar> > x = Thyra::createMember(*x_prod_space_);
235  Thyra::assign(x.ptr(), zero);
236  nominal.set_x(x);
237 
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);
243  }
244 
245  const int np = model_->Np();
246  for (int i=0; i<np; ++i)
247  nominal.set_p(i, me_nominal.get_p(i));
248 
249  return nominal;
250 }
251 
252 template <typename Scalar>
256 {
257  return prototypeOutArgs_;
258 }
259 
260 template <typename Scalar>
261 void
264  const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
265 {
266  typedef Thyra::ModelEvaluatorBase MEB;
267  using Teuchos::RCP;
268  using Teuchos::rcp_dynamic_cast;
269 
270  // Note: adjoint_model computes the transposed W (either explicitly or
271  // implicitly. Thus we need to always call adjoint_model->evalModel()
272  // whenever computing the adjoint operator, and subsequent calls to apply()
273  // do not transpose it.
274 
275  // Interpolate forward solution at supplied time, reusing previous
276  // interpolation if possible
277  TEUCHOS_ASSERT(sh_ != Teuchos::null);
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);
283  else
284  sh_->interpolateState(forward_t, forward_state_.get());
285  t_interp_ = t;
286  }
287 
288  // setup input arguments for model
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) &&
292  inArgs.get_x_dot() != Teuchos::null)
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));
299 
300  // compute W
301  RCP<Thyra::LinearOpBase<Scalar> > op = outArgs.get_W_op();
302  if (op != Teuchos::null) {
303  if (me_inArgs.supports(MEB::IN_ARG_alpha))
304  me_inArgs.set_alpha(inArgs.get_alpha());
305  if (me_inArgs.supports(MEB::IN_ARG_beta))
306  me_inArgs.set_beta(inArgs.get_beta());
307 
308  // Adjoint W
309  RCP<Thyra::DefaultBlockedLinearOp<Scalar> > block_op =
310  rcp_dynamic_cast<Thyra::DefaultBlockedLinearOp<Scalar> >(op,true);
311  RCP<Thyra::MultiVectorLinearOp<Scalar> > mv_adjoint_op =
312  rcp_dynamic_cast<Thyra::MultiVectorLinearOp<Scalar> >(
313  block_op->getNonconstBlock(0,0),true);
314  RCP<Thyra::LinearOpBase<Scalar> > adjoint_op =
315  mv_adjoint_op->getNonconstLinearOp();
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);
319 
320  // g W
321  RCP<Thyra::ScaledIdentityLinearOpWithSolve<Scalar> > si_op =
323  block_op->getNonconstBlock(1,1),true);
324  si_op->setScale(inArgs.get_alpha());
325  }
326 
327  // Compute adjoint residual F(y):
328  // * For implicit form, F(y) = d/dt( df/dx_dot^T*y ) + df/dx^T*y
329  // * For explict form, F(y) = -df/dx^T*y
330  // For implicit form, we assume df/dx_dot is constant w.r.t. x, x_dot, and t,
331  // so the residual becomes F(y) = df/dx_dot^T*y_dot + df/dx^T*y
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();
339 
340  RCP<DPV> prod_f = rcp_dynamic_cast<DPV>(f,true);
341  RCP<Thyra::VectorBase<Scalar> > adjoint_f =
342  prod_f->getNonconstVectorBlock(0);
343  RCP<Thyra::MultiVectorBase<Scalar> > adjoint_f_mv =
344  rcp_dynamic_cast<DMVPV>(adjoint_f,true)->getNonconstMultiVector();
345 
346  MEB::OutArgs<Scalar> adj_me_outArgs = adjoint_model_->createOutArgs();
347 
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);
356 
357  // Explicit form residual F(y) = -df/dx^T*y
358  my_dfdx_->apply(Thyra::NOTRANS, *adjoint_x_mv, adjoint_f_mv.ptr(),
359  Scalar(-1.0), Scalar(0.0));
360 
361  // Implicit form residual df/dx_dot^T*y_dot + df/dx^T*y using the second
362  // scalar argument to apply() to change the explicit term above
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 =
369  prod_x_dot->getVectorBlock(0);
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_) {
373  // F = -F + y_dot
374  Thyra::V_StVpV(adjoint_f_mv.ptr(), Scalar(-1.0), *adjoint_f_mv,
375  *adjoint_x_dot_mv);
376  }
377  else {
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);
385 
386  mass_matrix_is_computed_ = true;
387  }
388  my_dfdxdot_->apply(Thyra::NOTRANS, *adjoint_x_dot_mv,
389  adjoint_f_mv.ptr(), Scalar(1.0), Scalar(-1.0));
390  }
391  }
392  }
393 
394  // Compute g = z_dot - df/dp^T*y for computing the model parameter term
395  // in the adjoint sensitivity formula
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();
400 
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_));
409  trans = Thyra::CONJTRANS;
410  }
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_;
419  trans = Thyra::CONJTRANS;
420  }
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_;
429  trans = Thyra::CONJ;
430  }
431  else
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));
437 
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);
444  }
445  }
446 }
447 
448 template<class Scalar>
452 {
453  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
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);
458  return pl;
459 }
460 
461 } // namespace Tempus
462 
463 #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