126 TEMPUS_FUNC_TIME_MONITOR_DIFF(
"Tempus::IntegratorPseudoTransientAdjointSensitivity::advanceTime()", TEMPUS_PTAS_AT);
130 typedef Thyra::ModelEvaluatorBase MEB;
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;
339 using Teuchos::rcp_dynamic_cast;
340 using Thyra::VectorSpaceBase;
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));
351 const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
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);
582 using Teuchos::rcp_dynamic_cast;
583 using Teuchos::ParameterList;
585 using Thyra::VectorSpaceBase;
587 using Thyra::defaultProductVector;
588 typedef Thyra::DefaultProductVector<Scalar> DPV;
589 typedef Thyra::DefaultProductVectorSpace<Scalar> DPVS;
593 auto shPL = Teuchos::rcp_const_cast<Teuchos::ParameterList> (
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();
601 Teuchos::Array< RCP<const VectorSpaceBase<Scalar> > > spaces(2);
603 spaces[1] = adjoint_space;
604 RCP<const DPVS > prod_space = Thyra::productVectorSpace(spaces());
605 const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
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);