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>
34 template<
class Scalar>
38 std::string stepperType) :
40 force_W_update_(false)
48 template<
class Scalar>
52 force_W_update_(false)
58 template<
class Scalar>
67 bool state_status = state_integrator_->advanceTime();
70 sens_model_->setForwardSolutionState(state_integrator_->getCurrentState());
74 state_integrator_->getStepper()->getSolver() != Teuchos::null) {
75 sens_model_->setSolver(state_integrator_->getStepper()->getSolver(),
80 bool sens_status = sens_integrator_->advanceTime();
82 buildSolutionHistory();
84 return state_status && sens_status;
87 template<
class Scalar>
96 bool 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 = sens_integrator_->advanceTime(timeFinal);
111 buildSolutionHistory();
113 return state_status && sens_status;
116 template<
class Scalar>
121 return solutionHistory_->getCurrentTime();
124 template<
class Scalar>
129 return solutionHistory_->getCurrentIndex();
132 template<
class Scalar>
137 Status state_status = state_integrator_->getStatus();
138 Status sens_status = sens_integrator_->getStatus();
146 template<
class Scalar>
151 return state_integrator_->getStepper();
154 template<
class Scalar>
159 return state_integrator_->getTempusParameterList();
162 template<
class Scalar>
167 state_integrator_->setTempusParameterList(pl);
168 sens_integrator_->setTempusParameterList(pl);
171 template<
class Scalar>
176 return solutionHistory_;
179 template<
class Scalar>
184 return state_integrator_->getTimeStepControl();
187 template<
class Scalar>
192 return state_integrator_->getNonConstTimeStepControl();
195 template<
class Scalar>
206 using Teuchos::rcp_dynamic_cast;
209 using Thyra::createMember;
215 RCP< const VectorSpaceBase<Scalar> > space = sens_model_->get_x_space();
216 RCP<DMVPV> X = rcp_dynamic_cast<DMVPV>(createMember(space));
217 RCP<DMVPV> Xdot = rcp_dynamic_cast<DMVPV>(createMember(space));
218 RCP<DMVPV> Xdotdot = rcp_dynamic_cast<DMVPV>(createMember(space));
222 if (DxDp0 == Teuchos::null)
223 assign(X->getNonconstMultiVector().ptr(), zero);
225 assign(X->getNonconstMultiVector().ptr(), *DxDp0);
228 if (DxdotDp0 == Teuchos::null)
229 assign(Xdot->getNonconstMultiVector().ptr(), zero);
231 assign(Xdot->getNonconstMultiVector().ptr(), *DxdotDp0);
234 if (DxdotDp0 == Teuchos::null)
235 assign(Xdotdot->getNonconstMultiVector().ptr(), zero);
237 assign(Xdotdot->getNonconstMultiVector().ptr(), *DxdotdotDp0);
239 state_integrator_->initializeSolutionHistory(t0, x0, xdot0, xdotdot0);
240 sens_integrator_->initializeSolutionHistory(t0, X, Xdot, Xdotdot);
243 template<
class Scalar>
249 using Teuchos::rcp_dynamic_cast;
253 rcp_dynamic_cast<
const DMVPV>(solutionHistory_->getCurrentState()->getX());
257 template<
class Scalar>
263 using Teuchos::rcp_dynamic_cast;
267 rcp_dynamic_cast<
const DMVPV>(solutionHistory_->getCurrentState()->getX());
270 return X->getMultiVector()->subView(rng);
273 template<
class Scalar>
279 using Teuchos::rcp_dynamic_cast;
282 RCP<const DMVPV> Xdot =
283 rcp_dynamic_cast<
const DMVPV>(solutionHistory_->getCurrentState()->getXDot());
287 template<
class Scalar>
293 using Teuchos::rcp_dynamic_cast;
296 RCP<const DMVPV> Xdot =
297 rcp_dynamic_cast<
const DMVPV>(solutionHistory_->getCurrentState()->getXDot());
300 return Xdot->getMultiVector()->subView(rng);
303 template<
class Scalar>
309 using Teuchos::rcp_dynamic_cast;
312 RCP<const DMVPV> Xdotdot =
313 rcp_dynamic_cast<
const DMVPV>(solutionHistory_->getCurrentState()->getXDotDot());
317 template<
class Scalar>
323 using Teuchos::rcp_dynamic_cast;
326 RCP<const DMVPV> Xdotdot =
327 rcp_dynamic_cast<
const DMVPV>(solutionHistory_->getCurrentState()->getXDotDot());
328 const int num_param = Xdotdot->
getMultiVector()->domain()->dim()-1;
330 return Xdotdot->getMultiVector()->subView(rng);
333 template<
class Scalar>
338 std::string name =
"Tempus::IntegratorPseudoTransientForwardSensitivity";
342 template<
class Scalar>
349 auto out = Teuchos::fancyOStream( in_out.
getOStream() );
350 out->setOutputToRootOnly(0);
351 *out << description() <<
"::describe" << std::endl;
352 state_integrator_->describe(in_out, verbLevel);
353 sens_integrator_->describe(in_out, verbLevel);
356 template<
class Scalar>
361 state_integrator_->setParameterList(inputPL);
362 sens_integrator_->setParameterList(inputPL);
364 inputPL->
sublist(
"Sensitivities").
get(
"Reuse State Linear Solver",
false);
366 inputPL->
sublist(
"Sensitivities").
get(
"Force W Update",
false);
369 template<
class Scalar>
374 state_integrator_->unsetParameterList();
375 return sens_integrator_->unsetParameterList();
378 template<
class Scalar>
386 state_integrator_->getValidParameters();
391 pl->
sublist(
"Sensitivities").
set(
"Reuse State Linear Solver",
false);
392 pl->
sublist(
"Sensitivities").
set(
"Force W Update",
false);
397 template<
class Scalar>
402 return state_integrator_->getNonconstParameterList();
405 template <
class Scalar>
415 if (inputPL != Teuchos::null) {
416 *pl = inputPL->
sublist(
"Sensitivities");
418 reuse_solver_ = pl->
get(
"Reuse State Linear Solver",
false);
419 force_W_update_ = pl->
get(
"Force W Update",
true);
420 pl->
remove(
"Reuse State Linear Solver");
421 pl->
remove(
"Force W Update");
425 template<
class Scalar>
432 using Teuchos::rcp_dynamic_cast;
437 using Thyra::createMembers;
438 using Thyra::multiVectorProductVector;
444 RCP<ParameterList> shPL =
445 Teuchos::sublist(state_integrator_->getIntegratorParameterList(),
446 "Solution History",
true);
447 solutionHistory_ = createSolutionHistoryPL<Scalar>(shPL);
449 const int num_param =
450 rcp_dynamic_cast<
const DMVPV>(sens_integrator_->getX())->getMultiVector()->
domain()->dim();
451 RCP<const VectorSpaceBase<Scalar> > x_space = model_->get_x_space();
452 RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> > prod_space =
453 Thyra::multiVectorProductVectorSpace(x_space, num_param+1);
457 RCP<const SolutionHistory<Scalar> > state_solution_history =
458 state_integrator_->getSolutionHistory();
459 int num_states = state_solution_history->getNumStates();
460 for (
int i=0; i<num_states; ++i) {
461 RCP<const SolutionState<Scalar> > state = (*state_solution_history)[i];
464 RCP< MultiVectorBase<Scalar> > x_mv =
465 createMembers(x_space, num_param+1);
466 assign(x_mv->col(0).ptr(), *(state->getX()));
467 assign(x_mv->subView(rng).ptr(), zero);
468 RCP<VectorBase<Scalar> > x = multiVectorProductVector(prod_space, x_mv);
471 RCP<VectorBase<Scalar> > x_dot;
472 if (state->getXDot() != Teuchos::null) {
473 RCP< MultiVectorBase<Scalar> > x_dot_mv =
474 createMembers(x_space, num_param+1);
475 assign(x_dot_mv->col(0).ptr(), *(state->getXDot()));
476 assign(x_dot_mv->subView(rng).ptr(), zero);
477 x_dot = multiVectorProductVector(prod_space, x_dot_mv);
481 RCP<VectorBase<Scalar> > x_dot_dot;
482 if (state->getXDotDot() != Teuchos::null) {
483 RCP< MultiVectorBase<Scalar> > x_dot_dot_mv =
484 createMembers(x_space, num_param+1);
485 assign(x_dot_dot_mv->col(0).ptr(), *(state->getXDotDot()));
486 assign(x_dot_dot_mv->subView(rng).ptr(), zero);
487 x_dot_dot = multiVectorProductVector(prod_space, x_dot_dot_mv);
490 RCP<SolutionState<Scalar> > prod_state = state->clone();
492 prod_state->setXDot(x_dot);
493 prod_state->setXDotDot(x_dot_dot);
494 solutionHistory_->addState(prod_state);
497 RCP<const VectorBase<Scalar> > frozen_x =
498 state_solution_history->getCurrentState()->getX();
499 RCP<const VectorBase<Scalar> > frozen_x_dot =
500 state_solution_history->getCurrentState()->getXDot();
501 RCP<const VectorBase<Scalar> > frozen_x_dot_dot =
502 state_solution_history->getCurrentState()->getXDotDot();
503 RCP<const SolutionHistory<Scalar> > sens_solution_history =
504 sens_integrator_->getSolutionHistory();
505 num_states = sens_solution_history->getNumStates();
506 for (
int i=0; i<num_states; ++i) {
507 RCP<const SolutionState<Scalar> > state = (*sens_solution_history)[i];
510 RCP< MultiVectorBase<Scalar> > x_mv =
511 createMembers(x_space, num_param+1);
512 RCP<const MultiVectorBase<Scalar> > dxdp =
513 rcp_dynamic_cast<
const DMVPV>(state->getX())->getMultiVector();
514 assign(x_mv->col(0).ptr(), *(frozen_x));
515 assign(x_mv->subView(rng).ptr(), *dxdp);
516 RCP<VectorBase<Scalar> > x = multiVectorProductVector(prod_space, x_mv);
519 RCP<VectorBase<Scalar> > x_dot;
520 if (state->getXDot() != Teuchos::null) {
521 RCP< MultiVectorBase<Scalar> > x_dot_mv =
522 createMembers(x_space, num_param+1);
523 RCP<const MultiVectorBase<Scalar> > dxdotdp =
524 rcp_dynamic_cast<
const DMVPV>(state->getXDot())->getMultiVector();
525 assign(x_dot_mv->col(0).ptr(), *(frozen_x_dot));
526 assign(x_dot_mv->subView(rng).ptr(), *dxdotdp);
527 x_dot = multiVectorProductVector(prod_space, x_dot_mv);
531 RCP<VectorBase<Scalar> > x_dot_dot;
532 if (state->getXDotDot() != Teuchos::null) {
533 RCP< MultiVectorBase<Scalar> > x_dot_dot_mv =
534 createMembers(x_space, num_param+1);
535 RCP<const MultiVectorBase<Scalar> > dxdotdotdp =
536 rcp_dynamic_cast<
const DMVPV>(state->getXDotDot())->getMultiVector();
537 assign(x_dot_dot_mv->col(0).ptr(), *(frozen_x_dot_dot));
538 assign(x_dot_dot_mv->subView(rng).ptr(), *dxdotdotdp);
539 x_dot_dot = multiVectorProductVector(prod_space, x_dot_dot_mv);
542 RCP<SolutionState<Scalar> > prod_state = state->clone();
544 prod_state->setXDot(x_dot);
545 prod_state->setXDotDot(x_dot_dot);
546 solutionHistory_->addState(prod_state,
false);
551 template<
class Scalar>
563 template<
class Scalar>
567 std::string stepperType)
575 template<
class Scalar>
585 #endif // Tempus_IntegratorPseudoTransientForwardSensitivity_impl_hpp
virtual Teuchos::RCP< Stepper< Scalar > > getStepper() const override
Get the Stepper.
virtual Status getStatus() const override
Get Status.
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &pl) override
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< Teuchos::ParameterList > getNonconstParameterList() override
Teuchos::RCP< SensitivityModelEvaluatorBase< Scalar > > wrapStaggeredFSAModelEvaluator(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, const Teuchos::RCP< const Teuchos::ParameterList > &pList=Teuchos::null)
virtual Teuchos::RCP< TimeStepControl< Scalar > > getNonConstTimeStepControl() override
static Teuchos::RCP< const Teuchos::ParameterList > getValidParameters()
Teuchos::RCP< Tempus::IntegratorPseudoTransientForwardSensitivity< Scalar > > integratorPseudoTransientForwardSensitivity(Teuchos::RCP< Teuchos::ParameterList > pList, const Teuchos::RCP< Thyra::ModelEvaluator< Scalar > > &model)
Nonmember constructor.
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)
Status
Status for the Integrator, the Stepper and the SolutionState.
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
Teuchos::RCP< IntegratorBasicOld< Scalar > > state_integrator_
virtual Teuchos::RCP< Teuchos::ParameterList > getTempusParameterList() override
Return a copy of the Tempus ParameterList.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const override
Teuchos::RCP< SensitivityModelEvaluatorBase< Scalar > > sens_model_
virtual Scalar getTime() const override
Get current time.
ParameterList & setParameters(const ParameterList &source)
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getDXDotDp() const
Teuchos::RCP< SensitivityModelEvaluatorBase< Scalar > > createSensitivityModel(const Teuchos::RCP< Thyra::ModelEvaluator< Scalar > > &model, const Teuchos::RCP< Teuchos::ParameterList > &inputPL)
RCP< std::basic_ostream< char_type, traits_type > > getOStream()
IntegratorPseudoTransientForwardSensitivity()
Destructor.
Time integrator suitable for pseudotransient forward sensitivity analysis.
virtual RCP< const VectorSpaceBase< Scalar > > domain() const
virtual void setTempusParameterList(Teuchos::RCP< Teuchos::ParameterList > pl) override
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
void buildSolutionHistory()
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList() override
Teuchos::RCP< IntegratorBasicOld< Scalar > > sens_integrator_
virtual Teuchos::RCP< const TimeStepControl< Scalar > > getTimeStepControl() const override
Get the TimeStepControl.
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getDXDotDotDp() const
Teuchos::RCP< Thyra::ModelEvaluator< Scalar > > model_
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getXDot() const
Get current the time derivative of the solution, xdot.
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getX() const
Get current the solution, x.