9 #ifndef Tempus_IntegratorPseudoTransientAdjointSensitivity_impl_hpp
10 #define Tempus_IntegratorPseudoTransientAdjointSensitivity_impl_hpp
12 #include "Thyra_DefaultMultiVectorProductVector.hpp"
13 #include "Thyra_VectorStdOps.hpp"
14 #include "Thyra_MultiVectorStdOps.hpp"
20 template<
class Scalar>
29 adjoint_residual_model_ = adjoint_residual_model;
30 adjoint_solve_model_ = adjoint_solve_model;
31 state_integrator_ = createIntegratorBasic<Scalar>(inputPL, model_);
32 sens_model_ = createSensitivityModel(model_, adjoint_residual_model_,
33 adjoint_solve_model_, inputPL);
34 sens_integrator_ = createIntegratorBasic<Scalar>(inputPL, sens_model_);
36 do_forward_integration_ =
true;
37 do_adjoint_integration_ =
true;
40 template<
class Scalar>
46 std::string stepperType)
49 adjoint_residual_model_ = adjoint_residual_model;
50 adjoint_solve_model_ = adjoint_solve_model;
51 state_integrator_ = createIntegratorBasic<Scalar>(model_, stepperType);
52 sens_model_ = createSensitivityModel(model_, adjoint_residual_model_,
53 adjoint_solve_model_, Teuchos::null);
54 sens_integrator_ = createIntegratorBasic<Scalar>(sens_model_, stepperType);
56 do_forward_integration_ =
true;
57 do_adjoint_integration_ =
true;
60 template<
class Scalar>
67 inputPL, model, adjoint_model, adjoint_model)
71 template<
class Scalar>
76 std::string stepperType) :
78 model, adjoint_model, adjoint_model, stepperType)
82 template<
class Scalar>
92 template<
class Scalar>
96 std::string stepperType) :
102 template<
class Scalar>
106 state_integrator_ = createIntegratorBasic<Scalar>();
107 sens_integrator_ = createIntegratorBasic<Scalar>();
111 template<
class Scalar>
116 const Scalar tfinal =
117 state_integrator_->getTimeStepControl()->getFinalTime();
118 return advanceTime(tfinal);
121 template<
class Scalar>
126 TEMPUS_FUNC_TIME_MONITOR_DIFF(
"Tempus::IntegratorPseudoTransientAdjointSensitivity::advanceTime()", TEMPUS_PTAS_AT);
132 bool state_status =
true;
133 if (do_forward_integration_) {
134 TEMPUS_FUNC_TIME_MONITOR_DIFF(
"Tempus::IntegratorPseudoTransientAdjointSensitivity::advanceTime::state", TEMPUS_PTAS_AT_FWD);
138 state_status = state_integrator_->advanceTime(timeFinal);
141 bool sens_status =
true;
142 if (do_adjoint_integration_) {
143 TEMPUS_FUNC_TIME_MONITOR_DIFF(
"Tempus::IntegratorPseudoTransientAdjointSensitivity::advanceTime::adjoint", TEMPUS_PTAS_AT_ADJ);
149 sens_model_->setFinalTime(state_integrator_->getTime());
152 sens_model_->setForwardSolutionHistory(
153 state_integrator_->getSolutionHistory());
157 sens_status = sens_integrator_->advanceTime(timeFinal);
161 MEB::InArgs<Scalar> inargs = sens_model_->getNominalValues();
162 MEB::OutArgs<Scalar> outargs = sens_model_->createOutArgs();
163 inargs.set_t(sens_integrator_->getTime());
164 inargs.set_x(sens_integrator_->getX());
165 if (inargs.supports(MEB::IN_ARG_x_dot))
166 inargs.set_x_dot(sens_integrator_->getXDot());
167 if (inargs.supports(MEB::IN_ARG_x_dot_dot))
168 inargs.set_x_dot_dot(sens_integrator_->getXDotDot());
169 RCP<VectorBase<Scalar> > G = dgdp_;
170 if (G == Teuchos::null) {
171 G = Thyra::createMember(sens_model_->get_g_space(0));
172 dgdp_ = Teuchos::rcp_dynamic_cast<
DMVPV>(G);
174 if (g_ == Teuchos::null)
175 g_ = Thyra::createMember(sens_model_->get_g_space(1));
177 outargs.set_g(1, g_);
178 sens_model_->evalModel(inargs, outargs);
180 buildSolutionHistory();
183 return state_status && sens_status;
186 template<
class Scalar>
191 return solutionHistory_->getCurrentTime();
194 template<
class Scalar>
199 return solutionHistory_->getCurrentIndex();
202 template<
class Scalar>
207 Status state_status = state_integrator_->getStatus();
208 Status sens_status = sens_integrator_->getStatus();
216 template <
class Scalar>
219 state_integrator_->setStatus(st);
220 sens_integrator_->setStatus(st);
223 template<
class Scalar>
228 return state_integrator_->getStepper();
231 template<
class Scalar>
236 return state_integrator_->getStepper();
239 template<
class Scalar>
244 return sens_integrator_->getStepper();
247 template<
class Scalar>
252 return solutionHistory_;
255 template<
class Scalar>
260 return state_integrator_->getSolutionHistory();
263 template<
class Scalar>
268 return sens_integrator_->getSolutionHistory();
271 template<
class Scalar>
276 return solutionHistory_;
279 template<
class Scalar>
284 return state_integrator_->getTimeStepControl();
287 template<
class Scalar>
292 return state_integrator_->getNonConstTimeStepControl();
295 template<
class Scalar>
300 return state_integrator_->getNonConstTimeStepControl();
303 template<
class Scalar>
308 return sens_integrator_->getNonConstTimeStepControl();
311 template<
class Scalar>
316 return state_integrator_->getObserver();
319 template<
class Scalar>
324 state_integrator_->setObserver(obs);
325 sens_integrator_->setObserver(obs);
328 template<
class Scalar>
339 using Teuchos::rcp_dynamic_cast;
342 using Thyra::createMember;
347 RCP< const VectorSpaceBase<Scalar> > space = sens_model_->get_x_space();
348 RCP<DMVPV> Y = rcp_dynamic_cast<
DMVPV>(createMember(space));
349 RCP<DMVPV> Ydot = rcp_dynamic_cast<
DMVPV>(createMember(space));
350 RCP<DMVPV> Ydotdot = rcp_dynamic_cast<
DMVPV>(createMember(space));
354 if (y0 == Teuchos::null)
355 assign(Y->getNonconstMultiVector().ptr(), zero);
357 assign(Y->getNonconstMultiVector().ptr(), *y0);
360 if (ydot0 == Teuchos::null)
361 assign(Ydot->getNonconstMultiVector().ptr(), zero);
363 assign(Ydot->getNonconstMultiVector().ptr(), *ydot0);
366 if (ydotdot0 == Teuchos::null)
367 assign(Ydotdot->getNonconstMultiVector().ptr(), zero);
369 assign(Ydotdot->getNonconstMultiVector().ptr(), *ydotdot0);
371 state_integrator_->initializeSolutionHistory(t0, x0, xdot0, xdotdot0);
372 sens_integrator_->initializeSolutionHistory(t0, Y, Ydot, Ydotdot);
375 template<
class Scalar>
380 return state_integrator_->getX();
383 template<
class Scalar>
388 return state_integrator_->getXDot();
391 template<
class Scalar>
396 return state_integrator_->getXDotDot();
399 template<
class Scalar>
405 using Teuchos::rcp_dynamic_cast;
406 RCP<const DMVPV> mvpv =
407 rcp_dynamic_cast<
const DMVPV>(sens_integrator_->getX());
411 template<
class Scalar>
417 using Teuchos::rcp_dynamic_cast;
418 RCP<const DMVPV> mvpv =
419 rcp_dynamic_cast<
const DMVPV>(sens_integrator_->getXDot());
423 template<
class Scalar>
429 using Teuchos::rcp_dynamic_cast;
430 RCP<const DMVPV> mvpv =
431 rcp_dynamic_cast<
const DMVPV>(sens_integrator_->getXDotDot());
435 template<
class Scalar>
443 template<
class Scalar>
448 return dgdp_->getMultiVector();
451 template<
class Scalar>
456 std::string name =
"Tempus::IntegratorPseudoTransientAdjointSensitivity";
460 template<
class Scalar>
467 auto l_out = Teuchos::fancyOStream( out.
getOStream() );
469 l_out->setOutputToRootOnly(0);
471 *l_out << description() <<
"::describe" << std::endl;
472 state_integrator_->describe(*l_out, verbLevel);
473 sens_integrator_->describe(*l_out, verbLevel);
476 template<
class Scalar>
485 state_integrator_->getStepper()->getModel());
486 auto tmp_state_integrator = createIntegratorBasic<Scalar>(inputPL, model);
487 state_integrator_->copy(tmp_state_integrator);
490 sens_integrator_->getStepper()->getModel());
491 auto tmp_sens_integrator = createIntegratorBasic<Scalar>(inputPL, model);
492 sens_integrator_->copy(tmp_sens_integrator);
495 template<
class Scalar>
503 auto tmp_state_integrator = createIntegratorBasic<Scalar>();
504 auto model = state_integrator_->getStepper()->getModel();
505 tmp_state_integrator->setModel(model);
506 state_integrator_->copy(tmp_state_integrator);
508 auto tmp_sens_integrator = createIntegratorBasic<Scalar>();
509 model = sens_integrator_->getStepper()->getModel();
510 tmp_sens_integrator->setModel(model);
511 sens_integrator_->copy(tmp_sens_integrator);
514 sens_integrator_->getValidParameters());
518 template<
class Scalar>
526 state_integrator_->getValidParameters();
535 template<
class Scalar>
541 state_integrator_->getValidParameters());
545 template<
class Scalar>
553 template <
class Scalar>
565 if (inputPL != Teuchos::null) {
566 *pl = inputPL->
sublist(
"Sensitivities");
568 const Scalar tinit = state_integrator_->getTimeStepControl()->getInitTime();
569 const Scalar tfinal = state_integrator_->getTimeStepControl()->getFinalTime();
571 model, adjoint_residual_model, adjoint_solve_model_,
572 tinit, tfinal,
true, pl));
575 template<
class Scalar>
582 using Teuchos::rcp_dynamic_cast;
587 using Thyra::defaultProductVector;
594 state_integrator_->getSolutionHistory()->getValidParameters());
595 solutionHistory_ = createSolutionHistoryPL<Scalar>(shPL);
597 RCP<const VectorSpaceBase<Scalar> > x_space =
598 model_->get_x_space();
599 RCP<const VectorSpaceBase<Scalar> > adjoint_space =
600 sens_model_->get_x_space();
603 spaces[1] = adjoint_space;
604 RCP<const DPVS > prod_space = Thyra::productVectorSpace(spaces());
607 RCP<const SolutionHistory<Scalar> > state_solution_history =
608 state_integrator_->getSolutionHistory();
609 int num_states = state_solution_history->getNumStates();
610 for (
int i=0; i<num_states; ++i) {
611 RCP<const SolutionState<Scalar> > state = (*state_solution_history)[i];
614 RCP<DPV> x = defaultProductVector(prod_space);
615 assign(x->getNonconstVectorBlock(0).ptr(), *(state->getX()));
616 assign(x->getNonconstVectorBlock(1).ptr(), zero);
617 RCP<VectorBase<Scalar> > x_b = x;
620 RCP<VectorBase<Scalar> > x_dot_b;
621 if (state->getXDot() != Teuchos::null) {
622 RCP<DPV> x_dot = defaultProductVector(prod_space);
623 assign(x_dot->getNonconstVectorBlock(0).ptr(), *(state->getXDot()));
624 assign(x_dot->getNonconstVectorBlock(1).ptr(), zero);
629 RCP<VectorBase<Scalar> > x_dot_dot_b;
630 if (state->getXDotDot() != Teuchos::null) {
631 RCP<DPV> x_dot_dot = defaultProductVector(prod_space);
632 assign(x_dot_dot->getNonconstVectorBlock(0).ptr(),*(state->getXDotDot()));
633 assign(x_dot_dot->getNonconstVectorBlock(1).ptr(), zero);
634 x_dot_dot_b = x_dot_dot;
637 RCP<SolutionState<Scalar> > prod_state = state->clone();
638 prod_state->setX(x_b);
639 prod_state->setXDot(x_dot_b);
640 prod_state->setXDotDot(x_dot_dot_b);
641 solutionHistory_->addState(prod_state);
644 RCP<const VectorBase<Scalar> > frozen_x =
645 state_solution_history->getCurrentState()->getX();
646 RCP<const VectorBase<Scalar> > frozen_x_dot =
647 state_solution_history->getCurrentState()->getXDot();
648 RCP<const VectorBase<Scalar> > frozen_x_dot_dot =
649 state_solution_history->getCurrentState()->getXDotDot();
650 RCP<const SolutionHistory<Scalar> > sens_solution_history =
651 sens_integrator_->getSolutionHistory();
652 num_states = sens_solution_history->getNumStates();
653 for (
int i=0; i<num_states; ++i) {
654 RCP<const SolutionState<Scalar> > state = (*sens_solution_history)[i];
657 RCP<DPV> x = defaultProductVector(prod_space);
658 assign(x->getNonconstVectorBlock(0).ptr(), *frozen_x);
659 assign(x->getNonconstVectorBlock(1).ptr(), *(state->getX()));
660 RCP<VectorBase<Scalar> > x_b = x;
663 RCP<VectorBase<Scalar> > x_dot_b;
664 if (state->getXDot() != Teuchos::null) {
665 RCP<DPV> x_dot = defaultProductVector(prod_space);
666 assign(x_dot->getNonconstVectorBlock(0).ptr(), *frozen_x_dot);
667 assign(x_dot->getNonconstVectorBlock(1).ptr(), *(state->getXDot()));
672 RCP<VectorBase<Scalar> > x_dot_dot_b;
673 if (state->getXDotDot() != Teuchos::null) {
674 RCP<DPV> x_dot_dot = defaultProductVector(prod_space);
675 assign(x_dot_dot->getNonconstVectorBlock(0).ptr(), *frozen_x_dot_dot);
676 assign(x_dot_dot->getNonconstVectorBlock(1).ptr(),*(state->getXDotDot()));
677 x_dot_dot_b = x_dot_dot;
680 RCP<SolutionState<Scalar> > prod_state = state->clone();
681 prod_state->setX(x_b);
682 prod_state->setXDot(x_dot_b);
683 prod_state->setXDotDot(x_dot_dot_b);
684 solutionHistory_->addState(prod_state,
false);
689 template<
class Scalar>
701 template<
class Scalar>
705 std::string stepperType)
713 template<
class Scalar>
726 template<
class Scalar>
731 std::string stepperType)
739 template<
class Scalar>
753 template<
class Scalar>
759 std::string stepperType)
767 template<
class Scalar>
777 #endif // Tempus_IntegratorPseudoTransientAdjointSensitivity_impl_hpp
Teuchos::RCP< Stepper< Scalar > > getStateStepper() const
static Teuchos::RCP< const Teuchos::ParameterList > getValidParameters()
Teuchos::RCP< TimeStepControl< Scalar > > getSensNonConstTimeStepControl()
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &pl) override
virtual Scalar getTime() const override
Get current time.
Teuchos::RCP< const SolutionHistory< Scalar > > getStateSolutionHistory() const
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const override
virtual Teuchos::RCP< IntegratorObserver< Scalar > > getObserver()
Get the Observer.
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getY() const
Get the current adjoint solution, y.
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getXDot() const
Get the current time derivative of the solution, xdot.
Teuchos::RCP< Tempus::IntegratorPseudoTransientAdjointSensitivity< Scalar > > integratorPseudoTransientAdjointSensitivity(Teuchos::RCP< Teuchos::ParameterList > pList, const Teuchos::RCP< Thyra::ModelEvaluator< Scalar > > &model)
Nonmember constructor.
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList() override
virtual void setObserver(Teuchos::RCP< IntegratorObserver< Scalar > > obs=Teuchos::null)
Set the Observer.
std::string description() const override
virtual Teuchos::RCP< SolutionHistory< Scalar > > getNonConstSolutionHistory() override
Get the SolutionHistory.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const override
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getDgDp() const
Return adjoint sensitivity stored in gradient format.
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getYDot() const
Get the current time derivative of the adjoint solution, ydot.
Time integrator suitable for pseudotransient adjoint sensitivity analysis.
SensitivityStepMode getStepMode() const
What mode the current time integration step is in.
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getG() const
Return response function g.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual Teuchos::RCP< TimeStepControl< Scalar > > getNonConstTimeStepControl() override
Status
Status for the Integrator, the Stepper and the SolutionState.
IntegratorObserver class for time integrators.
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getXDotDot() const
Get the current second time derivative of the solution, xdotdot.
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 > > y0=Teuchos::null, Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > ydot0=Teuchos::null, Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > ydotdot0=Teuchos::null)
Set the initial state from Thyra::VectorBase(s)
RCP< const MultiVectorBase< Scalar > > getMultiVector() const
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getYDotDot() const
Get the current second time derivative of the adjoint solution, ydotdot.
IntegratorPseudoTransientAdjointSensitivity()
Destructor.
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList() override
ParameterList & setParameters(const ParameterList &source)
RCP< ImplicitAdjointModelEvaluator< Scalar > > implicitAdjointModelEvaluator(const RCP< const ModelEvaluator< Scalar > > &model)
Teuchos::RCP< Stepper< Scalar > > getSensStepper() const
Teuchos::RCP< TimeStepControl< Scalar > > getStateNonConstTimeStepControl()
RCP< std::basic_ostream< char_type, traits_type > > getOStream()
virtual void setStatus(const Status st) override
Set Status.
virtual int getIndex() const override
Get current index.
virtual Teuchos::RCP< Stepper< Scalar > > getStepper() const override
Get the Stepper.
virtual Teuchos::RCP< const SolutionHistory< Scalar > > getSolutionHistory() const override
Get the SolutionHistory.
virtual Teuchos::RCP< const TimeStepControl< Scalar > > getTimeStepControl() const override
Get the TimeStepControl.
ModelEvaluator for forming adjoint sensitivity equations.
Teuchos::RCP< const SolutionHistory< Scalar > > getSensSolutionHistory() const
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
virtual bool advanceTime()
Advance the solution to timeMax, and return true if successful.
void buildSolutionHistory()
Teuchos::RCP< AdjointSensitivityModelEvaluator< Scalar > > createSensitivityModel(const Teuchos::RCP< Thyra::ModelEvaluator< Scalar > > &model, const Teuchos::RCP< Thyra::ModelEvaluator< Scalar > > &adjoint_residual_model, const Teuchos::RCP< Thyra::ModelEvaluator< Scalar > > &adjoint_solve_model, const Teuchos::RCP< Teuchos::ParameterList > &inputPL)
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getX() const
Get the current solution, x.
virtual Status getStatus() const override
Get Status.