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>
30  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > & model,
31  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > & adjoint_residual_model,
32  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > & adjoint_solve_model,
33  const Scalar& t_init,
34  const Scalar& t_final,
35  const bool is_pseudotransient,
37  model_(model),
38  adjoint_residual_model_(adjoint_residual_model),
39  adjoint_solve_model_(adjoint_solve_model),
40  t_init_(t_init),
41  t_final_(t_final),
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())
47 {
48  typedef Thyra::ModelEvaluatorBase MEB;
49 
50  // Set parameters
53  if (pList != Teuchos::null)
54  *pl = *pList;
56  mass_matrix_is_constant_ = pl->get<bool>("Mass Matrix Is Constant");
57  mass_matrix_is_identity_ = pl->get<bool>("Mass Matrix Is Identity");
58  p_index_ = pl->get<int>("Sensitivity Parameter Index", 0);
59  g_index_ = pl->get<int>("Response Function Index", 0);
60  num_adjoint_ = model_->get_g_space(g_index_)->dim();
61 
62  // We currently do not support a non-constant mass matrix
64  mass_matrix_is_constant_ == false, std::logic_error,
65  "AdjointSensitivityModelEvaluator currently does not support " <<
66  "non-constant mass matrix df/dx_dot!");
67 
69  Thyra::multiVectorProductVectorSpace(model_->get_f_space(), num_adjoint_);
71  Thyra::multiVectorProductVectorSpace(model_->get_x_space(), num_adjoint_);
73  Thyra::multiVectorProductVectorSpace(model_->get_p_space(p_index_),
74  num_adjoint_);
75 
76  // forward and adjoint models must support same InArgs
77  MEB::InArgs<Scalar> me_inArgs = model_->createInArgs();
78  me_inArgs.assertSameSupport(adjoint_residual_model_->createInArgs());
79  me_inArgs.assertSameSupport(adjoint_solve_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  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);
91 
92  // Support additional parameters for x and xdot
93  inArgs.set_Np(me_inArgs.Np());
94  prototypeInArgs_ = inArgs;
95 
96  MEB::OutArgs<Scalar> me_outArgs = model_->createOutArgs();
97  MEB::OutArgs<Scalar> adj_mer_outArgs = adjoint_residual_model_->createOutArgs();
98  MEB::OutArgs<Scalar> adj_mes_outArgs = adjoint_solve_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);
105  prototypeOutArgs_ = outArgs;
106 
107  // Adjoint residual ME must support W_op to define adjoint ODE/DAE.
108  // Must support alpha, beta if it suports x_dot
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));
114  }
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));
121 }
122 
123 template <typename Scalar>
124 void
126 setFinalTime(const Scalar t_final)
127 {
128  t_final_ = t_final;
129 }
130 
131 template <typename Scalar>
132 void
136 {
137  sh_ = sh;
138  if (is_pseudotransient_)
139  forward_state_ = sh_->getCurrentState();
140  else {
142  forward_state_ = Teuchos::null;
143  }
144 
145  // Reset computation flags because we have done a new forward integration
146  mass_matrix_is_computed_ = false;
147  jacobian_matrix_is_computed_ = false;
148  response_gradient_is_computed_ = false;
149 }
150 
151 template <typename Scalar>
154 get_p_space(int p) const
155 {
156  TEUCHOS_ASSERT(p < model_->Np());
157  return model_->get_p_space(p);
158 }
159 
160 template <typename Scalar>
163 get_p_names(int p) const
164 {
165  TEUCHOS_ASSERT(p < model_->Np());
166  return model_->get_p_names(p);
167 }
168 
169 template <typename Scalar>
172 get_x_space() const
173 {
174  return adjoint_space_;
175 }
176 
177 template <typename Scalar>
180 get_f_space() const
181 {
182  return residual_space_;
183 }
184 
185 template <typename Scalar>
188 get_g_space(int j) const
189 {
190  TEUCHOS_ASSERT(j == 0 || j == 1);
191  if (j == 0)
192  return response_space_;
193  return model_->get_g_space(g_index_);
194 }
195 
196 template <typename Scalar>
199 create_W_op() const
200 {
202  adjoint_solve_model_->create_W_op();
203  if (adjoint_op == Teuchos::null)
204  return Teuchos::null;
205 
206  return Thyra::nonconstMultiVectorLinearOp(adjoint_op, num_adjoint_);
207 }
208 
209 template <typename Scalar>
213 {
214  using Teuchos::RCP;
215  using Teuchos::rcp_dynamic_cast;
217 
218  RCP<const LOWSFB> alowsfb = adjoint_solve_model_->get_W_factory();
219  if (alowsfb == Teuchos::null)
220  return Teuchos::null; // adjoint_solve_model_ doesn't support W_factory
221 
222  return Thyra::multiVectorLinearOpWithSolveFactory(
223  alowsfb, residual_space_, adjoint_space_);
224 }
225 
226 template <typename Scalar>
230 {
231  return prototypeInArgs_;
232 }
233 
234 template <typename Scalar>
238 {
239  typedef Thyra::ModelEvaluatorBase MEB;
240  using Teuchos::RCP;
241  using Teuchos::rcp_dynamic_cast;
242 
243  MEB::InArgs<Scalar> me_nominal = model_->getNominalValues();
244  MEB::InArgs<Scalar> nominal = this->createInArgs();
245 
246  const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
247 
248  // Set initial x, x_dot
249  RCP< Thyra::VectorBase<Scalar> > x = Thyra::createMember(*adjoint_space_);
250  Thyra::assign(x.ptr(), zero);
251  nominal.set_x(x);
252 
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);
258  }
259 
260  const int np = model_->Np();
261  for (int i=0; i<np; ++i)
262  nominal.set_p(i, me_nominal.get_p(i));
263 
264  return nominal;
265 }
266 
267 template <typename Scalar>
271 {
272  return prototypeOutArgs_;
273 }
274 
275 template <typename Scalar>
276 void
279  const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
280 {
281  typedef Thyra::ModelEvaluatorBase MEB;
282  using Teuchos::RCP;
283  using Teuchos::rcp_dynamic_cast;
284 
285  TEMPUS_FUNC_TIME_MONITOR_DIFF("Tempus::AdjointSensitivityModelEvaluator::evalModel()", TEMPUS_EVAL);
286 
287  // Note: adjoint models compute the transposed W (either explicitly or
288  // implicitly. Thus we need to always call their evalModel() functions
289  // whenever computing the adjoint operators, and subsequent calls to apply()
290  // do not transpose them.
291 
292  // Interpolate forward solution at supplied time, reusing previous
293  // interpolation if possible
294  TEUCHOS_ASSERT(sh_ != Teuchos::null);
295  const Scalar t = inArgs.get_t();
296  Scalar forward_t;
297  if (is_pseudotransient_)
298  forward_t = forward_state_->getTime();
299  else {
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);
304  else
305  sh_->interpolateState(forward_t, forward_state_.get());
306  t_interp_ = t;
307  }
308  }
309 
310  // setup input arguments for model
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());
316  else {
317  if (is_pseudotransient_) {
318  // For pseudo-transient, we have to always use the same form of the
319  // residual in order to reuse df/dx, df/dx_dot,..., so we force
320  // the underlying ME to always compute the implicit form with x_dot == 0
321  // if it wasn't provided.
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));
325  }
326  me_inArgs.set_x_dot(my_x_dot_);
327  }
328  else {
329  // clear out xdot if it was set in nominalValues to get to ensure we
330  // get the explicit form
331  me_inArgs.set_x_dot(Teuchos::null);
332  }
333  }
334  }
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));
340 
341  // compute adjoint W == model W
342  // It would be nice to not reevaluate W in the psuedo-transient case, but
343  // it isn't clear how to do this in a clean way. Probably just need to
344  // control that with the nonlinear solver.
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()); // Implicit form (see below)
355  else
356  me_inArgs.set_beta(-inArgs.get_beta()); // Explicit form (see below)
357  }
358 
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);
366  }
367 
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) {
374  adjoint_x =
375  inArgs.get_x().assert_not_null();
376  adjoint_x_mv =
377  rcp_dynamic_cast<const DMVPV>(adjoint_x,true)->getMultiVector();
378  }
379 
380  // Compute adjoint residual F(y):
381  // * For implicit form, F(y) = d/dt( df/dx_dot^T*y ) + df/dx^T*y - dg/dx^T
382  // * For explict form, F(y) = -df/dx^T*y + dg/dx^T
383  // For implicit form, we assume df/dx_dot is constant w.r.t. x, x_dot, and t,
384  // so the residual becomes F(y) = df/dx_dot^T*y_dot + df/dx^T*y - dg/dx^T
385  if (adjoint_f != Teuchos::null) {
386  TEMPUS_FUNC_TIME_MONITOR_DIFF("Tempus::AdjointSensitivityModelEvaluator::evalModel::f", TEMPUS_EVAL_F);
387 
388  RCP<Thyra::MultiVectorBase<Scalar> > adjoint_f_mv =
389  rcp_dynamic_cast<DMVPV>(adjoint_f,true)->getNonconstMultiVector();
390 
391  MEB::OutArgs<Scalar> me_outArgs = model_->createOutArgs();
392  MEB::OutArgs<Scalar> adj_me_outArgs = adjoint_residual_model_->createOutArgs();
393 
394  // dg/dx^T
395  // Don't re-evaluate dg/dx for pseudotransient
396  {
397  TEMPUS_FUNC_TIME_MONITOR_DIFF("Tempus::AdjointSensitivityModelEvaluator::evalModel::dg/dx", TEMPUS_EVAL_DGDX);
398  if (my_dgdx_mv_ == Teuchos::null)
399  my_dgdx_mv_ =
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;
410  }
411  Thyra::assign(adjoint_f_mv.ptr(), *my_dgdx_mv_);
412  }
413 
414  // Explicit form of the residual F(y) = -df/dx^T*y + dg/dx^T
415  // Don't re-evaluate df/dx for pseudotransient
416  {
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;
429  }
430  my_dfdx_->apply(Thyra::NOTRANS, *adjoint_x_mv, adjoint_f_mv.ptr(),
431  Scalar(-1.0), Scalar(1.0));
432  }
433 
434  // Implicit form residual F(y) df/dx_dot^T*y_dot + df/dx^T*y - dg/dx^T
435  // using the second scalar argument to apply() to change the explicit term
436  // above.
437  // Don't re-evaluate df/dx_dot for pseudotransient
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_) {
445  // F = -F + y_dot
446  Thyra::V_StVpV(adjoint_f_mv.ptr(), Scalar(-1.0), *adjoint_f_mv,
447  *adjoint_x_dot_mv);
448  }
449  else {
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;
459  }
460  my_dfdxdot_->apply(Thyra::NOTRANS, *adjoint_x_dot_mv,
461  adjoint_f_mv.ptr(), Scalar(1.0), Scalar(-1.0));
462  }
463  }
464  }
465  }
466 
467  // Compute g = dg/dp^T - df/dp^T*y for computing the model parameter term in
468  // the adjoint sensitivity formula.
469  // We don't add pseudotransient logic here because this part is only
470  // evaluated once in that case anyway.
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();
475 
476  MEB::OutArgs<Scalar> me_outArgs = model_->createOutArgs();
477 
478  // dg/dp
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);
486  }
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);
501  }
502  else
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>());
506 
507  // dg/dp - df/dp^T*y
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_));
515  trans = Thyra::CONJTRANS;
516  }
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;
526  }
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_;
535  trans = Thyra::CONJ;
536  }
537  else
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));
543  }
544 
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);
549  }
550 }
551 
552 template<class Scalar>
556 {
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);
562  return pl;
563 }
564 
565 } // namespace Tempus
566 
567 #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