Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tempus_CombinedForwardSensitivityModelEvaluator_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_CombinedForwardSensitivityModelEvaluator_impl_hpp
10 #define Tempus_CombinedForwardSensitivityModelEvaluator_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,
28  const Teuchos::RCP<MultiVector>& dxdp_init,
29  const Teuchos::RCP<MultiVector>& dx_dotdp_init,
30  const Teuchos::RCP<MultiVector>& dx_dotdotdp_init) :
31  model_(model),
32  sens_residual_model_(sens_residual_model),
33  sens_solve_model_(sens_solve_model),
34  dxdp_init_(dxdp_init),
35  dx_dotdp_init_(dx_dotdp_init),
36  dx_dotdotdp_init_(dx_dotdotdp_init),
37  p_index_(0),
38  g_index_(-1),
39  x_tangent_index_(1),
40  xdot_tangent_index_(2),
41  xdotdot_tangent_index_(3),
42  use_dfdp_as_tangent_(false),
43  use_dgdp_as_tangent_(false),
44  num_param_(0),
45  num_response_(0),
46  g_offset_(0)
47 {
48  typedef Thyra::ModelEvaluatorBase MEB;
49 
50  // Set parameters
53  if (pList != Teuchos::null)
54  *pl = *pList;
56  use_dfdp_as_tangent_ = pl->get<bool>("Use DfDp as Tangent");
57  use_dgdp_as_tangent_ = pl->get<bool>("Use DgDp as Tangent");
58  p_index_ = pl->get<int>("Sensitivity Parameter Index");
59  g_index_ = pl->get<int>("Response Function Index");
60  x_tangent_index_ = pl->get<int>("Sensitivity X Tangent Index");
61  xdot_tangent_index_ = pl->get<int>("Sensitivity X-Dot Tangent Index");
62  xdotdot_tangent_index_ = pl->get<int>("Sensitivity X-Dot-Dot Tangent Index");
63 
64  num_param_ = model_->get_p_space(p_index_)->dim();
65  if (g_index_ >= 0) {
66  num_response_ = model_->get_g_space(g_index_)->dim();
67  g_offset_ = 2;
68  }
69  dxdp_space_ =
70  Thyra::multiVectorProductVectorSpace(model_->get_x_space(), num_param_);
72  Thyra::multiVectorProductVectorSpace(model_->get_x_space(), num_param_+1);
73  dfdp_space_ =
74  Thyra::multiVectorProductVectorSpace(model_->get_f_space(), num_param_);
76  Thyra::multiVectorProductVectorSpace(model_->get_f_space(), num_param_+1);
77  if (g_index_ >= 0)
78  dgdp_space_ =
79  Thyra::multiVectorProductVectorSpace(model_->get_g_space(g_index_),
80  num_param_);
81 
82  // forward and sensitivity models must support same InArgs
83  MEB::InArgs<Scalar> me_inArgs = model_->createInArgs();
84  me_inArgs.assertSameSupport(sens_residual_model_->createInArgs());
85  me_inArgs.assertSameSupport(sens_solve_model_->createInArgs());
86 
87  MEB::InArgsSetup<Scalar> inArgs;
88  inArgs.setModelEvalDescription(this->description());
89  inArgs.setSupports(MEB::IN_ARG_x);
90  if (me_inArgs.supports(MEB::IN_ARG_x_dot))
91  inArgs.setSupports(MEB::IN_ARG_x_dot);
92  if (me_inArgs.supports(MEB::IN_ARG_t))
93  inArgs.setSupports(MEB::IN_ARG_t);
94  if (me_inArgs.supports(MEB::IN_ARG_alpha))
95  inArgs.setSupports(MEB::IN_ARG_alpha);
96  if (me_inArgs.supports(MEB::IN_ARG_beta))
97  inArgs.setSupports(MEB::IN_ARG_beta);
98  if (me_inArgs.supports(MEB::IN_ARG_W_x_dot_dot_coeff))
99  inArgs.setSupports(MEB::IN_ARG_W_x_dot_dot_coeff);
100 
101  // Support additional parameters for x and xdot
102  inArgs.set_Np(me_inArgs.Np());
103  prototypeInArgs_ = inArgs;
104 
105  MEB::OutArgs<Scalar> me_outArgs = model_->createOutArgs();
106  MEB::OutArgs<Scalar> sens_mer_outArgs = sens_residual_model_->createOutArgs();
107  MEB::OutArgs<Scalar> sens_mes_outArgs = sens_solve_model_->createOutArgs();
108  MEB::OutArgsSetup<Scalar> outArgs;
109  outArgs.setModelEvalDescription(this->description());
110  outArgs.set_Np_Ng(me_outArgs.Np(),me_outArgs.Ng()+g_offset_);
111  outArgs.setSupports(MEB::OUT_ARG_f);
112  if (sens_mer_outArgs.supports(MEB::OUT_ARG_W_op) ||
113  sens_mes_outArgs.supports(MEB::OUT_ARG_W_op))
114  outArgs.setSupports(MEB::OUT_ARG_W_op);
115 
116  // Response 0 is the reduced dg/dp (no sensitivities supported)
117  // Response 1 is the reduced g (no sensitivities supported)
118  for (int j=0; j<me_outArgs.Ng(); ++j) {
119  outArgs.setSupports(MEB::OUT_ARG_DgDx_dot, j+g_offset_,
120  me_outArgs.supports(MEB::OUT_ARG_DgDx_dot, j));
121  outArgs.setSupports(MEB::OUT_ARG_DgDx, j+g_offset_,
122  me_outArgs.supports(MEB::OUT_ARG_DgDx, j));
123  for (int l=0; l<me_outArgs.Np(); ++l) {
124  outArgs.setSupports(MEB::OUT_ARG_DgDp, j+g_offset_, l,
125  me_outArgs.supports(MEB::OUT_ARG_DgDp, j, l));
126  }
127  }
128  prototypeOutArgs_ = outArgs;
129 
130  // Sensitivity residual ME must support W_op to define adjoint ODE/DAE.
131  // Must support alpha, beta if it suports x_dot
132  TEUCHOS_ASSERT(me_inArgs.supports(MEB::IN_ARG_x));
134  TEUCHOS_ASSERT(sens_mer_outArgs.supports(MEB::OUT_ARG_W_op));
135  if (me_inArgs.supports(MEB::IN_ARG_x_dot)) {
136  TEUCHOS_ASSERT(me_inArgs.supports(MEB::IN_ARG_alpha));
137  TEUCHOS_ASSERT(me_inArgs.supports(MEB::IN_ARG_beta));
138  }
139  TEUCHOS_ASSERT(me_outArgs.supports(MEB::OUT_ARG_DfDp, p_index_).supports(MEB::DERIV_MV_JACOBIAN_FORM));
140 }
141 
142 template <typename Scalar>
145 get_p_space(int p) const
146 {
147  return model_->get_p_space(p);
148 }
149 
150 template <typename Scalar>
153 get_p_names(int p) const
154 {
155  return model_->get_p_names(p);
156 }
157 
158 template <typename Scalar>
161 get_x_space() const
162 {
163  return x_dxdp_space_;
164 }
165 
166 template <typename Scalar>
169 get_f_space() const
170 {
171  return f_dfdp_space_;
172 }
173 
174 template <typename Scalar>
177 get_g_space(int j) const
178 {
179  if (g_index_ >= 0) {
180  if (j == 0)
181  return dgdp_space_;
182  else if (j == 1)
183  return model_->get_g_space(g_index_);
184  }
185  return model_->get_g_space(j-g_offset_);
186 }
187 
188 template <typename Scalar>
191 get_g_names(int j) const
192 {
193  if (g_index_ >= 0) {
194  if (j == 0) {
195  Teuchos::Array<std::string> names = model_->get_g_names(g_index_);
196  for (int i=0; i<names.size(); ++i)
197  names[i] = names[i] + "_reduced sensitivity";
198  return names();
199  }
200  else if (j == 1)
201  return model_->get_g_names(g_index_);
202  }
203  return model_->get_g_names(j-g_offset_);
204 }
205 
206 template <typename Scalar>
209 create_W_op() const
210 {
212  sens_solve_model_->create_W_op();
213  return Thyra::nonconstMultiVectorLinearOp(op, num_param_+1);
214 }
215 
216 template <typename Scalar>
219 create_DgDx_dot_op(int j) const
220 {
221  return model_->create_DgDx_dot_op(j-g_offset_);
222 }
223 
224 template <typename Scalar>
227 create_DgDx_op(int j) const
228 {
229  return model_->create_DgDx_op(j-g_offset_);
230 }
231 
232 template <typename Scalar>
235 create_DgDp_op(int j, int l) const
236 {
237  return model_->create_DgDp_op(j-g_offset_,l);
238 }
239 
240 template <typename Scalar>
244 {
246  sens_solve_model_->get_W_factory();
247  if (factory == Teuchos::null)
248  return Teuchos::null; // model_ doesn't support W_factory
249 
250  return Thyra::multiVectorLinearOpWithSolveFactory(factory,
251  f_dfdp_space_,
252  x_dxdp_space_);
253 }
254 
255 template <typename Scalar>
259 {
260  return prototypeInArgs_;
261 }
262 
263 template <typename Scalar>
267 {
268  typedef Thyra::ModelEvaluatorBase MEB;
270  using Teuchos::RCP;
271  using Teuchos::rcp_dynamic_cast;
272  using Teuchos::Range1D;
273 
274  MEB::InArgs<Scalar> me_nominal = model_->getNominalValues();
275  MEB::InArgs<Scalar> nominal = this->createInArgs();
276 
277  const Teuchos::Range1D rng(1,num_param_);
278  const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
279 
280  // Set initial x. If dxdp_init == null, set initial dx/dp = 0
281  RCP< const Thyra::VectorBase<Scalar> > me_x = me_nominal.get_x();
282  if (me_x != Teuchos::null) {
283  RCP< Thyra::VectorBase<Scalar> > x = Thyra::createMember(*x_dxdp_space_);
284  RCP<DMVPV> x_dxdp = rcp_dynamic_cast<DMVPV>(x,true);
285  Thyra::assign(x_dxdp->getNonconstMultiVector()->col(0).ptr(), *me_x);
286  if (dxdp_init_ == Teuchos::null)
287  Thyra::assign(x_dxdp->getNonconstMultiVector()->subView(rng).ptr(),
288  zero);
289  else
290  Thyra::assign(x_dxdp->getNonconstMultiVector()->subView(rng).ptr(),
291  *dxdp_init_);
292  nominal.set_x(x);
293  }
294 
295  // Set initial xdot. If dx_dotdp_init == null, set initial dxdot/dp = 0
296  RCP< const Thyra::VectorBase<Scalar> > me_xdot;
297  if (me_nominal.supports(MEB::IN_ARG_x_dot))
298  me_xdot = me_nominal.get_x_dot();
299  if (me_xdot != Teuchos::null) {
300  RCP< Thyra::VectorBase<Scalar> > xdot = Thyra::createMember(*x_dxdp_space_);
301  RCP<DMVPV> xdot_dxdp = rcp_dynamic_cast<DMVPV>(xdot,true);
302  Thyra::assign(xdot_dxdp->getNonconstMultiVector()->col(0).ptr(), *me_xdot);
303  if (dx_dotdp_init_ == Teuchos::null)
304  Thyra::assign(xdot_dxdp->getNonconstMultiVector()->subView(rng).ptr(),
305  zero);
306  else
307  Thyra::assign(xdot_dxdp->getNonconstMultiVector()->subView(rng).ptr(),
308  *dx_dotdp_init_);
309  nominal.set_x_dot(xdot);
310  }
311 
312  // Set initial xdotdot. If dx_dotdotdp_init == null, set initial dxdotdot/dp = 0
313  RCP< const Thyra::VectorBase<Scalar> > me_xdotdot;
314  if (me_nominal.supports(MEB::IN_ARG_x_dot_dot))
315  me_xdotdot = me_nominal.get_x_dot_dot();
316  if (me_xdotdot != Teuchos::null) {
317  RCP< Thyra::VectorBase<Scalar> > xdotdot =
318  Thyra::createMember(*x_dxdp_space_);
319  RCP<DMVPV> xdotdot_dxdp = rcp_dynamic_cast<DMVPV>(xdotdot,true);
320  Thyra::assign(xdotdot_dxdp->getNonconstMultiVector()->col(0).ptr(), *me_xdotdot);
321  if (dx_dotdotdp_init_ == Teuchos::null)
322  Thyra::assign(xdotdot_dxdp->getNonconstMultiVector()->subView(rng).ptr(),
323  zero);
324  else
325  Thyra::assign(xdotdot_dxdp->getNonconstMultiVector()->subView(rng).ptr(),
326  *dx_dotdotdp_init_);
327  nominal.set_x_dot_dot(xdotdot);
328  }
329 
330  const int np = model_->Np();
331  for (int i=0; i<np; ++i)
332  nominal.set_p(i, me_nominal.get_p(i));
333 
334  return nominal;
335 }
336 
337 template <typename Scalar>
341 {
342  return prototypeOutArgs_;
343 }
344 
345 template <typename Scalar>
346 void
349  const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
350 {
351  typedef Thyra::ModelEvaluatorBase MEB;
353  using Teuchos::RCP;
354  using Teuchos::rcp_dynamic_cast;
355  using Teuchos::Range1D;
356 
357  const bool use_tangent = use_dfdp_as_tangent_ || use_dgdp_as_tangent_;
358 
359  // setup input arguments for model
360  RCP< const Thyra::VectorBase<Scalar> > x, xdot, xdotdot;
361  RCP< const Thyra::MultiVectorBase<Scalar> > dxdp, dxdotdp, dxdotdotdp;
362  RCP< const Thyra::VectorBase<Scalar> > dxdp_vec, dxdotdp_vec, dxdotdotdp_vec;
363  MEB::InArgs<Scalar> me_inArgs = model_->getNominalValues();
364  RCP<const DMVPV> x_dxdp = rcp_dynamic_cast<const DMVPV>(inArgs.get_x(),true);
365  x = x_dxdp->getMultiVector()->col(0);
366  dxdp = x_dxdp->getMultiVector()->subView(Range1D(1,num_param_));
367  me_inArgs.set_x(x);
368  if (use_tangent) {
369  dxdp_vec = Thyra::multiVectorProductVector(dxdp_space_, dxdp);
370  if (use_dfdp_as_tangent_)
371  me_inArgs.set_p(x_tangent_index_, dxdp_vec);
372  }
373  if (me_inArgs.supports(MEB::IN_ARG_x_dot)) {
374  if (inArgs.get_x_dot() != Teuchos::null) {
375  RCP<const DMVPV> xdot_dxdotdp =
376  rcp_dynamic_cast<const DMVPV>(inArgs.get_x_dot(),true);
377  xdot = xdot_dxdotdp->getMultiVector()->col(0);
378  dxdotdp = xdot_dxdotdp->getMultiVector()->subView(Range1D(1,num_param_));
379  me_inArgs.set_x_dot(xdot);
380  if (use_tangent) {
381  dxdotdp_vec = Thyra::multiVectorProductVector(dxdp_space_, dxdotdp);
382  if (use_dfdp_as_tangent_)
383  me_inArgs.set_p(xdot_tangent_index_, dxdotdp_vec);
384  }
385  }
386  else // clear out xdot if it was set in nominalValues
387  me_inArgs.set_x_dot(Teuchos::null);
388  }
389  if (me_inArgs.supports(MEB::IN_ARG_x_dot_dot)) {
390  if (inArgs.get_x_dot_dot() != Teuchos::null) {
391  RCP<const DMVPV> xdotdot_dxdotdotdp =
392  rcp_dynamic_cast<const DMVPV>(inArgs.get_x_dot_dot(),true);
393  xdotdot = xdotdot_dxdotdotdp->getMultiVector()->col(0);
394  dxdotdotdp = xdotdot_dxdotdotdp->getMultiVector()->subView(Range1D(1,num_param_));
395  me_inArgs.set_x_dot_dot(xdotdot);
396  if (use_tangent) {
397  dxdotdotdp_vec =
398  Thyra::multiVectorProductVector(dxdp_space_, dxdotdotdp);
399  if (use_dfdp_as_tangent_)
400  me_inArgs.set_p(xdotdot_tangent_index_, dxdotdotdp_vec);
401  }
402  }
403  else // clear out xdotdot if it was set in nominalValues
404  me_inArgs.set_x_dot_dot(Teuchos::null);
405  }
406  if (me_inArgs.supports(MEB::IN_ARG_t))
407  me_inArgs.set_t(inArgs.get_t());
408  if (me_inArgs.supports(MEB::IN_ARG_alpha))
409  me_inArgs.set_alpha(inArgs.get_alpha());
410  if (me_inArgs.supports(MEB::IN_ARG_beta))
411  me_inArgs.set_beta(inArgs.get_beta());
412  if (me_inArgs.supports(MEB::IN_ARG_W_x_dot_dot_coeff))
413  me_inArgs.set_W_x_dot_dot_coeff(inArgs.get_W_x_dot_dot_coeff());
414 
415  // Set parameters -- be careful to only set ones that were set in our
416  // inArgs to not null out any specified through nominalValues or
417  // dx/dp above
418  const int np = me_inArgs.Np();
419  for (int i=0; i<np; ++i) {
420  if (inArgs.get_p(i) != Teuchos::null)
421  if (!use_tangent ||
422  (use_tangent && i != x_tangent_index_ &&
423  i != xdot_tangent_index_ &&
424  i != xdotdot_tangent_index_ ))
425  me_inArgs.set_p(i, inArgs.get_p(i));
426  }
427 
428  // setup output arguments for model
429  RCP< Thyra::VectorBase<Scalar> > f;
430  RCP< Thyra::MultiVectorBase<Scalar> > dfdp;
431  MEB::OutArgs<Scalar> me_outArgs = model_->createOutArgs();
432  if (outArgs.get_f() != Teuchos::null) {
433  RCP<DMVPV> f_dfdp = rcp_dynamic_cast<DMVPV>(outArgs.get_f(),true);
434  f = f_dfdp->getNonconstMultiVector()->col(0);
435  dfdp = f_dfdp->getNonconstMultiVector()->subView(Range1D(1,num_param_));
436  me_outArgs.set_f(f);
437  me_outArgs.set_DfDp(p_index_, dfdp);
438  }
439  if (outArgs.supports(MEB::OUT_ARG_W_op) &&
440  outArgs.get_W_op() != Teuchos::null &&
441  model_.ptr() == sens_solve_model_.ptr()) {
442  RCP< Thyra::LinearOpBase<Scalar> > op = outArgs.get_W_op();
443  RCP< Thyra::MultiVectorLinearOp<Scalar> > mv_op =
444  rcp_dynamic_cast< Thyra::MultiVectorLinearOp<Scalar> >(op,true);
445  me_outArgs.set_W_op(mv_op->getNonconstLinearOp());
446  }
447  for (int j=g_offset_; j<outArgs.Ng(); ++j) {
448  if (!me_outArgs.supports(MEB::OUT_ARG_DgDx_dot,j-g_offset_).none())
449  me_outArgs.set_DgDx_dot(j-g_offset_, outArgs.get_DgDx_dot(j));
450  if (!me_outArgs.supports(MEB::OUT_ARG_DgDx,j-g_offset_).none())
451  me_outArgs.set_DgDx(j-g_offset_, outArgs.get_DgDx(j));
452  for (int l=0; l<outArgs.Np(); ++l)
453  if (!me_outArgs.supports(MEB::OUT_ARG_DgDp,j-g_offset_,l).none())
454  me_outArgs.set_DgDp(j-g_offset_, l, outArgs.get_DgDp(j,l));
455  }
456  if (g_index_ >= 0 && outArgs.get_g(1) != Teuchos::null)
457  me_outArgs.set_g(g_index_, outArgs.get_g(1));
458 
459  // build residual and jacobian
460  model_->evalModel(me_inArgs, me_outArgs);
461 
462  // Compute W_op separately if we have a separate sensitivity solve ME
463  if (outArgs.supports(MEB::OUT_ARG_W_op) &&
464  outArgs.get_W_op() != Teuchos::null &&
465  model_.ptr() != sens_solve_model_.ptr()) {
466  MEB::OutArgs<Scalar> sens_me_outArgs = sens_solve_model_->createOutArgs();
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  sens_me_outArgs.set_W_op(mv_op->getNonconstLinearOp());
471  sens_solve_model_->evalModel(me_inArgs, sens_me_outArgs);
472  }
473 
474  // Compute (df/dx) * (dx/dp) + (df/dxdot) * (dxdot/dp) + (df/dxdotdot) * (dxdotdot/dp) + (df/dp)
475  // if the underlying ME doesn't already do this. This requires computing
476  // df/dx, df/dxdot, df/dxdotdot as separate operators
477  if (!use_dfdp_as_tangent_) {
478  if (dxdp != Teuchos::null && dfdp != Teuchos::null) {
479  if (my_dfdx_ == Teuchos::null)
480  my_dfdx_ = sens_residual_model_->create_W_op();
481  MEB::OutArgs<Scalar> meo = sens_residual_model_->createOutArgs();
482  meo.set_W_op(my_dfdx_);
483  if (me_inArgs.supports(MEB::IN_ARG_alpha))
484  me_inArgs.set_alpha(0.0);
485  if (me_inArgs.supports(MEB::IN_ARG_beta))
486  me_inArgs.set_beta(1.0);
487  if (me_inArgs.supports(MEB::IN_ARG_W_x_dot_dot_coeff))
488  me_inArgs.set_W_x_dot_dot_coeff(0.0);
489  sens_residual_model_->evalModel(me_inArgs, meo);
490  my_dfdx_->apply(Thyra::NOTRANS, *dxdp, dfdp.ptr(), Scalar(1.0), Scalar(1.0));
491  }
492  if (dxdotdp != Teuchos::null && dfdp != Teuchos::null) {
493  if (my_dfdxdot_ == Teuchos::null)
494  my_dfdxdot_ = sens_residual_model_->create_W_op();
495  MEB::OutArgs<Scalar> meo = sens_residual_model_->createOutArgs();
496  meo.set_W_op(my_dfdxdot_);
497  if (me_inArgs.supports(MEB::IN_ARG_alpha))
498  me_inArgs.set_alpha(1.0);
499  if (me_inArgs.supports(MEB::IN_ARG_beta))
500  me_inArgs.set_beta(0.0);
501  if (me_inArgs.supports(MEB::IN_ARG_W_x_dot_dot_coeff))
502  me_inArgs.set_W_x_dot_dot_coeff(0.0);
503  sens_residual_model_->evalModel(me_inArgs, meo);
504  my_dfdxdot_->apply(Thyra::NOTRANS, *dxdotdp, dfdp.ptr(), Scalar(1.0), Scalar(1.0));
505  }
506  if (dxdotdotdp != Teuchos::null && dfdp != Teuchos::null) {
507  if (my_dfdxdotdot_ == Teuchos::null)
508  my_dfdxdotdot_ = sens_residual_model_->create_W_op();
509  MEB::OutArgs<Scalar> meo = sens_residual_model_->createOutArgs();
510  meo.set_W_op(my_dfdxdotdot_);
511  if (me_inArgs.supports(MEB::IN_ARG_alpha))
512  me_inArgs.set_alpha(0.0);
513  if (me_inArgs.supports(MEB::IN_ARG_beta))
514  me_inArgs.set_beta(0.0);
515  if (me_inArgs.supports(MEB::IN_ARG_W_x_dot_dot_coeff))
516  me_inArgs.set_W_x_dot_dot_coeff(1.0);
517  sens_residual_model_->evalModel(me_inArgs, meo);
518  my_dfdxdotdot_->apply(Thyra::NOTRANS, *dxdotdotdp, dfdp.ptr(), Scalar(1.0), Scalar(1.0));
519  }
520  }
521 
522  // Reduced dg/dp = dg/dx*dx/dp + dg/dp
523  if (g_index_ >= 0 && outArgs.get_g(0) != Teuchos::null) {
524  MEB::OutArgs<Scalar> meo = model_->createOutArgs();
525  RCP<Thyra::MultiVectorBase<Scalar> > dgdp =
526  rcp_dynamic_cast<DMVPV>(outArgs.get_g(0),true)->getNonconstMultiVector();
527  RCP<Thyra::MultiVectorBase<Scalar> > dgdp_trans;
528  MEB::DerivativeSupport dgdp_support =
529  meo.supports(MEB::OUT_ARG_DgDp, g_index_, p_index_);
530  if (dgdp_support.supports(MEB::DERIV_MV_JACOBIAN_FORM))
531  meo.set_DgDp(g_index_, p_index_,
532  MEB::Derivative<Scalar>(dgdp, MEB::DERIV_MV_JACOBIAN_FORM));
533  else if (dgdp_support.supports(MEB::DERIV_MV_GRADIENT_FORM)) {
534  dgdp_trans = createMembers(model_->get_p_space(p_index_),
535  num_response_);
536  meo.set_DgDp(g_index_, p_index_,
537  MEB::Derivative<Scalar>(dgdp_trans, MEB::DERIV_MV_GRADIENT_FORM));
538  }
539  else
541  true, std::logic_error,
542  "Operator form of dg/dp not supported for reduced sensitivity");
543 
544  if (!use_dgdp_as_tangent_ || dxdp == Teuchos::null) {
545  MEB::DerivativeSupport dgdx_support =
546  meo.supports(MEB::OUT_ARG_DgDx, g_index_);
547  if (dgdx_support.supports(MEB::DERIV_LINEAR_OP)) {
548  if (my_dgdx_ == Teuchos::null)
549  my_dgdx_ = model_->create_DgDx_op(g_index_);
550  meo.set_DgDx(g_index_, my_dgdx_);
551  }
552  else if (dgdx_support.supports(MEB::DERIV_MV_GRADIENT_FORM)) {
553  if (my_dgdx_mv_ == Teuchos::null)
554  my_dgdx_mv_ = createMembers(model_->get_g_space(g_index_), num_param_);
555  meo.set_DgDx(g_index_,
556  MEB::Derivative<Scalar>(my_dgdx_mv_,
557  MEB::DERIV_MV_GRADIENT_FORM));
558  }
559  else
561  true, std::logic_error,
562  "Jacobian form of dg/dx not supported for reduced sensitivity");
563  }
564 
565  // Clear dx/dp, dxdot/dp, dxdotdot/dp from inArgs if set and we are not
566  // using dg/dp as the tangent
567  if (!use_dgdp_as_tangent_ && use_dfdp_as_tangent_ &&
568  dxdp_vec != Teuchos::null)
569  me_inArgs.set_p(x_tangent_index_, Teuchos::null);
570  if (!use_dgdp_as_tangent_ && use_dfdp_as_tangent_ &&
571  dxdotdp_vec != Teuchos::null)
572  me_inArgs.set_p(xdot_tangent_index_, Teuchos::null);
573  if (!use_dgdp_as_tangent_ && use_dfdp_as_tangent_ &&
574  dxdotdotdp_vec != Teuchos::null)
575  me_inArgs.set_p(xdotdot_tangent_index_, Teuchos::null);
576 
577  // Set dx/dp, dxdot/dp, dxdodot/dp if necessary
578  if (use_dgdp_as_tangent_ && !use_dfdp_as_tangent_ &&
579  dxdp_vec != Teuchos::null)
580  me_inArgs.set_p(x_tangent_index_, dxdp_vec);
581  if (use_dgdp_as_tangent_ && !use_dfdp_as_tangent_ &&
582  dxdotdp_vec != Teuchos::null)
583  me_inArgs.set_p(xdot_tangent_index_, dxdotdp_vec);
584  if (use_dgdp_as_tangent_ && !use_dfdp_as_tangent_ &&
585  dxdotdotdp_vec != Teuchos::null)
586  me_inArgs.set_p(xdotdot_tangent_index_, dxdotdotdp_vec);
587 
588  model_->evalModel(me_inArgs, meo);
589 
590  // transpose reduced dg/dp if necessary
591  if (dgdp_trans != Teuchos::null) {
593  Thyra::DetachedMultiVectorView<Scalar> dgdp_trans_view(*dgdp_trans);
594  for (int i=0; i<num_param_; ++i)
595  for (int j=0; j<num_response_; ++j)
596  dgdp_view(j,i) = dgdp_trans_view(i,j);
597  }
598 
599  // Compute (dg/dx) * (dx/dp) + (dg/dp) if the underlying ME doesn't already
600  // do this.
601  if (!use_dgdp_as_tangent_ && dxdp != Teuchos::null) {
602  MEB::DerivativeSupport dgdx_support =
603  meo.supports(MEB::OUT_ARG_DgDx, g_index_);
604  if (dgdx_support.supports(MEB::DERIV_LINEAR_OP))
605  my_dgdx_->apply(Thyra::NOTRANS, *dxdp, dgdp.ptr(), Scalar(1.0), Scalar(1.0));
606  else if (dgdx_support.supports(MEB::DERIV_MV_GRADIENT_FORM))
607  my_dgdx_mv_->apply(Thyra::TRANS, *dxdp, dgdp.ptr(), Scalar(1.0), Scalar(1.0));
608  else
610  true, std::logic_error,
611  "Jacobian form of dg/dx not supported for reduced sensitivity");
612  }
613  }
614 }
615 
616 template<class Scalar>
620 {
621  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
622  pl->set<bool>("Use DfDp as Tangent", false);
623  pl->set<bool>("Use DgDp as Tangent", false);
624  pl->set<int>("Sensitivity Parameter Index", 0);
625  pl->set<int>("Response Function Index", -1);
626  pl->set<int>("Sensitivity X Tangent Index", 1);
627  pl->set<int>("Sensitivity X-Dot Tangent Index", 2);
628  pl->set<int>("Sensitivity X-Dot-Dot Tangent Index", 3);
629  return pl;
630 }
631 
632 } // namespace Tempus
633 
634 #endif
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_p_space(int p) const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_g_space(int j) const
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_DgDx_dot_op(int j) const
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)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int p) const
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_DgDp_op(int j, int l) const
Evaluation< VectorBase< Scalar > > get_g(int j) const
Evaluation< VectorBase< Scalar > > get_f() const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_DgDx_op(int j) const
bool supports(EOutArgsMembers arg) const
Derivative< Scalar > get_DgDp(int j, int l) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
CombinedForwardSensitivityModelEvaluator(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 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.
void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
virtual std::string description() const
void validateParametersAndSetDefaults(ParameterList const &validParamList, int const depth=1000)
Teuchos::RCP< const Thyra::LinearOpWithSolveFactoryBase< Scalar > > get_W_factory() const
Implicit concrete LinearOpBase subclass that takes a flattended out multi-vector and performs a multi...
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_f_space() const
size_type size() const
Derivative< Scalar > get_DgDx_dot(int j) const
RCP< const VectorBase< Scalar > > get_x_dot_dot() const
#define TEUCHOS_ASSERT(assertion_test)
RCP< LinearOpBase< Scalar > > get_W_op() const
RCP< const VectorBase< Scalar > > get_x() const
RCP< const VectorBase< Scalar > > get_p(int l) const
Teuchos::Range1D Range1D
Derivative< Scalar > get_DgDx(int j) const