Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tempus_StaggeredForwardSensitivityModelEvaluator_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_StaggeredForwardSensitivityModelEvaluator_impl_hpp
10 #define Tempus_StaggeredForwardSensitivityModelEvaluator_impl_hpp
11 
12 #include "Thyra_DefaultMultiVectorProductVector.hpp"
16 #include "Thyra_VectorStdOps.hpp"
17 #include "Thyra_MultiVectorStdOps.hpp"
18 
19 namespace Tempus {
20 
21 template <typename Scalar>
24  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > & model,
25  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > & sens_residual_model,
26  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > & sens_solve_model,
27  const bool is_pseudotransient,
29  const Teuchos::RCP<MultiVector>& dxdp_init,
30  const Teuchos::RCP<MultiVector>& dx_dotdp_init,
31  const Teuchos::RCP<MultiVector>& dx_dotdotdp_init) :
32  model_(model),
33  sens_residual_model_(sens_residual_model),
34  sens_solve_model_(sens_solve_model),
35  dxdp_init_(dxdp_init),
36  dx_dotdp_init_(dx_dotdp_init),
37  dx_dotdotdp_init_(dx_dotdotdp_init),
38  p_index_(0),
39  g_index_(-1),
40  x_tangent_index_(1),
41  xdot_tangent_index_(2),
42  xdotdot_tangent_index_(3),
43  use_dfdp_as_tangent_(false),
44  use_dgdp_as_tangent_(false),
45  num_param_(0),
46  num_response_(0),
47  g_offset_(0),
48  is_pseudotransient_(is_pseudotransient),
49  mass_matrix_is_computed_(false),
50  jacobian_matrix_is_computed_(false),
51  acceleration_matrix_is_computed_(false),
52  residual_sensitivity_is_computed_(false),
53  t_interp_(Teuchos::ScalarTraits<Scalar>::rmax())
54 {
55  typedef Thyra::ModelEvaluatorBase MEB;
56 
57  // Set parameters
60  if (pList != Teuchos::null)
61  *pl = *pList;
63  use_dfdp_as_tangent_ = pl->get<bool>("Use DfDp as Tangent");
64  use_dgdp_as_tangent_ = pl->get<bool>("Use DgDp as Tangent");
65  p_index_ = pl->get<int>("Sensitivity Parameter Index");
66  g_index_ = pl->get<int>("Response Function Index", -1);
67  x_tangent_index_ = pl->get<int>("Sensitivity X Tangent Index");
68  xdot_tangent_index_ = pl->get<int>("Sensitivity X-Dot Tangent Index");
69  xdotdot_tangent_index_ = pl->get<int>("Sensitivity X-Dot-Dot Tangent Index");
70 
71  num_param_ = model_->get_p_space(p_index_)->dim();
72  if (g_index_ >= 0) {
73  num_response_ = model_->get_g_space(g_index_)->dim();
74  g_offset_ = 2;
75  }
76  dxdp_space_ =
77  Thyra::multiVectorProductVectorSpace(model_->get_x_space(), num_param_);
78  dfdp_space_ =
79  Thyra::multiVectorProductVectorSpace(model_->get_f_space(), num_param_);
80  if (g_index_ >= 0)
81  dgdp_space_ =
82  Thyra::multiVectorProductVectorSpace(model_->get_g_space(g_index_),
83  num_param_);
84 
85  // forward and sensitivity models must support same InArgs
86  MEB::InArgs<Scalar> me_inArgs = model_->createInArgs();
87  me_inArgs.assertSameSupport(sens_residual_model_->createInArgs());
88  me_inArgs.assertSameSupport(sens_solve_model_->createInArgs());
89 
90  MEB::InArgsSetup<Scalar> inArgs;
91  inArgs.setModelEvalDescription(this->description());
92  inArgs.setSupports(MEB::IN_ARG_x);
93  if (me_inArgs.supports(MEB::IN_ARG_x_dot))
94  inArgs.setSupports(MEB::IN_ARG_x_dot);
95  if (me_inArgs.supports(MEB::IN_ARG_t))
96  inArgs.setSupports(MEB::IN_ARG_t);
97  if (me_inArgs.supports(MEB::IN_ARG_alpha))
98  inArgs.setSupports(MEB::IN_ARG_alpha);
99  if (me_inArgs.supports(MEB::IN_ARG_beta))
100  inArgs.setSupports(MEB::IN_ARG_beta);
101  if (me_inArgs.supports(MEB::IN_ARG_W_x_dot_dot_coeff))
102  inArgs.setSupports(MEB::IN_ARG_W_x_dot_dot_coeff);
103 
104  // Support additional parameters for x and xdot
105  inArgs.set_Np(me_inArgs.Np());
106  prototypeInArgs_ = inArgs;
107 
108  MEB::OutArgs<Scalar> me_outArgs = model_->createOutArgs();
109  MEB::OutArgs<Scalar> sens_mer_outArgs = sens_residual_model_->createOutArgs();
110  MEB::OutArgs<Scalar> sens_mes_outArgs = sens_solve_model_->createOutArgs();
111  MEB::OutArgsSetup<Scalar> outArgs;
112  outArgs.setModelEvalDescription(this->description());
113  outArgs.set_Np_Ng(me_inArgs.Np(),me_outArgs.Ng()+g_offset_);
114  outArgs.setSupports(MEB::OUT_ARG_f);
115  if (sens_mer_outArgs.supports(MEB::OUT_ARG_W_op) ||
116  sens_mes_outArgs.supports(MEB::OUT_ARG_W_op))
117  outArgs.setSupports(MEB::OUT_ARG_W_op);
118 
119  // Response 0 is the reduced dg/dp (no sensitivities supported)
120  // Response 1 is the reduced g (no sensitivities supported)
121  for (int j=0; j<me_outArgs.Ng(); ++j) {
122  outArgs.setSupports(MEB::OUT_ARG_DgDx_dot, j+g_offset_,
123  me_outArgs.supports(MEB::OUT_ARG_DgDx_dot, j));
124  outArgs.setSupports(MEB::OUT_ARG_DgDx, j+g_offset_,
125  me_outArgs.supports(MEB::OUT_ARG_DgDx, j));
126  for (int l=0; l<me_outArgs.Np(); ++l) {
127  outArgs.setSupports(MEB::OUT_ARG_DgDp, j+g_offset_, l,
128  me_outArgs.supports(MEB::OUT_ARG_DgDp, j, l));
129  }
130  }
131  prototypeOutArgs_ = outArgs;
132 
133  // Sensitivity residual ME must support W_op to define adjoint ODE/DAE.
134  // Must support alpha, beta if it suports x_dot
135  TEUCHOS_ASSERT(me_inArgs.supports(MEB::IN_ARG_x));
137  TEUCHOS_ASSERT(sens_mer_outArgs.supports(MEB::OUT_ARG_W_op));
138  if (me_inArgs.supports(MEB::IN_ARG_x_dot)) {
139  TEUCHOS_ASSERT(me_inArgs.supports(MEB::IN_ARG_alpha));
140  TEUCHOS_ASSERT(me_inArgs.supports(MEB::IN_ARG_beta));
141  }
142  TEUCHOS_ASSERT(me_outArgs.supports(MEB::OUT_ARG_DfDp, p_index_).supports(MEB::DERIV_MV_JACOBIAN_FORM));
143 }
144 
145 template <typename Scalar>
146 void
150 {
151  sh_ = sh;
153 }
154 
155 template <typename Scalar>
156 void
160 {
161  sh_ = Teuchos::null;
162  forward_state_ = s;
163 }
164 
165 template <typename Scalar>
168 get_p_space(int p) const
169 {
170  TEUCHOS_ASSERT(p < model_->Np());
171  return model_->get_p_space(p);
172 }
173 
174 template <typename Scalar>
177 get_p_names(int p) const
178 {
179  TEUCHOS_ASSERT(p < model_->Np());
180  return model_->get_p_names(p);
181 }
182 
183 template <typename Scalar>
186 get_x_space() const
187 {
188  return dxdp_space_;
189 }
190 
191 template <typename Scalar>
194 get_f_space() const
195 {
196  return dfdp_space_;
197 }
198 
199 template <typename Scalar>
202 get_g_space(int j) const
203 {
204  if (g_index_ >= 0) {
205  if (j == 0)
206  return dgdp_space_;
207  else if (j == 1)
208  return model_->get_g_space(g_index_);
209  }
210  return model_->get_g_space(j-g_offset_);
211 }
212 
213 template <typename Scalar>
216 get_g_names(int j) const
217 {
218  if (g_index_ >= 0) {
219  if (j == 0) {
220  Teuchos::Array<std::string> names = model_->get_g_names(g_index_);
221  for (int i=0; i<names.size(); ++i)
222  names[i] = names[i] + "_reduced sensitivity";
223  return names();
224  }
225  else if (j == 1)
226  return model_->get_g_names(g_index_);
227  }
228  return model_->get_g_names(j-g_offset_);
229 }
230 
231 template <typename Scalar>
234 create_W_op() const
235 {
237  if (lo_ != Teuchos::null)
238  op = lo_;
239  else
240  op = sens_solve_model_->create_W_op();
241  return Thyra::nonconstMultiVectorLinearOp(op, num_param_);
242 }
243 
244 template <typename Scalar>
247 create_DgDx_dot_op(int j) const
248 {
249  return model_->create_DgDx_dot_op(j-g_offset_);
250 }
251 
252 template <typename Scalar>
255 create_DgDx_op(int j) const
256 {
257  return model_->create_DgDx_op(j-g_offset_);
258 }
259 
260 template <typename Scalar>
263 create_DgDp_op(int j, int l) const
264 {
265  return model_->create_DgDp_op(j-g_offset_,l);
266 }
267 
268 template <typename Scalar>
272 {
273  using Teuchos::RCP;
274  using Teuchos::rcp_dynamic_cast;
275 
278  = sens_solve_model_->get_W_factory();
279  if (model_factory == Teuchos::null)
280  return Teuchos::null; // model_ doesn't support W_factory
281  if (po_ != Teuchos::null) {
282  factory =
283  Thyra::reuseLinearOpWithSolveFactory<Scalar>(model_factory, po_);
284  }
285  else
286  factory = model_factory;
287  return Thyra::multiVectorLinearOpWithSolveFactory(factory,
288  dfdp_space_,
289  dxdp_space_);
290 }
291 
292 template <typename Scalar>
296 {
297  return prototypeInArgs_;
298 }
299 
300 template <typename Scalar>
304 {
305  typedef Thyra::ModelEvaluatorBase MEB;
307  using Teuchos::RCP;
308  using Teuchos::rcp_dynamic_cast;
309 
310  MEB::InArgs<Scalar> me_nominal = model_->getNominalValues();
311  MEB::InArgs<Scalar> nominal = this->createInArgs();
312 
313  const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
314 
315  // Set initial x. If dxdp_init == null, set initial dx/dp = 0
316  RCP< Thyra::VectorBase<Scalar> > x = Thyra::createMember(*dxdp_space_);
317  RCP<DMVPV> dxdp = rcp_dynamic_cast<DMVPV>(x,true);
318  if (dxdp_init_ == Teuchos::null)
319  Thyra::assign(dxdp->getNonconstMultiVector().ptr(), zero);
320  else
321  Thyra::assign(dxdp->getNonconstMultiVector().ptr(), *dxdp_init_);
322  nominal.set_x(x);
323 
324  // Set initial xdot. If dx_dotdp_init == null, set initial dxdot/dp = 0
325  if (me_nominal.supports(MEB::IN_ARG_x_dot)) {
326  RCP< Thyra::VectorBase<Scalar> > xdot = Thyra::createMember(*dxdp_space_);
327  RCP<DMVPV> dxdotdp = rcp_dynamic_cast<DMVPV>(xdot,true);
328  if (dx_dotdp_init_ == Teuchos::null)
329  Thyra::assign(dxdotdp->getNonconstMultiVector().ptr(), zero);
330  else
331  Thyra::assign(dxdotdp->getNonconstMultiVector().ptr(), *dx_dotdp_init_);
332  nominal.set_x_dot(xdot);
333  }
334 
335  // Set initial xdotdot. If dx_dotdotdp_init == null, set initial dxdotdot/dp = 0
336  if (me_nominal.supports(MEB::IN_ARG_x_dot_dot)) {
337  RCP< Thyra::VectorBase<Scalar> > xdotdot =
338  Thyra::createMember(*dxdp_space_);
339  RCP<DMVPV> dxdotdotdp = rcp_dynamic_cast<DMVPV>(xdotdot,true);
340  if (dx_dotdotdp_init_ == Teuchos::null)
341  Thyra::assign(dxdotdotdp->getNonconstMultiVector().ptr(), zero);
342  else
343  Thyra::assign(dxdotdotdp->getNonconstMultiVector().ptr(),
344  *dx_dotdotdp_init_);
345  nominal.set_x_dot_dot(xdotdot);
346  }
347 
348  const int np = model_->Np();
349  for (int i=0; i<np; ++i)
350  nominal.set_p(i, me_nominal.get_p(i));
351  return nominal;
352 }
353 
354 template <typename Scalar>
358 {
359  return prototypeOutArgs_;
360 }
361 
362 template <typename Scalar>
363 void
366  const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
367 {
369  typedef Thyra::ModelEvaluatorBase MEB;
370  using Teuchos::RCP;
371  using Teuchos::rcp_dynamic_cast;
372 
373  // Interpolate forward solution at supplied time, reusing previous
374  // interpolation if possible
375  Scalar forward_t;
376  if (sh_ != Teuchos::null) {
377  forward_t = inArgs.get_t();
378  if (t_interp_ != forward_t) {
379  if (nc_forward_state_ == Teuchos::null)
380  nc_forward_state_ = sh_->interpolateState(forward_t);
381  else
382  sh_->interpolateState(forward_t, nc_forward_state_.get());
383  forward_state_ = nc_forward_state_;
384  t_interp_ = forward_t;
385  }
386  }
387  else {
388  TEUCHOS_ASSERT(forward_state_ != Teuchos::null);
389  forward_t = forward_state_->getTime();
390  }
391 
392  const bool use_tangent = use_dfdp_as_tangent_ || use_dgdp_as_tangent_;
393 
394  // setup input arguments for model
395  RCP< const Thyra::MultiVectorBase<Scalar> > dxdp, dxdotdp, dxdotdotdp;
396  MEB::InArgs<Scalar> me_inArgs = model_->getNominalValues();
397  dxdp = rcp_dynamic_cast<const DMVPV>(inArgs.get_x(),true)->getMultiVector();
398  me_inArgs.set_x(forward_state_->getX());
399  if (use_dfdp_as_tangent_)
400  me_inArgs.set_p(x_tangent_index_, inArgs.get_x());
401  if (me_inArgs.supports(MEB::IN_ARG_x_dot)) {
402  if (inArgs.get_x_dot() != Teuchos::null) {
403  dxdotdp =
404  rcp_dynamic_cast<const DMVPV>(inArgs.get_x_dot(),true)->getMultiVector();
405  me_inArgs.set_x_dot(forward_state_->getXDot());
406  if (use_dfdp_as_tangent_)
407  me_inArgs.set_p(xdot_tangent_index_, inArgs.get_x_dot());
408  }
409  else {
410  // clear out xdot if it was set in nominalValues to get to ensure we
411  // get the explicit form
412  me_inArgs.set_x_dot(Teuchos::null);
413  }
414  }
415  if (me_inArgs.supports(MEB::IN_ARG_x_dot_dot)) {
416  if (inArgs.get_x_dot_dot() != Teuchos::null) {
417  dxdotdotdp =
418  rcp_dynamic_cast<const DMVPV>(inArgs.get_x_dot_dot(),true)->getMultiVector();
419  me_inArgs.set_x_dot_dot(forward_state_->getXDotDot());
420  if (use_dfdp_as_tangent_)
421  me_inArgs.set_p(xdotdot_tangent_index_, inArgs.get_x_dot_dot());
422  }
423  else // clear out xdotdot if it was set in nominalValues
424  me_inArgs.set_x_dot_dot(Teuchos::null);
425  }
426  if (me_inArgs.supports(MEB::IN_ARG_t))
427  me_inArgs.set_t(forward_t);
428  if (me_inArgs.supports(MEB::IN_ARG_alpha))
429  me_inArgs.set_alpha(inArgs.get_alpha());
430  if (me_inArgs.supports(MEB::IN_ARG_beta))
431  me_inArgs.set_beta(inArgs.get_beta());
432  if (me_inArgs.supports(MEB::IN_ARG_W_x_dot_dot_coeff))
433  me_inArgs.set_W_x_dot_dot_coeff(inArgs.get_W_x_dot_dot_coeff());
434 
435  // Set parameters -- be careful to only set ones that were set in our
436  // inArgs to not null out any specified through nominalValues or
437  // dx/dp above
438  const int np = me_inArgs.Np();
439  for (int i=0; i<np; ++i) {
440  if (inArgs.get_p(i) != Teuchos::null)
441  if (!use_tangent ||
442  (use_tangent && i != x_tangent_index_ &&
443  i != xdot_tangent_index_ &&
444  i != xdotdot_tangent_index_ ))
445  me_inArgs.set_p(i, inArgs.get_p(i));
446  }
447 
448  // setup output arguments for model
449  RCP< Thyra::MultiVectorBase<Scalar> > dfdp;
450  MEB::OutArgs<Scalar> me_outArgs = model_->createOutArgs();
451  if (outArgs.get_f() != Teuchos::null) {
452  dfdp = rcp_dynamic_cast<DMVPV>(outArgs.get_f(),true)->getNonconstMultiVector();
453  if (!residual_sensitivity_is_computed_) {
454  if (!use_dfdp_as_tangent_ && is_pseudotransient_) {
455  if (my_dfdp_ == Teuchos::null)
456  my_dfdp_ = Thyra::createMembers(model_->get_f_space(),
457  model_->get_p_space(p_index_)->dim());
458  me_outArgs.set_DfDp(p_index_, my_dfdp_);
459  }
460  else
461  me_outArgs.set_DfDp(p_index_, dfdp);
462  }
463  }
464  if (outArgs.supports(MEB::OUT_ARG_W_op) && lo_ == Teuchos::null &&
465  outArgs.get_W_op() != Teuchos::null &&
466  model_.ptr() == sens_solve_model_.ptr()) {
467  RCP<Thyra::LinearOpBase<Scalar> > op = outArgs.get_W_op();
468  RCP<Thyra::MultiVectorLinearOp<Scalar> > mv_op =
469  rcp_dynamic_cast<Thyra::MultiVectorLinearOp<Scalar> >(op,true);
470  me_outArgs.set_W_op(mv_op->getNonconstLinearOp());
471  }
472  for (int j=g_offset_; j<outArgs.Ng(); ++j) {
473  me_outArgs.set_g(j-g_offset_, outArgs.get_g(j));
474  if (!me_outArgs.supports(MEB::OUT_ARG_DgDx_dot,j-g_offset_).none())
475  me_outArgs.set_DgDx_dot(j-g_offset_, outArgs.get_DgDx_dot(j));
476  if (!me_outArgs.supports(MEB::OUT_ARG_DgDx,j-g_offset_).none())
477  me_outArgs.set_DgDx(j-g_offset_, outArgs.get_DgDx(j));
478  for (int l=0; l<outArgs.Np(); ++l)
479  if (!me_outArgs.supports(MEB::OUT_ARG_DgDp,j-g_offset_,l).none())
480  me_outArgs.set_DgDp(j-g_offset_, l, outArgs.get_DgDp(j,l));
481  }
482  if (g_index_ >= 0 && outArgs.get_g(1) != Teuchos::null)
483  me_outArgs.set_g(g_index_, outArgs.get_g(1));
484 
485  // build residual and jacobian
486  model_->evalModel(me_inArgs, me_outArgs);
487 
488  // Compute W_op separately if we have a separate sensitivity solve ME
489  if (outArgs.supports(MEB::OUT_ARG_W_op) && lo_ == Teuchos::null &&
490  outArgs.get_W_op() != Teuchos::null &&
491  model_.ptr() != sens_solve_model_.ptr()) {
492  MEB::OutArgs<Scalar> sens_me_outArgs = sens_solve_model_->createOutArgs();
493  RCP<Thyra::LinearOpBase<Scalar> > op = outArgs.get_W_op();
494  RCP<Thyra::MultiVectorLinearOp<Scalar> > mv_op =
495  rcp_dynamic_cast<Thyra::MultiVectorLinearOp<Scalar> >(op,true);
496  sens_me_outArgs.set_W_op(mv_op->getNonconstLinearOp());
497  sens_solve_model_->evalModel(me_inArgs, sens_me_outArgs);
498  }
499 
500  // Compute (df/dx) * (dx/dp) + (df/dxdot) * (dxdot/dp) + (df/dxdotdot) * (dxdotdot/dp) + (df/dp)
501  // if the underlying ME doesn't already do this. This requires computing
502  // df/dx, df/dxdot, df/dxdotdot as separate operators.
503  // For pseudo-transient, we would like to reuse these operators, but this is
504  // complicated when steppers use both implicit and explicit forms.
505  if (!use_dfdp_as_tangent_) {
506  if (dfdp != Teuchos::null && is_pseudotransient_) {
507  residual_sensitivity_is_computed_ = true;
508  Thyra::assign(dfdp.ptr(), *my_dfdp_);
509  }
510  if (dxdp != Teuchos::null && dfdp != Teuchos::null) {
511  if (my_dfdx_ == Teuchos::null)
512  my_dfdx_ = sens_residual_model_->create_W_op();
513  if (!jacobian_matrix_is_computed_) {
514  MEB::OutArgs<Scalar> meo = sens_residual_model_->createOutArgs();
515  meo.set_W_op(my_dfdx_);
516  if (me_inArgs.supports(MEB::IN_ARG_alpha))
517  me_inArgs.set_alpha(0.0);
518  if (me_inArgs.supports(MEB::IN_ARG_beta))
519  me_inArgs.set_beta(1.0);
520  if (me_inArgs.supports(MEB::IN_ARG_W_x_dot_dot_coeff))
521  me_inArgs.set_W_x_dot_dot_coeff(0.0);
522  sens_residual_model_->evalModel(me_inArgs, meo);
523  if (is_pseudotransient_)
524  jacobian_matrix_is_computed_ = true;
525  }
526  my_dfdx_->apply(Thyra::NOTRANS, *dxdp, dfdp.ptr(), Scalar(1.0), Scalar(1.0));
527  }
528  if (dxdotdp != Teuchos::null && dfdp != Teuchos::null) {
529  if (my_dfdxdot_ == Teuchos::null)
530  my_dfdxdot_ = sens_residual_model_->create_W_op();
531  if (!mass_matrix_is_computed_) {
532  MEB::OutArgs<Scalar> meo = sens_residual_model_->createOutArgs();
533  meo.set_W_op(my_dfdxdot_);
534  if (me_inArgs.supports(MEB::IN_ARG_alpha))
535  me_inArgs.set_alpha(1.0);
536  if (me_inArgs.supports(MEB::IN_ARG_beta))
537  me_inArgs.set_beta(0.0);
538  if (me_inArgs.supports(MEB::IN_ARG_W_x_dot_dot_coeff))
539  me_inArgs.set_W_x_dot_dot_coeff(0.0);
540  sens_residual_model_->evalModel(me_inArgs, meo);
541  if (is_pseudotransient_)
542  mass_matrix_is_computed_ = true;
543  }
544  my_dfdxdot_->apply(Thyra::NOTRANS, *dxdotdp, dfdp.ptr(), Scalar(1.0), Scalar(1.0));
545  }
546  if (dxdotdotdp != Teuchos::null && dfdp != Teuchos::null) {
547  if (my_dfdxdotdot_ == Teuchos::null)
548  my_dfdxdotdot_ = sens_residual_model_->create_W_op();
549  if (!acceleration_matrix_is_computed_) {
550  MEB::OutArgs<Scalar> meo = sens_residual_model_->createOutArgs();
551  meo.set_W_op(my_dfdxdotdot_);
552  if (me_inArgs.supports(MEB::IN_ARG_alpha))
553  me_inArgs.set_alpha(0.0);
554  if (me_inArgs.supports(MEB::IN_ARG_beta))
555  me_inArgs.set_beta(0.0);
556  if (me_inArgs.supports(MEB::IN_ARG_W_x_dot_dot_coeff))
557  me_inArgs.set_W_x_dot_dot_coeff(1.0);
558  sens_residual_model_->evalModel(me_inArgs, meo);
559  if (is_pseudotransient_)
560  acceleration_matrix_is_computed_ = true;
561  }
562  my_dfdxdotdot_->apply(Thyra::NOTRANS, *dxdotdotdp, dfdp.ptr(), Scalar(1.0), Scalar(1.0));
563  }
564  }
565 
566  if (g_index_ >= 0 && outArgs.get_g(0) != Teuchos::null) {
567  MEB::OutArgs<Scalar> meo = model_->createOutArgs();
568  RCP<Thyra::MultiVectorBase<Scalar> > dgdp =
569  rcp_dynamic_cast<DMVPV>(outArgs.get_g(0),true)->getNonconstMultiVector();
570  RCP<Thyra::MultiVectorBase<Scalar> > dgdp_trans;
571  MEB::DerivativeSupport dgdp_support =
572  meo.supports(MEB::OUT_ARG_DgDp, g_index_, p_index_);
573  if (dgdp_support.supports(MEB::DERIV_MV_JACOBIAN_FORM))
574  meo.set_DgDp(g_index_, p_index_,
575  MEB::Derivative<Scalar>(dgdp, MEB::DERIV_MV_JACOBIAN_FORM));
576  else if (dgdp_support.supports(MEB::DERIV_MV_GRADIENT_FORM)) {
577  dgdp_trans = createMembers(model_->get_p_space(p_index_),
578  num_response_);
579  meo.set_DgDp(g_index_, p_index_,
580  MEB::Derivative<Scalar>(dgdp_trans, MEB::DERIV_MV_GRADIENT_FORM));
581  }
582  else
584  true, std::logic_error,
585  "Operator form of dg/dp not supported for reduced sensitivity");
586 
587  if (!use_dgdp_as_tangent_ || dxdp == Teuchos::null) {
588  MEB::DerivativeSupport dgdx_support =
589  meo.supports(MEB::OUT_ARG_DgDx, g_index_);
590  if (dgdx_support.supports(MEB::DERIV_LINEAR_OP)) {
591  if (my_dgdx_ == Teuchos::null)
592  my_dgdx_ = model_->create_DgDx_op(g_index_);
593  meo.set_DgDx(g_index_, my_dgdx_);
594  }
595  else if (dgdx_support.supports(MEB::DERIV_MV_GRADIENT_FORM)) {
596  if (my_dgdx_mv_ == Teuchos::null)
597  my_dgdx_mv_ = createMembers(model_->get_g_space(g_index_), num_param_);
598  meo.set_DgDx(g_index_,
599  MEB::Derivative<Scalar>(my_dgdx_mv_,
600  MEB::DERIV_MV_GRADIENT_FORM));
601  }
602  else
604  true, std::logic_error,
605  "Jacobian form of dg/dx not supported for reduced sensitivity");
606  }
607 
608  // Clear dx/dp, dxdot/dp, dxdotdot/dp from inArgs if set and we are not
609  // using dg/dp as the tangent. Note, even though dxdp, dxdotdp, and
610  // dxdotdotdp are not used here, they are non-null if x/xdot/xdotdot are
611  // supported and supplied.
612  if (!use_dgdp_as_tangent_ && use_dfdp_as_tangent_ &&
613  dxdp != Teuchos::null)
614  me_inArgs.set_p(x_tangent_index_, Teuchos::null);
615  if (!use_dgdp_as_tangent_ && use_dfdp_as_tangent_ &&
616  dxdotdp != Teuchos::null)
617  me_inArgs.set_p(xdot_tangent_index_, Teuchos::null);
618  if (!use_dgdp_as_tangent_ && use_dfdp_as_tangent_ &&
619  dxdotdotdp != Teuchos::null)
620  me_inArgs.set_p(xdotdot_tangent_index_, Teuchos::null);
621 
622  // Set dx/dp, dxdot/dp, dxdodot/dp if necessary
623  if (use_dgdp_as_tangent_ && !use_dfdp_as_tangent_ &&
624  dxdp != Teuchos::null)
625  me_inArgs.set_p(x_tangent_index_, inArgs.get_x());
626  if (use_dgdp_as_tangent_ && !use_dfdp_as_tangent_ &&
627  dxdotdp != Teuchos::null)
628  me_inArgs.set_p(xdot_tangent_index_, inArgs.get_x_dot());
629  if (use_dgdp_as_tangent_ && !use_dfdp_as_tangent_ &&
630  dxdotdotdp != Teuchos::null)
631  me_inArgs.set_p(xdotdot_tangent_index_, inArgs.get_x_dot_dot());
632 
633  model_->evalModel(me_inArgs, meo);
634 
635  // transpose reduced dg/dp if necessary
636  if (dgdp_trans != Teuchos::null) {
638  Thyra::DetachedMultiVectorView<Scalar> dgdp_trans_view(*dgdp_trans);
639  for (int i=0; i<num_param_; ++i)
640  for (int j=0; j<num_response_; ++j)
641  dgdp_view(j,i) = dgdp_trans_view(i,j);
642  }
643 
644  // Compute (dg/dx) * (dx/dp) + (dg/dp) if the underlying ME doesn't already
645  // do this.
646  if (!use_dgdp_as_tangent_ && dxdp != Teuchos::null) {
647  MEB::DerivativeSupport dgdx_support =
648  me_outArgs.supports(MEB::OUT_ARG_DgDx, g_index_);
649  if (dgdx_support.supports(MEB::DERIV_LINEAR_OP)) {
650  my_dgdx_->apply(Thyra::NOTRANS, *dxdp, dgdp.ptr(), Scalar(1.0), Scalar(1.0));
651  }
652  else if (dgdx_support.supports(MEB::DERIV_MV_GRADIENT_FORM)) {
653  my_dgdx_mv_->apply(Thyra::TRANS, *dxdp, dgdp.ptr(), Scalar(1.0), Scalar(1.0));
654  }
655  else
657  true, std::logic_error,
658  "Jacobian form of dg/dx not supported for reduced sensitivity");
659  }
660  }
661 }
662 
663 template<class Scalar>
667 {
668  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
669  pl->set<bool>("Use DfDp as Tangent", false);
670  pl->set<bool>("Use DgDp as Tangent", false);
671  pl->set<int>("Sensitivity Parameter Index", 0);
672  pl->set<int>("Response Function Index", -1);
673  pl->set<int>("Sensitivity X Tangent Index", 1);
674  pl->set<int>("Sensitivity X-Dot Tangent Index", 2);
675  pl->set<int>("Sensitivity X-Dot-Dot Tangent Index", 3);
676  return pl;
677 }
678 
679 } // namespace Tempus
680 
681 #endif
RCP< const VectorBase< Scalar > > get_x_dot() const
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
void setForwardSolutionHistory(const Teuchos::RCP< const Tempus::SolutionHistory< Scalar > > &sh)
Set solution history from forward state evaluation (for interpolation)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
virtual void setForwardSolutionState(const Teuchos::RCP< const Tempus::SolutionState< Scalar > > &s)
Set solution state from forward state evaluation (for frozen state)
Evaluation< VectorBase< Scalar > > get_g(int j) const
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int p) const
Evaluation< VectorBase< Scalar > > get_f() const
bool supports(EOutArgsMembers arg) const
static magnitudeType rmax()
Derivative< Scalar > get_DgDp(int j, int l) const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_p_space(int p) const
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_DgDx_dot_op(int j) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual std::string description() const
void validateParametersAndSetDefaults(ParameterList const &validParamList, int const depth=1000)
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_DgDp_op(int j, int l) const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
Implicit concrete LinearOpBase subclass that takes a flattended out multi-vector and performs a multi...
size_type size() const
Teuchos::RCP< const Thyra::LinearOpWithSolveFactoryBase< Scalar > > get_W_factory() const
Derivative< Scalar > get_DgDx_dot(int j) const
RCP< const VectorBase< Scalar > > get_x_dot_dot() const
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_DgDx_op(int j) const
#define TEUCHOS_ASSERT(assertion_test)
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_f_space() const
RCP< LinearOpBase< Scalar > > get_W_op() const
RCP< const VectorBase< Scalar > > get_x() const
StaggeredForwardSensitivityModelEvaluator(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &sens_residual_model, const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &sens_solve_model, const bool is_pseudotransient, const Teuchos::RCP< const Teuchos::ParameterList > &pList=Teuchos::null, const Teuchos::RCP< MultiVector > &dxdp_init=Teuchos::null, const Teuchos::RCP< MultiVector > &dx_dotdp_init=Teuchos::null, const Teuchos::RCP< MultiVector > &dx_dotdot_dp_init=Teuchos::null)
Constructor.
Solution state for integrators and steppers.
RCP< const VectorBase< Scalar > > get_p(int l) const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_g_space(int j) const
void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
Derivative< Scalar > get_DgDx(int j) const