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"
21 template <
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)
39 template<
class Scalar>
43 force_W_update_(false),
50 template<
class Scalar>
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;
81 template<
class Scalar>
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;
120 template<
class Scalar>
125 return solutionHistory_->getCurrentTime();
128 template<
class Scalar>
133 return solutionHistory_->getCurrentIndex();
136 template<
class Scalar>
141 Status state_status = state_integrator_->getStatus();
142 Status sens_status = sens_integrator_->getStatus();
150 template <
class Scalar>
153 state_integrator_->setStatus(st);
154 sens_integrator_->setStatus(st);
157 template<
class Scalar>
162 return state_integrator_->getStepper();
165 template<
class Scalar>
170 return state_integrator_->getStepper();
173 template<
class Scalar>
178 return sens_integrator_->getStepper();
181 template<
class Scalar>
186 return solutionHistory_;
189 template<
class Scalar>
194 return state_integrator_->getSolutionHistory();
197 template<
class Scalar>
202 return sens_integrator_->getSolutionHistory();
205 template<
class Scalar>
210 return solutionHistory_;
213 template<
class Scalar>
218 return state_integrator_->getTimeStepControl();
221 template<
class Scalar>
226 return state_integrator_->getNonConstTimeStepControl();
229 template<
class Scalar>
234 return state_integrator_->getNonConstTimeStepControl();
237 template<
class Scalar>
242 return sens_integrator_->getNonConstTimeStepControl();
245 template<
class Scalar>
250 return state_integrator_->getObserver();
253 template<
class Scalar>
258 state_integrator_->setObserver(obs);
259 sens_integrator_->setObserver(obs);
262 template<
class Scalar>
273 using Teuchos::rcp_dynamic_cast;
276 using Thyra::createMember;
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));
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);
310 template<
class Scalar>
315 return state_integrator_->getX();
318 template<
class Scalar>
324 using Teuchos::rcp_dynamic_cast;
328 rcp_dynamic_cast<
const DMVPV>(sens_integrator_->getX());
332 template<
class Scalar>
337 return state_integrator_->getXDot();
340 template<
class Scalar>
346 using Teuchos::rcp_dynamic_cast;
349 RCP<const DMVPV> Xdot =
350 rcp_dynamic_cast<
const DMVPV>(sens_integrator_->getXDot());
354 template<
class Scalar>
359 return state_integrator_->getXDotDot();
362 template<
class Scalar>
368 using Teuchos::rcp_dynamic_cast;
371 RCP<const DMVPV> Xdotdot =
372 rcp_dynamic_cast<
const DMVPV>(sens_integrator_->getXDotDot());
376 template<
class Scalar>
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());
395 Thyra::createMember(sens_model_->get_g_space(1));
398 sens_model_->evalModel(inargs, outargs);
402 template<
class Scalar>
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());
422 Thyra::createMember(sens_model_->get_g_space(0));
426 sens_model_->evalModel(inargs, outargs);
427 return dgdp->getMultiVector();
430 template<
class Scalar>
435 std::string name =
"Tempus::IntegratorPseudoTransientForwardSensitivity";
439 template<
class Scalar>
446 auto l_out = Teuchos::fancyOStream( out.
getOStream() );
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);
455 template<
class Scalar>
463 template<
class Scalar>
470 using Teuchos::rcp_dynamic_cast;
475 using Thyra::createMembers;
476 using Thyra::multiVectorProductVector;
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);
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);
590 template<
class Scalar>
599 auto fwd_integrator = createIntegratorBasic<Scalar>(pList, model);
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);
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);
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);
642 template<
class Scalar>
652 #endif // Tempus_IntegratorPseudoTransientForwardSensitivity_impl_hpp
virtual Teuchos::RCP< Stepper< Scalar > > getStepper() const override
Get the Stepper.
virtual void setObserver(Teuchos::RCP< IntegratorObserver< Scalar > > obs=Teuchos::null)
Set the Observer.
virtual void setStatus(const Status st) override
Set Status.
Teuchos::RCP< const SolutionHistory< Scalar > > getSensSolutionHistory() const
virtual Status getStatus() const override
Get Status.
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(Teuchos::RCP< Teuchos::ParameterList > pList, const Teuchos::RCP< Thyra::ModelEvaluator< Scalar > > &model, const Teuchos::RCP< Thyra::ModelEvaluator< Scalar > > &sens_residual_model, const Teuchos::RCP< Thyra::ModelEvaluator< Scalar > > &sens_solve_model)
Nonmember constructor.
T & get(const std::string &name, T def_value)
virtual Teuchos::RCP< const SolutionHistory< Scalar > > getSolutionHistory() const override
Get the SolutionHistory.
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getDxDp() const
virtual int getIndex() const override
Get current index.
Teuchos::RCP< TimeStepControl< Scalar > > getStateNonConstTimeStepControl()
Teuchos::RCP< TimeStepControl< Scalar > > getSensNonConstTimeStepControl()
virtual Teuchos::RCP< TimeStepControl< Scalar > > getNonConstTimeStepControl() override
static Teuchos::RCP< const Teuchos::ParameterList > getValidParameters()
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getG() const
Return response function g.
Teuchos::RCP< IntegratorBasic< Scalar > > state_integrator_
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 bool advanceTime()
Advance the solution to timeMax, and return true if successful.
bool remove(std::string const &name, bool throwIfNotExists=true)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
SensitivityStepMode getStepMode() const
What mode the current time integration step is in.
Status
Status for the Integrator, the Stepper and the SolutionState.
Teuchos::RCP< const SolutionHistory< Scalar > > getStateSolutionHistory() const
IntegratorObserver class for time integrators.
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getDgDp() const
Return forward sensitivity stored in Jacobian format.
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getXDotDot() const
Get current the second time derivative of the solution, xdotdot.
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const override
std::string description() const override
RCP< const MultiVectorBase< Scalar > > getMultiVector() const
virtual Teuchos::RCP< IntegratorObserver< Scalar > > getObserver()
Get the Observer.
virtual Scalar getTime() const override
Get current time.
ParameterList & setParameters(const ParameterList &source)
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getDXDotDp() const
RCP< std::basic_ostream< char_type, traits_type > > getOStream()
IntegratorPseudoTransientForwardSensitivity()
Destructor.
Time integrator suitable for pseudotransient forward sensitivity analysis.
ParameterList & setParametersNotAlreadySet(const ParameterList &source)
virtual RCP< const VectorSpaceBase< Scalar > > domain() const
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
Teuchos::RCP< Stepper< Scalar > > getSensStepper() const
void buildSolutionHistory()
virtual Teuchos::RCP< const TimeStepControl< Scalar > > getTimeStepControl() const override
Get the TimeStepControl.
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getDXDotDotDp() const
virtual Teuchos::RCP< SolutionHistory< Scalar > > getNonConstSolutionHistory() override
Get the SolutionHistory.
Teuchos::RCP< IntegratorBasic< Scalar > > sens_integrator_
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getXDot() const
Get current the time derivative of the solution, xdot.
A ModelEvaluator decorator for sensitivity analysis.
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getX() const
Get current the solution, x.
Teuchos::RCP< Stepper< Scalar > > getStateStepper() const