9#ifndef Tempus_IntegratorPseudoTransientForwardSensitivity_impl_hpp
10#define Tempus_IntegratorPseudoTransientForwardSensitivity_impl_hpp
12#include "Thyra_DefaultMultiVectorProductVector.hpp"
13#include "Thyra_VectorStdOps.hpp"
14#include "Thyra_MultiVectorStdOps.hpp"
21template <
class Scalar>
27 const bool reuse_solver,
28 const bool force_W_update)
30 , sens_model_(sens_model)
31 , state_integrator_(fwd_integrator)
32 , sens_integrator_(sens_integrator)
33 , reuse_solver_(reuse_solver)
34 , force_W_update_(force_W_update)
43 force_W_update_(false),
60 bool state_status = state_integrator_->advanceTime();
63 sens_model_->setForwardSolutionState(state_integrator_->getCurrentState());
67 state_integrator_->getStepper()->getSolver() != Teuchos::null) {
68 sens_model_->setSolver(state_integrator_->getStepper()->getSolver(),
74 bool sens_status = sens_integrator_->advanceTime();
76 buildSolutionHistory();
78 return state_status && sens_status;
86 TEMPUS_FUNC_TIME_MONITOR_DIFF(
"Tempus::IntegratorPseudoTransientForwardSensitivity::advanceTime()", TEMPUS_PTFS_AT);
92 bool state_status =
true;
94 TEMPUS_FUNC_TIME_MONITOR_DIFF(
"Tempus::IntegratorPseudoTransientForwardSensitivity::advanceTime::state", TEMPUS_PTFS_AT_FWD);
95 state_status = state_integrator_->advanceTime(timeFinal);
99 sens_model_->setForwardSolutionState(state_integrator_->getCurrentState());
103 state_integrator_->getStepper()->getSolver() != Teuchos::null) {
104 sens_model_->setSolver(state_integrator_->getStepper()->getSolver(),
109 bool sens_status =
true;
111 TEMPUS_FUNC_TIME_MONITOR_DIFF(
"Tempus::IntegratorPseudoTransientForwardSensitivity::advanceTime::sensitivity", TEMPUS_PTFS_AT_SEN);
112 sens_status = sens_integrator_->advanceTime(timeFinal);
115 buildSolutionHistory();
117 return state_status && sens_status;
120template<
class Scalar>
125 return solutionHistory_->getCurrentTime();
128template<
class Scalar>
133 return solutionHistory_->getCurrentIndex();
136template<
class Scalar>
141 Status state_status = state_integrator_->getStatus();
142 Status sens_status = sens_integrator_->getStatus();
150template <
class Scalar>
153 state_integrator_->setStatus(st);
154 sens_integrator_->setStatus(st);
157template<
class Scalar>
158Teuchos::RCP<Stepper<Scalar> >
162 return state_integrator_->getStepper();
165template<
class Scalar>
166Teuchos::RCP<Stepper<Scalar> >
170 return state_integrator_->getStepper();
173template<
class Scalar>
174Teuchos::RCP<Stepper<Scalar> >
178 return sens_integrator_->getStepper();
181template<
class Scalar>
182Teuchos::RCP<const SolutionHistory<Scalar> >
186 return solutionHistory_;
189template<
class Scalar>
190Teuchos::RCP<const SolutionHistory<Scalar> >
194 return state_integrator_->getSolutionHistory();
197template<
class Scalar>
198Teuchos::RCP<const SolutionHistory<Scalar> >
202 return sens_integrator_->getSolutionHistory();
205template<
class Scalar>
206Teuchos::RCP<SolutionHistory<Scalar> >
210 return solutionHistory_;
213template<
class Scalar>
214Teuchos::RCP<const TimeStepControl<Scalar> >
218 return state_integrator_->getTimeStepControl();
221template<
class Scalar>
222Teuchos::RCP<TimeStepControl<Scalar> >
226 return state_integrator_->getNonConstTimeStepControl();
229template<
class Scalar>
230Teuchos::RCP<TimeStepControl<Scalar> >
234 return state_integrator_->getNonConstTimeStepControl();
237template<
class Scalar>
238Teuchos::RCP<TimeStepControl<Scalar> >
242 return sens_integrator_->getNonConstTimeStepControl();
245template<
class Scalar>
246Teuchos::RCP<IntegratorObserver<Scalar> >
250 return state_integrator_->getObserver();
253template<
class Scalar>
258 state_integrator_->setObserver(obs);
259 sens_integrator_->setObserver(obs);
262template<
class Scalar>
273 using Teuchos::rcp_dynamic_cast;
274 using Thyra::VectorSpaceBase;
276 using Thyra::createMember;
277 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
282 RCP< const VectorSpaceBase<Scalar> > space = sens_model_->get_x_space();
283 RCP<DMVPV> X = rcp_dynamic_cast<DMVPV>(createMember(space));
284 RCP<DMVPV> Xdot = rcp_dynamic_cast<DMVPV>(createMember(space));
285 RCP<DMVPV> Xdotdot = rcp_dynamic_cast<DMVPV>(createMember(space));
286 const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
289 if (DxDp0 == Teuchos::null)
290 assign(X->getNonconstMultiVector().ptr(), zero);
292 assign(X->getNonconstMultiVector().ptr(), *DxDp0);
295 if (DxdotDp0 == Teuchos::null)
296 assign(Xdot->getNonconstMultiVector().ptr(), zero);
298 assign(Xdot->getNonconstMultiVector().ptr(), *DxdotDp0);
301 if (DxdotDp0 == Teuchos::null)
302 assign(Xdotdot->getNonconstMultiVector().ptr(), zero);
304 assign(Xdotdot->getNonconstMultiVector().ptr(), *DxdotdotDp0);
306 state_integrator_->initializeSolutionHistory(t0, x0, xdot0, xdotdot0);
307 sens_integrator_->initializeSolutionHistory(t0, X, Xdot, Xdotdot);
310template<
class Scalar>
311Teuchos::RCP<const Thyra::VectorBase<Scalar> >
315 return state_integrator_->getX();
318template<
class Scalar>
319Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
324 using Teuchos::rcp_dynamic_cast;
325 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
328 rcp_dynamic_cast<const DMVPV>(sens_integrator_->getX());
329 return X->getMultiVector();
332template<
class Scalar>
333Teuchos::RCP<const Thyra::VectorBase<Scalar> >
337 return state_integrator_->getXDot();
340template<
class Scalar>
341Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
346 using Teuchos::rcp_dynamic_cast;
347 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
349 RCP<const DMVPV> Xdot =
350 rcp_dynamic_cast<const DMVPV>(sens_integrator_->getXDot());
351 return Xdot->getMultiVector();
354template<
class Scalar>
355Teuchos::RCP<const Thyra::VectorBase<Scalar> >
359 return state_integrator_->getXDotDot();
362template<
class Scalar>
363Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
368 using Teuchos::rcp_dynamic_cast;
369 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
371 RCP<const DMVPV> Xdotdot =
372 rcp_dynamic_cast<const DMVPV>(sens_integrator_->getXDotDot());
373 return Xdotdot->getMultiVector();
376template<
class Scalar>
377Teuchos::RCP<const Thyra::VectorBase<Scalar> >
381 typedef Thyra::ModelEvaluatorBase MEB;
385 MEB::InArgs<Scalar> inargs = sens_model_->getNominalValues();
386 MEB::OutArgs<Scalar> outargs = sens_model_->createOutArgs();
387 inargs.set_t(sens_integrator_->getTime());
388 inargs.set_x(sens_integrator_->getX());
389 if (inargs.supports(MEB::IN_ARG_x_dot))
390 inargs.set_x_dot(sens_integrator_->getXDot());
391 if (inargs.supports(MEB::IN_ARG_x_dot_dot))
392 inargs.set_x_dot_dot(sens_integrator_->getXDotDot());
394 Teuchos::RCP<Thyra::VectorBase<Scalar> > g =
395 Thyra::createMember(sens_model_->get_g_space(1));
398 sens_model_->evalModel(inargs, outargs);
402template<
class Scalar>
403Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
407 typedef Thyra::ModelEvaluatorBase MEB;
408 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
412 MEB::InArgs<Scalar> inargs = sens_model_->getNominalValues();
413 MEB::OutArgs<Scalar> outargs = sens_model_->createOutArgs();
414 inargs.set_t(sens_integrator_->getTime());
415 inargs.set_x(sens_integrator_->getX());
416 if (inargs.supports(MEB::IN_ARG_x_dot))
417 inargs.set_x_dot(sens_integrator_->getXDot());
418 if (inargs.supports(MEB::IN_ARG_x_dot_dot))
419 inargs.set_x_dot_dot(sens_integrator_->getXDotDot());
421 Teuchos::RCP<Thyra::VectorBase<Scalar> > G =
422 Thyra::createMember(sens_model_->get_g_space(0));
423 Teuchos::RCP<DMVPV> dgdp = Teuchos::rcp_dynamic_cast<DMVPV>(G);
426 sens_model_->evalModel(inargs, outargs);
427 return dgdp->getMultiVector();
430template<
class Scalar>
435 std::string name =
"Tempus::IntegratorPseudoTransientForwardSensitivity";
439template<
class Scalar>
443 Teuchos::FancyOStream &out,
444 const Teuchos::EVerbosityLevel verbLevel)
const
446 auto l_out = Teuchos::fancyOStream( out.getOStream() );
447 Teuchos::OSTab ostab(*l_out, 2, this->description());
448 l_out->setOutputToRootOnly(0);
450 *l_out << description() <<
"::describe" << std::endl;
451 state_integrator_->describe(*l_out, verbLevel);
452 sens_integrator_->describe(*l_out, verbLevel);
455template<
class Scalar>
463template<
class Scalar>
470 using Teuchos::rcp_dynamic_cast;
471 using Teuchos::ParameterList;
474 using Thyra::VectorSpaceBase;
475 using Thyra::createMembers;
476 using Thyra::multiVectorProductVector;
478 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
484 RCP<ParameterList> shPL;
486 solutionHistory_ = createSolutionHistoryPL<Scalar>(shPL);
488 const int num_param =
489 rcp_dynamic_cast<const DMVPV>(sens_integrator_->getX())->getMultiVector()->domain()->dim();
490 RCP<const VectorSpaceBase<Scalar> > x_space = model_->get_x_space();
491 RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> > prod_space =
492 Thyra::multiVectorProductVectorSpace(x_space, num_param+1);
493 const Teuchos::Range1D rng(1,num_param);
494 const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
496 RCP<const SolutionHistory<Scalar> > state_solution_history =
497 state_integrator_->getSolutionHistory();
498 int num_states = state_solution_history->getNumStates();
499 for (
int i=0; i<num_states; ++i) {
500 RCP<const SolutionState<Scalar> > state = (*state_solution_history)[i];
503 RCP< MultiVectorBase<Scalar> > x_mv =
504 createMembers(x_space, num_param+1);
505 assign(x_mv->col(0).ptr(), *(state->getX()));
506 assign(x_mv->subView(rng).ptr(), zero);
507 RCP<VectorBase<Scalar> > x = multiVectorProductVector(prod_space, x_mv);
510 RCP<VectorBase<Scalar> > x_dot;
511 if (state->getXDot() != Teuchos::null) {
512 RCP< MultiVectorBase<Scalar> > x_dot_mv =
513 createMembers(x_space, num_param+1);
514 assign(x_dot_mv->col(0).ptr(), *(state->getXDot()));
515 assign(x_dot_mv->subView(rng).ptr(), zero);
516 x_dot = multiVectorProductVector(prod_space, x_dot_mv);
520 RCP<VectorBase<Scalar> > x_dot_dot;
521 if (state->getXDotDot() != Teuchos::null) {
522 RCP< MultiVectorBase<Scalar> > x_dot_dot_mv =
523 createMembers(x_space, num_param+1);
524 assign(x_dot_dot_mv->col(0).ptr(), *(state->getXDotDot()));
525 assign(x_dot_dot_mv->subView(rng).ptr(), zero);
526 x_dot_dot = multiVectorProductVector(prod_space, x_dot_dot_mv);
529 RCP<SolutionState<Scalar> > prod_state = state->clone();
531 prod_state->setXDot(x_dot);
532 prod_state->setXDotDot(x_dot_dot);
533 solutionHistory_->addState(prod_state);
536 RCP<const VectorBase<Scalar> > frozen_x =
537 state_solution_history->getCurrentState()->getX();
538 RCP<const VectorBase<Scalar> > frozen_x_dot =
539 state_solution_history->getCurrentState()->getXDot();
540 RCP<const VectorBase<Scalar> > frozen_x_dot_dot =
541 state_solution_history->getCurrentState()->getXDotDot();
542 RCP<const SolutionHistory<Scalar> > sens_solution_history =
543 sens_integrator_->getSolutionHistory();
544 num_states = sens_solution_history->getNumStates();
545 for (
int i=0; i<num_states; ++i) {
546 RCP<const SolutionState<Scalar> > state = (*sens_solution_history)[i];
549 RCP< MultiVectorBase<Scalar> > x_mv =
550 createMembers(x_space, num_param+1);
551 RCP<const MultiVectorBase<Scalar> > dxdp =
552 rcp_dynamic_cast<const DMVPV>(state->getX())->getMultiVector();
553 assign(x_mv->col(0).ptr(), *(frozen_x));
554 assign(x_mv->subView(rng).ptr(), *dxdp);
555 RCP<VectorBase<Scalar> > x = multiVectorProductVector(prod_space, x_mv);
558 RCP<VectorBase<Scalar> > x_dot;
559 if (state->getXDot() != Teuchos::null) {
560 RCP< MultiVectorBase<Scalar> > x_dot_mv =
561 createMembers(x_space, num_param+1);
562 RCP<const MultiVectorBase<Scalar> > dxdotdp =
563 rcp_dynamic_cast<const DMVPV>(state->getXDot())->getMultiVector();
564 assign(x_dot_mv->col(0).ptr(), *(frozen_x_dot));
565 assign(x_dot_mv->subView(rng).ptr(), *dxdotdp);
566 x_dot = multiVectorProductVector(prod_space, x_dot_mv);
570 RCP<VectorBase<Scalar> > x_dot_dot;
571 if (state->getXDotDot() != Teuchos::null) {
572 RCP< MultiVectorBase<Scalar> > x_dot_dot_mv =
573 createMembers(x_space, num_param+1);
574 RCP<const MultiVectorBase<Scalar> > dxdotdotdp =
575 rcp_dynamic_cast<const DMVPV>(state->getXDotDot())->getMultiVector();
576 assign(x_dot_dot_mv->col(0).ptr(), *(frozen_x_dot_dot));
577 assign(x_dot_dot_mv->subView(rng).ptr(), *dxdotdotdp);
578 x_dot_dot = multiVectorProductVector(prod_space, x_dot_dot_mv);
581 RCP<SolutionState<Scalar> > prod_state = state->clone();
583 prod_state->setXDot(x_dot);
584 prod_state->setXDotDot(x_dot_dot);
585 solutionHistory_->addState(prod_state,
false);
590template<
class Scalar>
591Teuchos::RCP<Tempus::IntegratorPseudoTransientForwardSensitivity<Scalar> >
593 Teuchos::RCP<Teuchos::ParameterList> pList,
599 auto fwd_integrator = createIntegratorBasic<Scalar>(pList, model);
600 Teuchos::RCP<SensitivityModelEvaluatorBase<Scalar> > sens_model;
601 Teuchos::RCP<IntegratorBasic<Scalar> > sens_integrator;
604 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::rcp(
new Teuchos::ParameterList);
605 Teuchos::RCP<const Teuchos::ParameterList> integrator_pl = fwd_integrator->getValidParameters();
606 Teuchos::RCP<const Teuchos::ParameterList> sensitivity_pl =
608 pl->setParameters(*integrator_pl);
609 pl->sublist(
"Sensitivities").setParameters(*sensitivity_pl);
610 pl->sublist(
"Sensitivities").set(
"Reuse State Linear Solver",
false);
611 pl->sublist(
"Sensitivities").set(
"Force W Update",
false);
612 pl->sublist(
"Sensitivities").set(
"Cache Matrices",
false);
613 pList->setParametersNotAlreadySet(*pl);
616 bool reuse_solver = pList->sublist(
"Sensitivities").get(
"Reuse State Linear Solver",
false);
617 bool force_W_update = pList->sublist(
"Sensitivities").get(
"Force W Update",
false);
618 bool cache_matrices = pList->sublist(
"Sensitivities").get(
"Cache Matrices",
false);
621 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
622 if (pList!= Teuchos::null)
624 *pl = pList->sublist(
"Sensitivities");
626 pl->remove(
"Reuse State Linear Solver");
627 pl->remove(
"Force W Update");
628 pl->remove(
"Cache Matrices");
630 model, sens_residual_model, sens_solve_model, cache_matrices, pl);
631 sens_integrator = createIntegratorBasic<Scalar>(pList, sens_model);
634 Teuchos::RCP<Tempus::IntegratorPseudoTransientForwardSensitivity<Scalar>> integrator =
642template<
class Scalar>
643Teuchos::RCP<Tempus::IntegratorPseudoTransientForwardSensitivity<Scalar> >
646 Teuchos::RCP<Tempus::IntegratorPseudoTransientForwardSensitivity<Scalar> > integrator =
IntegratorObserver class for time integrators.
Time integrator suitable for pseudotransient forward sensitivity analysis.
Teuchos::RCP< const SolutionHistory< Scalar > > getSensSolutionHistory() const
virtual Status getStatus() const override
Get Status.
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getXDotDot() const
Get current the second time derivative of the solution, xdotdot.
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getG() const
Return response function g.
virtual Teuchos::RCP< TimeStepControl< Scalar > > getNonConstTimeStepControl() override
std::string description() const override
IntegratorPseudoTransientForwardSensitivity()
Destructor.
SensitivityStepMode getStepMode() const
What mode the current time integration step is in.
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 void setStatus(const Status st) override
Set Status.
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getX() const
Get current the solution, x.
virtual Teuchos::RCP< SolutionHistory< Scalar > > getNonConstSolutionHistory() override
Get the SolutionHistory.
Teuchos::RCP< Stepper< Scalar > > getSensStepper() const
virtual Teuchos::RCP< IntegratorObserver< Scalar > > getObserver()
Get the Observer.
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getDgDp() const
Return forward sensitivity stored in Jacobian format.
virtual int getIndex() const override
Get current index.
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getXDot() const
Get current the time derivative of the solution, xdot.
Teuchos::RCP< TimeStepControl< Scalar > > getStateNonConstTimeStepControl()
void buildSolutionHistory()
virtual Teuchos::RCP< const TimeStepControl< Scalar > > getTimeStepControl() const override
Get the TimeStepControl.
Teuchos::RCP< IntegratorBasic< Scalar > > state_integrator_
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getDXDotDp() const
Teuchos::RCP< IntegratorBasic< Scalar > > sens_integrator_
Teuchos::RCP< const SolutionHistory< Scalar > > getStateSolutionHistory() const
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getDXDotDotDp() const
virtual Scalar getTime() const override
Get current time.
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const override
virtual Teuchos::RCP< Stepper< Scalar > > getStepper() const override
Get the Stepper.
virtual Teuchos::RCP< const SolutionHistory< Scalar > > getSolutionHistory() const override
Get the SolutionHistory.
virtual void setObserver(Teuchos::RCP< IntegratorObserver< Scalar > > obs=Teuchos::null)
Set the Observer.
virtual bool advanceTime()
Advance the solution to timeMax, and return true if successful.
Teuchos::RCP< TimeStepControl< Scalar > > getSensNonConstTimeStepControl()
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getDxDp() const
Teuchos::RCP< Stepper< Scalar > > getStateStepper() const
A ModelEvaluator decorator for sensitivity analysis.
static Teuchos::RCP< const Teuchos::ParameterList > getValidParameters()
Status
Status for the Integrator, the Stepper and the SolutionState.
Teuchos::RCP< SensitivityModelEvaluatorBase< Scalar > > wrapStaggeredFSAModelEvaluator(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)
Teuchos::RCP< Tempus::IntegratorPseudoTransientForwardSensitivity< Scalar > > createIntegratorPseudoTransientForwardSensitivity()
Nonmember constructor.