Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tempus_AdjointSensitivityModelEvaluator_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_AdjointSensitivityModelEvaluator_impl_hpp
10 #define Tempus_AdjointSensitivityModelEvaluator_impl_hpp
11 
12 #include "Teuchos_TimeMonitor.hpp"
13 #include "Tempus_config.hpp"
14 
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"
24 
25 namespace Tempus {
26 
27 template <typename Scalar>
29  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
31  adjoint_residual_model,
33  adjoint_solve_model,
34  const Scalar& t_init, const Scalar& t_final, const bool is_pseudotransient,
36  : model_(model),
37  adjoint_residual_model_(adjoint_residual_model),
38  adjoint_solve_model_(adjoint_solve_model),
39  t_init_(t_init),
40  t_final_(t_final),
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())
46 {
47  typedef Thyra::ModelEvaluatorBase MEB;
48 
49  // Set parameters
52  if (pList != Teuchos::null) *pl = *pList;
54  mass_matrix_is_constant_ = pl->get<bool>("Mass Matrix Is Constant");
55  mass_matrix_is_identity_ = pl->get<bool>("Mass Matrix Is Identity");
56  p_index_ = pl->get<int>("Sensitivity Parameter Index", 0);
57  g_index_ = pl->get<int>("Response Function Index", 0);
58  num_adjoint_ = model_->get_g_space(g_index_)->dim();
59 
60  // We currently do not support a non-constant mass matrix
62  mass_matrix_is_constant_ == false, std::logic_error,
63  "AdjointSensitivityModelEvaluator currently does not support "
64  << "non-constant mass matrix df/dx_dot!");
65 
67  Thyra::multiVectorProductVectorSpace(model_->get_f_space(), num_adjoint_);
69  Thyra::multiVectorProductVectorSpace(model_->get_x_space(), num_adjoint_);
70  response_space_ = Thyra::multiVectorProductVectorSpace(
71  model_->get_p_space(p_index_), num_adjoint_);
72 
73  // forward and adjoint models must support same InArgs
74  MEB::InArgs<Scalar> me_inArgs = model_->createInArgs();
75  me_inArgs.assertSameSupport(adjoint_residual_model_->createInArgs());
76  me_inArgs.assertSameSupport(adjoint_solve_model_->createInArgs());
77 
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);
88 
89  // Support additional parameters for x and xdot
90  inArgs.set_Np(me_inArgs.Np());
91  prototypeInArgs_ = inArgs;
92 
93  MEB::OutArgs<Scalar> me_outArgs = model_->createOutArgs();
94  MEB::OutArgs<Scalar> adj_mer_outArgs =
95  adjoint_residual_model_->createOutArgs();
96  MEB::OutArgs<Scalar> adj_mes_outArgs = adjoint_solve_model_->createOutArgs();
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);
103  prototypeOutArgs_ = outArgs;
104 
105  // Adjoint residual ME must support W_op to define adjoint ODE/DAE.
106  // Must support alpha, beta if it suports x_dot
107  TEUCHOS_ASSERT(me_inArgs.supports(MEB::IN_ARG_x));
108  TEUCHOS_ASSERT(adj_mer_outArgs.supports(MEB::OUT_ARG_W_op));
109  if (me_inArgs.supports(MEB::IN_ARG_x_dot)) {
110  TEUCHOS_ASSERT(me_inArgs.supports(MEB::IN_ARG_alpha));
111  TEUCHOS_ASSERT(me_inArgs.supports(MEB::IN_ARG_beta));
112  }
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));
119 }
120 
121 template <typename Scalar>
123  const Scalar t_final)
124 {
125  t_final_ = t_final;
126 }
127 
128 template <typename Scalar>
131 {
132  sh_ = sh;
133  if (is_pseudotransient_)
134  forward_state_ = sh_->getCurrentState();
135  else {
137  forward_state_ = Teuchos::null;
138  }
139 
140  // Reset computation flags because we have done a new forward integration
141  mass_matrix_is_computed_ = false;
142  jacobian_matrix_is_computed_ = false;
143  response_gradient_is_computed_ = false;
144 }
145 
146 template <typename Scalar>
149 {
150  TEUCHOS_ASSERT(p < model_->Np());
151  return model_->get_p_space(p);
152 }
153 
154 template <typename Scalar>
157 {
158  TEUCHOS_ASSERT(p < model_->Np());
159  return model_->get_p_names(p);
160 }
161 
162 template <typename Scalar>
165 {
166  return adjoint_space_;
167 }
168 
169 template <typename Scalar>
172 {
173  return residual_space_;
174 }
175 
176 template <typename Scalar>
179 {
180  TEUCHOS_ASSERT(j == 0 || j == 1);
181  if (j == 0) return response_space_;
182  return model_->get_g_space(g_index_);
183 }
184 
185 template <typename Scalar>
188 {
190  adjoint_solve_model_->create_W_op();
191  if (adjoint_op == Teuchos::null) return Teuchos::null;
192 
193  return Thyra::nonconstMultiVectorLinearOp(adjoint_op, num_adjoint_);
194 }
195 
196 template <typename Scalar>
199 {
200  using Teuchos::RCP;
201  using Teuchos::rcp_dynamic_cast;
203 
204  RCP<const LOWSFB> alowsfb = adjoint_solve_model_->get_W_factory();
205  if (alowsfb == Teuchos::null)
206  return Teuchos::null; // adjoint_solve_model_ doesn't support W_factory
207 
208  return Thyra::multiVectorLinearOpWithSolveFactory(alowsfb, residual_space_,
209  adjoint_space_);
210 }
211 
212 template <typename Scalar>
215 {
216  return prototypeInArgs_;
217 }
218 
219 template <typename Scalar>
222 {
223  typedef Thyra::ModelEvaluatorBase MEB;
224  using Teuchos::RCP;
225  using Teuchos::rcp_dynamic_cast;
226 
227  MEB::InArgs<Scalar> me_nominal = model_->getNominalValues();
228  MEB::InArgs<Scalar> nominal = this->createInArgs();
229 
230  const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
231 
232  // Set initial x, x_dot
233  RCP<Thyra::VectorBase<Scalar> > x = Thyra::createMember(*adjoint_space_);
234  Thyra::assign(x.ptr(), zero);
235  nominal.set_x(x);
236 
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);
242  }
243 
244  const int np = model_->Np();
245  for (int i = 0; i < np; ++i) nominal.set_p(i, me_nominal.get_p(i));
246 
247  return nominal;
248 }
249 
250 template <typename Scalar>
253 {
254  return prototypeOutArgs_;
255 }
256 
257 template <typename Scalar>
260  const Thyra::ModelEvaluatorBase::OutArgs<Scalar>& outArgs) const
261 {
262  typedef Thyra::ModelEvaluatorBase MEB;
263  using Teuchos::RCP;
264  using Teuchos::rcp_dynamic_cast;
265 
266  TEMPUS_FUNC_TIME_MONITOR_DIFF(
267  "Tempus::AdjointSensitivityModelEvaluator::evalModel()", TEMPUS_EVAL);
268 
269  // Note: adjoint models compute the transposed W (either explicitly or
270  // implicitly. Thus we need to always call their evalModel() functions
271  // whenever computing the adjoint operators, and subsequent calls to apply()
272  // do not transpose them.
273 
274  // Interpolate forward solution at supplied time, reusing previous
275  // interpolation if possible
276  TEUCHOS_ASSERT(sh_ != Teuchos::null);
277  const Scalar t = inArgs.get_t();
278  Scalar forward_t;
279  if (is_pseudotransient_)
280  forward_t = forward_state_->getTime();
281  else {
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);
286  else
287  sh_->interpolateState(forward_t, forward_state_.get());
288  t_interp_ = t;
289  }
290  }
291 
292  // setup input arguments for model
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)) {
296  if (inArgs.get_x_dot() != Teuchos::null)
297  me_inArgs.set_x_dot(forward_state_->getXDot());
298  else {
299  if (is_pseudotransient_) {
300  // For pseudo-transient, we have to always use the same form of the
301  // residual in order to reuse df/dx, df/dx_dot,..., so we force
302  // the underlying ME to always compute the implicit form with x_dot == 0
303  // if it wasn't provided.
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));
307  }
308  me_inArgs.set_x_dot(my_x_dot_);
309  }
310  else {
311  // clear out xdot if it was set in nominalValues to get to ensure we
312  // get the explicit form
313  me_inArgs.set_x_dot(Teuchos::null);
314  }
315  }
316  }
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));
320 
321  // compute adjoint W == model W
322  // It would be nice to not reevaluate W in the psuedo-transient case, but
323  // it isn't clear how to do this in a clean way. Probably just need to
324  // control that with the nonlinear solver.
325  RCP<Thyra::LinearOpBase<Scalar> > op;
326  if (outArgs.supports(MEB::OUT_ARG_W_op)) op = outArgs.get_W_op();
327  if (op != Teuchos::null) {
328  TEMPUS_FUNC_TIME_MONITOR_DIFF(
329  "Tempus::AdjointSensitivityModelEvaluator::evalModel::W",
330  TEMPUS_EVAL_W);
331  if (me_inArgs.supports(MEB::IN_ARG_alpha))
332  me_inArgs.set_alpha(inArgs.get_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()); // Implicit form (see below)
336  else
337  me_inArgs.set_beta(-inArgs.get_beta()); // Explicit form (see below)
338  }
339 
340  RCP<Thyra::MultiVectorLinearOp<Scalar> > mv_adjoint_op =
341  rcp_dynamic_cast<Thyra::MultiVectorLinearOp<Scalar> >(op, true);
342  RCP<Thyra::LinearOpBase<Scalar> > adjoint_op =
343  mv_adjoint_op->getNonconstLinearOp();
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);
347  }
348 
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();
356  adjoint_x_mv =
357  rcp_dynamic_cast<const DMVPV>(adjoint_x, true)->getMultiVector();
358  }
359 
360  // Compute adjoint residual F(y):
361  // * For implicit form, F(y) = d/dt( df/dx_dot^T*y ) + df/dx^T*y - dg/dx^T
362  // * For explict form, F(y) = -df/dx^T*y + dg/dx^T
363  // For implicit form, we assume df/dx_dot is constant w.r.t. x, x_dot, and t,
364  // so the residual becomes F(y) = df/dx_dot^T*y_dot + df/dx^T*y - dg/dx^T
365  if (adjoint_f != Teuchos::null) {
366  TEMPUS_FUNC_TIME_MONITOR_DIFF(
367  "Tempus::AdjointSensitivityModelEvaluator::evalModel::f",
368  TEMPUS_EVAL_F);
369 
370  RCP<Thyra::MultiVectorBase<Scalar> > adjoint_f_mv =
371  rcp_dynamic_cast<DMVPV>(adjoint_f, true)->getNonconstMultiVector();
372 
373  MEB::OutArgs<Scalar> me_outArgs = model_->createOutArgs();
374  MEB::OutArgs<Scalar> adj_me_outArgs =
375  adjoint_residual_model_->createOutArgs();
376 
377  // dg/dx^T
378  // Don't re-evaluate dg/dx for pseudotransient
379  {
380  TEMPUS_FUNC_TIME_MONITOR_DIFF(
381  "Tempus::AdjointSensitivityModelEvaluator::evalModel::dg/dx",
382  TEMPUS_EVAL_DGDX);
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_) {
387  me_outArgs.set_DgDx(
388  g_index_,
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;
393  }
394  Thyra::assign(adjoint_f_mv.ptr(), *my_dgdx_mv_);
395  }
396 
397  // Explicit form of the residual F(y) = -df/dx^T*y + dg/dx^T
398  // Don't re-evaluate df/dx for pseudotransient
399  {
400  TEMPUS_FUNC_TIME_MONITOR_DIFF(
401  "Tempus::AdjointSensitivityModelEvaluator::evalModel::df/dx",
402  TEMPUS_EVAL_DFDX);
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;
411  }
412  my_dfdx_->apply(Thyra::NOTRANS, *adjoint_x_mv, adjoint_f_mv.ptr(),
413  Scalar(-1.0), Scalar(1.0));
414  }
415 
416  // Implicit form residual F(y) df/dx_dot^T*y_dot + df/dx^T*y - dg/dx^T
417  // using the second scalar argument to apply() to change the explicit term
418  // above.
419  // Don't re-evaluate df/dx_dot for pseudotransient
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)
428  ->getMultiVector();
429  if (mass_matrix_is_identity_) {
430  // F = -F + y_dot
431  Thyra::V_StVpV(adjoint_f_mv.ptr(), Scalar(-1.0), *adjoint_f_mv,
432  *adjoint_x_dot_mv);
433  }
434  else {
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;
444  }
445  my_dfdxdot_->apply(Thyra::NOTRANS, *adjoint_x_dot_mv,
446  adjoint_f_mv.ptr(), Scalar(1.0), Scalar(-1.0));
447  }
448  }
449  }
450  }
451 
452  // Compute g = dg/dp^T - df/dp^T*y for computing the model parameter term in
453  // the adjoint sensitivity formula.
454  // We don't add pseudotransient logic here because this part is only
455  // evaluated once in that case anyway.
456  if (adjoint_g != Teuchos::null) {
457  TEMPUS_FUNC_TIME_MONITOR_DIFF(
458  "Tempus::AdjointSensitivityModelEvaluator::evalModel::g",
459  TEMPUS_EVAL_G);
460  RCP<Thyra::MultiVectorBase<Scalar> > adjoint_g_mv =
461  rcp_dynamic_cast<DMVPV>(adjoint_g, true)->getNonconstMultiVector();
462 
463  MEB::OutArgs<Scalar> me_outArgs = model_->createOutArgs();
464 
465  // dg/dp
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)) {
469  me_outArgs.set_DgDp(
470  g_index_, p_index_,
471  MEB::Derivative<Scalar>(adjoint_g_mv, MEB::DERIV_MV_GRADIENT_FORM));
472  model_->evalModel(me_inArgs, me_outArgs);
473  }
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);
479  me_outArgs.set_DgDp(
480  g_index_, p_index_,
481  MEB::Derivative<Scalar>(dgdp_trans, MEB::DERIV_MV_JACOBIAN_FORM));
482  model_->evalModel(me_inArgs, me_outArgs);
483  Thyra::DetachedMultiVectorView<Scalar> dgdp_view(*adjoint_g_mv);
484  Thyra::DetachedMultiVectorView<Scalar> dgdp_trans_view(*dgdp_trans);
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);
487  }
488  else
489  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
490  "Invalid dg/dp support");
491  me_outArgs.set_DgDp(g_index_, p_index_, MEB::Derivative<Scalar>());
492 
493  // dg/dp - df/dp^T*y
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_));
501  trans = Thyra::CONJTRANS;
502  }
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());
507  me_outArgs.set_DfDp(
508  p_index_,
509  MEB::Derivative<Scalar>(my_dfdp_mv_, MEB::DERIV_MV_JACOBIAN_FORM));
510  my_dfdp_op_ = my_dfdp_mv_;
511  trans = Thyra::CONJTRANS;
512  }
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());
517  me_outArgs.set_DfDp(
518  p_index_,
519  MEB::Derivative<Scalar>(my_dfdp_mv_, MEB::DERIV_MV_GRADIENT_FORM));
520  my_dfdp_op_ = my_dfdp_mv_;
521  trans = Thyra::CONJ;
522  }
523  else
524  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
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),
528  Scalar(1.0));
529  }
530 
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);
535  }
536 }
537 
538 template <class Scalar>
541 {
542  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
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);
547  return pl;
548 }
549 
550 } // namespace Tempus
551 
552 #endif
static Teuchos::RCP< const Teuchos::ParameterList > getValidParameters()
Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > adjoint_solve_model_
EOpTransp
RCP< const VectorBase< Scalar > > get_x_dot() const
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
bool supports(EOutArgsMembers arg) const
Teuchos::RCP< const Thyra::LinearOpWithSolveFactoryBase< Scalar > > get_W_factory() const
static magnitudeType rmax()
void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
Thyra::ModelEvaluatorBase::OutArgs< Scalar > createOutArgsImpl() const
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