Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tempus_IntegratorAdjointSensitivity_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_IntegratorAdjointSensitivity_impl_hpp
10 #define Tempus_IntegratorAdjointSensitivity_impl_hpp
11 
13 #include "Thyra_DefaultMultiVectorProductVectorSpace.hpp"
14 #include "Thyra_DefaultMultiVectorProductVector.hpp"
15 #include "Thyra_DefaultProductVectorSpace.hpp"
16 #include "Thyra_DefaultProductVector.hpp"
17 #include "Thyra_VectorStdOps.hpp"
18 #include "Thyra_MultiVectorStdOps.hpp"
20 #include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
21 
22 namespace Tempus {
23 
24 template <class Scalar>
27  const Teuchos::RCP<IntegratorBasic<Scalar> > &state_integrator,
28  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> > &adjoint_model,
30  const Teuchos::RCP<IntegratorBasic<Scalar> > &adjoint_integrator,
31  const Teuchos::RCP<SolutionHistory<Scalar> > &solutionHistory,
32  const int p_index,
33  const int g_index,
34  const bool g_depends_on_p,
35  const bool f_depends_on_p,
36  const bool ic_depends_on_p,
37  const bool mass_matrix_is_identity)
38  : model_(model)
39  , state_integrator_(state_integrator)
40  , adjoint_model_(adjoint_model)
41  , adjoint_aux_model_(adjoint_aux_model)
42  , adjoint_integrator_(adjoint_integrator)
43  , solutionHistory_(solutionHistory)
44  , p_index_(p_index)
45  , g_index_(g_index)
46  , g_depends_on_p_(g_depends_on_p)
47  , f_depends_on_p_(f_depends_on_p)
48  , ic_depends_on_p_(ic_depends_on_p)
49  , mass_matrix_is_identity_(mass_matrix_is_identity)
50  , stepMode_(SensitivityStepMode::Forward)
51 {
52 
54  getStepper()->getUseFSAL(), std::logic_error,
55  "Error - IntegratorAdjointSensitivity(): Cannot use FSAL with\n"
56  " IntegratorAdjointSensitivity, because the state and adjoint\n"
57  " integrators require ModelEvaluator evaluation in the\n"
58  " constructor to make the initial conditions consistent.\n"
59  " For the adjoint integrator, this requires special construction\n"
60  " which has not been implemented yet.\n");
61 }
62 
63 template<class Scalar>
66 {
67  state_integrator_ = createIntegratorBasic<Scalar>();
68  adjoint_integrator_ = createIntegratorBasic<Scalar>();
69  stepMode_ = SensitivityStepMode::Forward;
70 }
71 
72 template<class Scalar>
73 bool
76 {
77  const Scalar tfinal =
78  state_integrator_->getTimeStepControl()->getFinalTime();
79  return advanceTime(tfinal);
80 }
81 
82 template<class Scalar>
83 bool
85 advanceTime(const Scalar timeFinal)
86 {
87  TEMPUS_FUNC_TIME_MONITOR_DIFF("Tempus::IntegratorAdjointSensitivity::advanceTime()", TEMPUS_AS_AT);
88 
89  using Teuchos::RCP;
90  using Teuchos::rcp_dynamic_cast;
91  using Thyra::VectorBase;
94  using Thyra::LinearOpBase;
99  using Thyra::createMember;
100  using Thyra::createMembers;
101  using Thyra::assign;
102  typedef Thyra::ModelEvaluatorBase MEB;
105 
106  // Get initial state for later
107  RCP<const SolutionHistory<Scalar> > state_solution_history =
108  state_integrator_->getSolutionHistory();
109  RCP<const SolutionState<Scalar> > initial_state =
110  (*state_solution_history)[0];
111 
112  // Run state integrator and get solution
113  bool state_status = true;
114  {
115  TEMPUS_FUNC_TIME_MONITOR_DIFF("Tempus::IntegratorAdjointSensitivity::advanceTime::state", TEMPUS_AS_AT_FWD);
116  stepMode_ = SensitivityStepMode::Forward;
117  state_status = state_integrator_->advanceTime(timeFinal);
118  }
119 
120  // For at least some time-stepping methods, the time of the last time step
121  // may not be timeFinal (e.g., it may be greater by at most delta_t).
122  // But since the adjoint model requires timeFinal in its formulation, reset
123  // it to the achieved final time.
124  adjoint_aux_model_->setFinalTime(state_integrator_->getTime());
125 
126  // Set solution history in adjoint stepper
127  adjoint_aux_model_->setForwardSolutionHistory(state_solution_history);
128 
129  // Compute dg/dx
130  RCP<const VectorSpaceBase<Scalar> > g_space = model_->get_g_space(g_index_);
131  RCP<const VectorSpaceBase<Scalar> > x_space = model_->get_x_space();
132  const int num_g = g_space->dim();
133  RCP<MultiVectorBase<Scalar> > dgdx = createMembers(x_space, num_g);
134  MEB::InArgs<Scalar> inargs = model_->getNominalValues();
135  RCP<const SolutionState<Scalar> > state =
136  state_solution_history->getCurrentState();
137  inargs.set_t(state->getTime());
138  inargs.set_x(state->getX());
139  inargs.set_x_dot(state->getXDot());
140  MEB::OutArgs<Scalar> outargs = model_->createOutArgs();
141  MEB::OutArgs<Scalar> adj_outargs = adjoint_model_->createOutArgs();
142  outargs.set_DgDx(g_index_,
143  MEB::Derivative<Scalar>(dgdx, MEB::DERIV_MV_GRADIENT_FORM));
144  model_->evalModel(inargs, outargs);
145  outargs.set_DgDx(g_index_, MEB::Derivative<Scalar>());
146 
147  // Compute ICs == [ (df/dx_dot)^{-T} (dg/dx)^T; 0 ]
148  // For explicit form, we are relying on the user to inform us the
149  // the mass matrix is the identity. It would be nice to be able to determine
150  // somehow automatically that we are using an explicit stepper.
151  RCP<DPV> adjoint_init =
152  rcp_dynamic_cast<DPV>(Thyra::createMember(adjoint_aux_model_->get_x_space()));
153  RCP<MultiVectorBase<Scalar> > adjoint_init_mv =
154  rcp_dynamic_cast<DMVPV>(adjoint_init->getNonconstVectorBlock(0))->getNonconstMultiVector();
155  assign(adjoint_init->getNonconstVectorBlock(1).ptr(),
157  if (mass_matrix_is_identity_)
158  assign(adjoint_init_mv.ptr(), *dgdx);
159  else {
160  inargs.set_alpha(1.0);
161  inargs.set_beta(0.0);
162  RCP<LinearOpWithSolveBase<Scalar> > W;
163  if (adj_outargs.supports(MEB::OUT_ARG_W)) {
164  // Model supports W
165  W = adjoint_model_->create_W();
166  adj_outargs.set_W(W);
167  adjoint_model_->evalModel(inargs, adj_outargs);
168  adj_outargs.set_W(Teuchos::null);
169  }
170  else {
171  // Otherwise model must support a W_op and W factory
172  RCP<const LinearOpWithSolveFactoryBase<Scalar> > lowsfb =
173  adjoint_model_->get_W_factory();
175  lowsfb == Teuchos::null, std::logic_error,
176  "Adjoint ME must support W out-arg or provide a W_factory for non-identity mass matrix");
177 
178  // Compute W_op (and W_prec if supported)
179  RCP<LinearOpBase<Scalar> > W_op = adjoint_model_->create_W_op();
180  adj_outargs.set_W_op(W_op);
181  RCP<PreconditionerFactoryBase<Scalar> > prec_factory =
182  lowsfb->getPreconditionerFactory();
183  RCP<PreconditionerBase<Scalar> > W_prec;
184  if (prec_factory != Teuchos::null)
185  W_prec = prec_factory->createPrec();
186  else if (adj_outargs.supports(MEB::OUT_ARG_W_prec)) {
187  W_prec = adjoint_model_->create_W_prec();
188  adj_outargs.set_W_prec(W_prec);
189  }
190  adjoint_model_->evalModel(inargs, adj_outargs);
191  adj_outargs.set_W_op(Teuchos::null);
192  if (adj_outargs.supports(MEB::OUT_ARG_W_prec))
193  adj_outargs.set_W_prec(Teuchos::null);
194 
195  // Create and initialize W
196  W = lowsfb->createOp();
197  if (W_prec != Teuchos::null) {
198  if (prec_factory != Teuchos::null)
199  prec_factory->initializePrec(
200  Thyra::defaultLinearOpSource<Scalar>(W_op), W_prec.get());
201  Thyra::initializePreconditionedOp<Scalar>(
202  *lowsfb, W_op, W_prec, W.ptr());
203  }
204  else
205  Thyra::initializeOp<Scalar>(*lowsfb, W_op, W.ptr());
206  }
208  W == Teuchos::null, std::logic_error,
209  "A null W has been encountered in Tempus::IntegratorAdjointSensitivity::advanceTime!\n");
210  // Initialize adjoint_init_mv to zero before solve for linear solvers that
211  // use what is passed in as the initial guess
212  assign(adjoint_init_mv.ptr(), Teuchos::ScalarTraits<Scalar>::zero());
213  W->solve(Thyra::NOTRANS, *dgdx, adjoint_init_mv.ptr());
214  }
215 
216  // Run sensitivity integrator and get solution
217  bool sens_status = true;
218  {
219  TEMPUS_FUNC_TIME_MONITOR_DIFF("Tempus::IntegratorAdjointSensitivity::advanceTime::adjoint", TEMPUS_AS_AT_ADJ);
220  stepMode_ = SensitivityStepMode::Adjoint;
221  const Scalar tinit = adjoint_integrator_->getTimeStepControl()->getInitTime();
222  adjoint_integrator_->initializeSolutionHistory(tinit, adjoint_init);
223  sens_status = adjoint_integrator_->advanceTime(timeFinal);
224  }
225  RCP<const SolutionHistory<Scalar> > adjoint_solution_history =
226  adjoint_integrator_->getSolutionHistory();
227 
228  // Compute dg/dp at final time T
229  RCP<const VectorSpaceBase<Scalar> > p_space = model_->get_p_space(p_index_);
230  dgdp_ = createMembers(p_space, num_g);
231  if (g_depends_on_p_) {
232  MEB::DerivativeSupport dgdp_support =
233  outargs.supports(MEB::OUT_ARG_DgDp, g_index_, p_index_);
234  if (dgdp_support.supports(MEB::DERIV_MV_GRADIENT_FORM)) {
235  outargs.set_DgDp(g_index_, p_index_,
236  MEB::Derivative<Scalar>(dgdp_,
237  MEB::DERIV_MV_GRADIENT_FORM));
238  model_->evalModel(inargs, outargs);
239  }
240  else if (dgdp_support.supports(MEB::DERIV_MV_JACOBIAN_FORM)) {
241  const int num_p = p_space->dim();
242  RCP<MultiVectorBase<Scalar> > dgdp_trans =
243  createMembers(g_space, num_p);
244  outargs.set_DgDp(g_index_, p_index_,
245  MEB::Derivative<Scalar>(dgdp_trans,
246  MEB::DERIV_MV_JACOBIAN_FORM));
247  model_->evalModel(inargs, outargs);
248  Thyra::DetachedMultiVectorView<Scalar> dgdp_view(*dgdp_);
249  Thyra::DetachedMultiVectorView<Scalar> dgdp_trans_view(*dgdp_trans);
250  for (int i=0; i<num_p; ++i)
251  for (int j=0; j<num_g; ++j)
252  dgdp_view(i,j) = dgdp_trans_view(j,i);
253  }
254  else
255  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
256  "Invalid dg/dp support");
257  outargs.set_DgDp(g_index_, p_index_, MEB::Derivative<Scalar>());
258  }
259  else
260  assign(dgdp_.ptr(), Scalar(0.0));
261 
262  // Add in initial condition term = (dx/dp^T(0))*(df/dx_dot^T(0))*y(0)
263  // If dxdp_init_ is null, assume it is zero
264  if (ic_depends_on_p_ && dxdp_init_ != Teuchos::null) {
265  RCP<const SolutionState<Scalar> > adjoint_state =
266  adjoint_solution_history->getCurrentState();
267  RCP<const VectorBase<Scalar> > adjoint_x =
268  rcp_dynamic_cast<const DPV>(adjoint_state->getX())->getVectorBlock(0);
269  RCP<const MultiVectorBase<Scalar> > adjoint_mv =
270  rcp_dynamic_cast<const DMVPV>(adjoint_x)->getMultiVector();
271  if (mass_matrix_is_identity_)
272  dxdp_init_->apply(Thyra::CONJTRANS, *adjoint_mv, dgdp_.ptr(), Scalar(1.0),
273  Scalar(1.0));
274  else {
275  inargs.set_t(initial_state->getTime());
276  inargs.set_x(initial_state->getX());
277  inargs.set_x_dot(initial_state->getXDot());
278  inargs.set_alpha(1.0);
279  inargs.set_beta(0.0);
280  RCP<LinearOpBase<Scalar> > W_op = adjoint_model_->create_W_op();
281  adj_outargs.set_W_op(W_op);
282  adjoint_model_->evalModel(inargs, adj_outargs);
283  adj_outargs.set_W_op(Teuchos::null);
284  RCP<MultiVectorBase<Scalar> > tmp = createMembers(x_space, num_g);
285  W_op->apply(Thyra::NOTRANS, *adjoint_mv, tmp.ptr(), Scalar(1.0),
286  Scalar(0.0));
287  dxdp_init_->apply(Thyra::CONJTRANS, *tmp, dgdp_.ptr(), Scalar(1.0),
288  Scalar(1.0));
289  }
290  }
291 
292  // Add in model parameter term = \int_0^T( (df/dp^T(t)*y(t) )dt which
293  // is computed during the adjoint integration as an auxiliary integral
294  // (2nd block of the solution vector)
295  if (f_depends_on_p_) {
296  RCP<const SolutionState<Scalar> > adjoint_state =
297  adjoint_solution_history->getCurrentState();
298  RCP<const VectorBase<Scalar> > z =
299  rcp_dynamic_cast<const DPV>(adjoint_state->getX())->getVectorBlock(1);
300  RCP<const MultiVectorBase<Scalar> > z_mv =
301  rcp_dynamic_cast<const DMVPV>(z)->getMultiVector();
302  Thyra::V_VmV(dgdp_.ptr(), *dgdp_, *z_mv);
303  }
304 
305  buildSolutionHistory(state_solution_history, adjoint_solution_history);
306 
307  return state_status && sens_status;
308 }
309 
310 template<class Scalar>
311 Scalar
313 getTime() const
314 {
315  return solutionHistory_->getCurrentTime();
316 }
317 
318 template<class Scalar>
319 int
321 getIndex() const
322 {
323  return solutionHistory_->getCurrentIndex();
324 }
325 
326 template<class Scalar>
327 Status
329 getStatus() const
330 {
331  Status state_status = state_integrator_->getStatus();
332  Status sens_status = adjoint_integrator_->getStatus();
333  if (state_status == FAILED || sens_status == FAILED)
334  return FAILED;
335  if (state_status == WORKING || sens_status == WORKING)
336  return WORKING;
337  return PASSED;
338 }
339 
340 template <class Scalar>
342  state_integrator_->setStatus(st);
343  adjoint_integrator_->setStatus(st);
344 }
345 
346 template<class Scalar>
349 getStepper() const
350 {
351  return state_integrator_->getStepper();
352 }
353 
354 template<class Scalar>
358 {
359  return solutionHistory_;
360 }
361 
362 template<class Scalar>
366 {
367  return state_integrator_->getSolutionHistory();
368 }
369 
370 template<class Scalar>
374 {
375  return adjoint_integrator_->getSolutionHistory();
376 }
377 
378 template<class Scalar>
382 {
383  return solutionHistory_;
384 }
385 
386 template<class Scalar>
390 {
391  return state_integrator_->getTimeStepControl();
392 }
393 
394 template<class Scalar>
398 {
399  return state_integrator_->getNonConstTimeStepControl();
400 }
401 
402 template<class Scalar>
406 {
407  return state_integrator_->getNonConstTimeStepControl();
408 }
409 
410 template<class Scalar>
414 {
415  return adjoint_integrator_->getNonConstTimeStepControl();
416 }
417 
418 template<class Scalar>
423  Teuchos::RCP<const Thyra::VectorBase<Scalar> > xdotdot0,
425  Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > /* DxdotDp0 */,
426  Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > /* DxdotdotDp0 */)
427 {
428  state_integrator_->initializeSolutionHistory(t0, x0, xdot0, xdotdot0);
429  dxdp_init_ = DxDp0;
430 }
431 
432 template<class Scalar>
436 {
437  return state_integrator_->getObserver();
438 }
439 
440 template<class Scalar>
443 {
444  state_integrator_->setObserver(obs);
445  // ETP 1/12/22 Disabling passing of the observer to the adjoint
446  // integrator to work around issues in Piro
447  //adjoint_integrator_->setObserver(obs);
448 }
449 
450 template<class Scalar>
453 {
454  state_integrator_->initialize();
455  adjoint_integrator_->initialize();
456 }
457 
458 template<class Scalar>
461 getX() const
462 {
463  return state_integrator_->getX();
464 }
465 
466 template<class Scalar>
469 getXDot() const
470 {
471  return state_integrator_->getXDot();
472 }
473 
474 template<class Scalar>
477 getXDotDot() const
478 {
479  return state_integrator_->getXDotDot();
480 }
481 
482 template<class Scalar>
485 getY() const
486 {
487  using Teuchos::RCP;
488  using Teuchos::rcp_dynamic_cast;
491  RCP<const DPV> pv =
492  rcp_dynamic_cast<const DPV>(adjoint_integrator_->getX());
493  RCP<const DMVPV> mvpv =
494  rcp_dynamic_cast<const DMVPV>(pv->getVectorBlock(0));
495  return mvpv->getMultiVector();
496 }
497 
498 template<class Scalar>
501 getYDot() const
502 {
503  using Teuchos::RCP;
504  using Teuchos::rcp_dynamic_cast;
507  RCP<const DPV> pv =
508  rcp_dynamic_cast<const DPV>(adjoint_integrator_->getXDot());
509  RCP<const DMVPV> mvpv =
510  rcp_dynamic_cast<const DMVPV>(pv->getVectorBlock(0));
511  return mvpv->getMultiVector();
512 }
513 
514 template<class Scalar>
517 getYDotDot() const
518 {
519  using Teuchos::RCP;
520  using Teuchos::rcp_dynamic_cast;
523  RCP<const DPV> pv =
524  rcp_dynamic_cast<const DPV>(adjoint_integrator_->getXDotDot());
525  RCP<const DMVPV> mvpv =
526  rcp_dynamic_cast<const DMVPV>(pv->getVectorBlock(0));
527  return mvpv->getMultiVector();
528 }
529 
530 template<class Scalar>
533 getDgDp() const
534 {
535  return dgdp_;
536 }
537 
538 template<class Scalar>
539 std::string
541 description() const
542 {
543  std::string name = "Tempus::IntegratorAdjointSensitivity";
544  return(name);
545 }
546 
547 template<class Scalar>
548 void
552  const Teuchos::EVerbosityLevel verbLevel) const
553 {
554  auto l_out = Teuchos::fancyOStream( out.getOStream() );
555  Teuchos::OSTab ostab(*l_out, 2, this->description());
556  l_out->setOutputToRootOnly(0);
557 
558  *l_out << description() << "::describe" << std::endl;
559  state_integrator_->describe(*l_out, verbLevel);
560  adjoint_integrator_->describe(*l_out, verbLevel);
561 }
562 
563 template<class Scalar>
566 getStepMode() const
567 {
568  return stepMode_;
569 }
570 
571 template <class Scalar>
576  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& adjoint_model,
578 {
580  using Teuchos::RCP;
581  using Teuchos::rcp;
582 
583  RCP<ParameterList> spl = Teuchos::parameterList();
584  if (inputPL != Teuchos::null)
585  {
586  *spl = inputPL->sublist("Sensitivities");
587  }
588  if (spl->isParameter("Response Depends on Parameters"))
589  spl->remove("Response Depends on Parameters");
590  if (spl->isParameter("Residual Depends on Parameters"))
591  spl->remove("Residual Depends on Parameters");
592  if (spl->isParameter("IC Depends on Parameters"))
593  spl->remove("IC Depends on Parameters");
594 
595  const Scalar tinit = state_integrator_->getTimeStepControl()->getInitTime();
596  const Scalar tfinal = state_integrator_->getTimeStepControl()->getFinalTime();
598  model, adjoint_model, tinit, tfinal, spl));
599 }
600 
601 template<class Scalar>
602 void
605  const Teuchos::RCP<const SolutionHistory<Scalar> >& state_solution_history,
606  const Teuchos::RCP<const SolutionHistory<Scalar> >& adjoint_solution_history)
607 {
608  using Teuchos::RCP;
609  using Teuchos::rcp;
610  using Teuchos::rcp_dynamic_cast;
612  using Thyra::VectorBase;
615  using Thyra::createMembers;
616  using Thyra::multiVectorProductVector;
617  using Thyra::assign;
620 
621  RCP<const VectorSpaceBase<Scalar> > x_space = model_->get_x_space();
622  RCP<const VectorSpaceBase<Scalar> > adjoint_space =
623  rcp_dynamic_cast<const DPVS>(adjoint_aux_model_->get_x_space())->getBlock(0);
625  spaces[0] = x_space;
626  spaces[1] = adjoint_space;
627  RCP<const DPVS > prod_space = Thyra::productVectorSpace(spaces());
628 
629  int num_states = state_solution_history->getNumStates();
630  const Scalar t_init = state_integrator_->getTimeStepControl()->getInitTime();
631  const Scalar t_final = state_integrator_->getTime();
632  for (int i=0; i<num_states; ++i) {
633  RCP<const SolutionState<Scalar> > forward_state =
634  (*state_solution_history)[i];
635  RCP<const SolutionState<Scalar> > adjoint_state =
636  adjoint_solution_history->findState(t_final+t_init-forward_state->getTime());
637 
638  // X
639  RCP<DPV> x = Thyra::defaultProductVector(prod_space);
640  RCP<const VectorBase<Scalar> > adjoint_x =
641  rcp_dynamic_cast<const DPV>(adjoint_state->getX())->getVectorBlock(0);
642  assign(x->getNonconstVectorBlock(0).ptr(), *(forward_state->getX()));
643  assign(x->getNonconstVectorBlock(1).ptr(), *(adjoint_x));
644  RCP<VectorBase<Scalar> > x_b = x;
645 
646  // X-Dot
647  RCP<DPV> x_dot = Thyra::defaultProductVector(prod_space);
648  RCP<const VectorBase<Scalar> > adjoint_x_dot =
649  rcp_dynamic_cast<const DPV>(adjoint_state->getXDot())->getVectorBlock(0);
650  assign(x_dot->getNonconstVectorBlock(0).ptr(), *(forward_state->getXDot()));
651  assign(x_dot->getNonconstVectorBlock(1).ptr(), *(adjoint_x_dot));
652  RCP<VectorBase<Scalar> > x_dot_b = x_dot;
653 
654  // X-Dot-Dot
655  RCP<DPV> x_dot_dot;
656  if (forward_state->getXDotDot() != Teuchos::null) {
657  x_dot_dot = Thyra::defaultProductVector(prod_space);
658  RCP<const VectorBase<Scalar> > adjoint_x_dot_dot =
659  rcp_dynamic_cast<const DPV>(
660  adjoint_state->getXDotDot())->getVectorBlock(0);
661  assign(x_dot_dot->getNonconstVectorBlock(0).ptr(),
662  *(forward_state->getXDotDot()));
663  assign(x_dot_dot->getNonconstVectorBlock(1).ptr(),
664  *(adjoint_x_dot_dot));
665  }
666  RCP<VectorBase<Scalar> > x_dot_dot_b = x_dot_dot;
667 
668  RCP<SolutionState<Scalar> > prod_state = forward_state->clone();
669  prod_state->setX(x_b);
670  prod_state->setXDot(x_dot_b);
671  prod_state->setXDotDot(x_dot_dot_b);
672  prod_state->setPhysicsState(Teuchos::null);
673  solutionHistory_->addState(prod_state);
674  }
675 }
676 
678 template<class Scalar>
683  ,const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& adjoint_model)
684 {
685  // set the parameters
686  Teuchos::RCP<Teuchos::ParameterList> spl = Teuchos::parameterList();
687  if (inputPL != Teuchos::null)
688  *spl = inputPL->sublist("Sensitivities");
689 
690  int p_index = spl->get<int>("Sensitivity Parameter Index", 0);
691  int g_index = spl->get<int>("Response Function Index", 0);
692  bool g_depends_on_p = spl->get<bool>("Response Depends on Parameters", true);
693  bool f_depends_on_p = spl->get<bool>("Residual Depends on Parameters", true);
694  bool ic_depends_on_p = spl->get<bool>("IC Depends on Parameters", true);
695  bool mass_matrix_is_identity = spl->get<bool>("Mass Matrix Is Identity", false);
696 
697  auto state_integrator = createIntegratorBasic<Scalar>(inputPL, model);
698 
699  // createAdjointModel
700  if (spl->isParameter("Response Depends on Parameters"))
701  spl->remove("Response Depends on Parameters");
702  if (spl->isParameter("Residual Depends on Parameters"))
703  spl->remove("Residual Depends on Parameters");
704  if (spl->isParameter("IC Depends on Parameters"))
705  spl->remove("IC Depends on Parameters");
706 
707  const Scalar tinit = state_integrator->getTimeStepControl()->getInitTime();
708  const Scalar tfinal = state_integrator->getTimeStepControl()->getFinalTime();
709  //auto adjoint_model = Teuchos::rcp(new AdjointAuxSensitivityModelEvaluator<Scalar>(model, tfinal, spl));
710  //TODO: where is the adjoint ME coming from?
711 
712  Teuchos::RCP<Thyra::ModelEvaluator<Scalar>> adjt_model = adjoint_model;
713  if (adjoint_model == Teuchos::null)
714  adjt_model = Thyra::implicitAdjointModelEvaluator(model);
715 
716  auto adjoint_aux_model = Teuchos::rcp(new AdjointAuxSensitivityModelEvaluator<Scalar>(model, adjt_model, tinit, tfinal, spl));
717 
718  // Create combined solution histories combining the forward and adjoint
719  // solutions. We do not include the auxiliary part from the adjoint solution.
720  auto integrator_name = inputPL->get<std::string>("Integrator Name");
721  auto integratorPL = Teuchos::sublist(inputPL, integrator_name, true);
722  auto shPL = Teuchos::sublist(integratorPL, "Solution History", true);
723  auto combined_solution_History = createSolutionHistoryPL<Scalar>(shPL);
724 
725  auto adjoint_integrator = createIntegratorBasic<Scalar>(inputPL, adjoint_aux_model);
726 
728  model, state_integrator, adjt_model, adjoint_aux_model, adjoint_integrator, combined_solution_History, p_index, g_index, g_depends_on_p,
729  f_depends_on_p, ic_depends_on_p, mass_matrix_is_identity));
730 
731  return (integrator);
732 }
733 
735 template<class Scalar>
738 {
741  return(integrator);
742 }
743 
744 } // namespace Tempus
745 #endif // Tempus_IntegratorAdjointSensitivity_impl_hpp
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const override
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getYDot() const
Get the current time derivative of the adjoint solution, ydot.
SensitivityStepMode getStepMode() const
What mode the current time integration step is in.
virtual Teuchos::RCP< TimeStepControl< Scalar > > getNonConstTimeStepControl() override
T & get(const std::string &name, T def_value)
Teuchos::RCP< IntegratorAdjointSensitivity< Scalar > > createIntegratorAdjointSensitivity(Teuchos::RCP< Teuchos::ParameterList > pList, const Teuchos::RCP< Thyra::ModelEvaluator< Scalar >> &model, const Teuchos::RCP< Thyra::ModelEvaluator< Scalar > > &adjoint_model=Teuchos::null)
Nonmember constructor.
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getXDotDot() const
Get the current second time derivative of the solution, xdotdot.
virtual void setStatus(const Status st) override
Set Status.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
virtual Teuchos::RCP< SolutionHistory< Scalar > > getNonConstSolutionHistory() override
Get the SolutionHistory.
ModelEvaluator for forming adjoint sensitivity equations.
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getXDot() const
Get the current time derivative of the solution, xdot.
void buildSolutionHistory(const Teuchos::RCP< const SolutionHistory< Scalar > > &state_solution_history, const Teuchos::RCP< const SolutionHistory< Scalar > > &adjoint_solution_history)
virtual Teuchos::RCP< IntegratorObserver< Scalar > > getObserver()
Get the Observer.
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getX() const
Get the current solution, x.
bool isParameter(const std::string &name) const
bool remove(std::string const &name, bool throwIfNotExists=true)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual Teuchos::RCP< const SolutionHistory< Scalar > > getSolutionHistory() const override
Get the SolutionHistory.
Status
Status for the Integrator, the Stepper and the SolutionState.
IntegratorObserver class for time integrators.
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getYDotDot() const
Get the current second time derivative of the adjoint solution, ydotdot.
virtual void initializeSolutionHistory(Scalar t0, Teuchos::RCP< const Thyra::VectorBase< Scalar > > x0, Teuchos::RCP< const Thyra::VectorBase< Scalar > > xdot0=Teuchos::null, Teuchos::RCP< const Thyra::VectorBase< Scalar > > xdotdot0=Teuchos::null, Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > DxDp0=Teuchos::null, Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > DxdotDp0=Teuchos::null, Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > DxdotdotDp0=Teuchos::null)
Set the initial state from Thyra::VectorBase(s)
virtual Teuchos::RCP< const TimeStepControl< Scalar > > getTimeStepControl() const override
Get the TimeStepControl.
RCP< const MultiVectorBase< Scalar > > getMultiVector() const
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
virtual bool advanceTime()
Advance the solution to timeMax, and return true if successful.
Teuchos::RCP< TimeStepControl< Scalar > > getSensNonConstTimeStepControl()
Teuchos::RCP< const SolutionHistory< Scalar > > getSensSolutionHistory() const
virtual void setObserver(Teuchos::RCP< IntegratorObserver< Scalar > > obs=Teuchos::null)
Set the Observer.
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getY() const
Get the current adjoint solution, y.
Teuchos::RCP< const SolutionHistory< Scalar > > getStateSolutionHistory() const
RCP< ImplicitAdjointModelEvaluator< Scalar > > implicitAdjointModelEvaluator(const RCP< const ModelEvaluator< Scalar > > &model)
RCP< std::basic_ostream< char_type, traits_type > > getOStream()
Time integrator suitable for adjoint sensitivity analysis.
virtual int getIndex() const override
Get current index.
Teuchos::RCP< TimeStepControl< Scalar > > getStateNonConstTimeStepControl()
Teuchos::RCP< AdjointAuxSensitivityModelEvaluator< Scalar > > createAdjointModel(const Teuchos::RCP< Thyra::ModelEvaluator< Scalar > > &model, const Teuchos::RCP< Thyra::ModelEvaluator< Scalar > > &adjoint_model, const Teuchos::RCP< Teuchos::ParameterList > &inputPL)
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
virtual Scalar getTime() const override
Get current time.
virtual Status getStatus() const override
Get Status.
virtual void initialize()
Initializes the Integrator after set* function calls.
virtual Teuchos::RCP< Stepper< Scalar > > getStepper() const override
Get the Stepper.
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getDgDp() const
Return adjoint sensitivity stored in gradient format.
IntegratorAdjointSensitivity()
Constructor that requires a subsequent setParameterList, setStepper, and initialize calls...