9 #ifndef Tempus_IntegratorPseudoTransientForwardSensitivity_impl_hpp
10 #define Tempus_IntegratorPseudoTransientForwardSensitivity_impl_hpp
13 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
14 #include "Thyra_DefaultMultiVectorProductVector.hpp"
15 #include "Thyra_VectorStdOps.hpp"
16 #include "Thyra_MultiVectorStdOps.hpp"
20 template<
class Scalar>
23 Teuchos::RCP<Teuchos::ParameterList> inputPL,
24 const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model) :
33 template<
class Scalar>
36 const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model,
37 std::string stepperType) :
39 force_W_update_(false)
47 template<
class Scalar>
51 force_W_update_(false)
57 template<
class Scalar>
63 using Thyra::VectorBase;
66 bool state_status = state_integrator_->advanceTime();
69 sens_model_->setForwardSolutionState(state_integrator_->getCurrentState());
73 state_integrator_->getStepper()->getSolver() != Teuchos::null) {
74 sens_model_->setSolver(state_integrator_->getStepper()->getSolver(),
79 bool sens_status = sens_integrator_->advanceTime();
81 buildSolutionHistory();
83 return state_status && sens_status;
86 template<
class Scalar>
92 using Thyra::VectorBase;
95 bool state_status = state_integrator_->advanceTime(timeFinal);
98 sens_model_->setForwardSolutionState(state_integrator_->getCurrentState());
102 state_integrator_->getStepper()->getSolver() != Teuchos::null) {
103 sens_model_->setSolver(state_integrator_->getStepper()->getSolver(),
108 bool sens_status = sens_integrator_->advanceTime(timeFinal);
110 buildSolutionHistory();
112 return state_status && sens_status;
115 template<
class Scalar>
120 return solutionHistory_->getCurrentTime();
123 template<
class Scalar>
128 return solutionHistory_->getCurrentIndex();
131 template<
class Scalar>
136 Status state_status = state_integrator_->getStatus();
137 Status sens_status = sens_integrator_->getStatus();
145 template<
class Scalar>
146 Teuchos::RCP<Stepper<Scalar> >
150 return state_integrator_->getStepper();
153 template<
class Scalar>
154 Teuchos::RCP<Teuchos::ParameterList>
158 return state_integrator_->getTempusParameterList();
161 template<
class Scalar>
166 state_integrator_->setTempusParameterList(pl);
167 sens_integrator_->setTempusParameterList(pl);
170 template<
class Scalar>
171 Teuchos::RCP<const SolutionHistory<Scalar> >
175 return solutionHistory_;
178 template<
class Scalar>
179 Teuchos::RCP<const TimeStepControl<Scalar> >
183 return state_integrator_->getTimeStepControl();
186 template<
class Scalar>
187 Teuchos::RCP<TimeStepControl<Scalar> >
191 return state_integrator_->getNonConstTimeStepControl();
194 template<
class Scalar>
197 Teuchos::RCP<
const Thyra::VectorBase<Scalar> > x0,
198 Teuchos::RCP<
const Thyra::VectorBase<Scalar> > xdot0,
199 Teuchos::RCP<
const Thyra::VectorBase<Scalar> > xdotdot0,
200 Teuchos::RCP<
const Thyra::MultiVectorBase<Scalar> > DxDp0,
201 Teuchos::RCP<
const Thyra::MultiVectorBase<Scalar> > DxdotDp0,
202 Teuchos::RCP<
const Thyra::MultiVectorBase<Scalar> > DxdotdotDp0)
205 using Teuchos::rcp_dynamic_cast;
206 using Thyra::VectorSpaceBase;
208 using Thyra::createMember;
209 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
214 RCP< const VectorSpaceBase<Scalar> > space = sens_model_->get_x_space();
215 RCP<DMVPV> X = rcp_dynamic_cast<DMVPV>(createMember(space));
216 RCP<DMVPV> Xdot = rcp_dynamic_cast<DMVPV>(createMember(space));
217 RCP<DMVPV> Xdotdot = rcp_dynamic_cast<DMVPV>(createMember(space));
218 const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
221 if (DxDp0 == Teuchos::null)
222 assign(X->getNonconstMultiVector().ptr(), zero);
224 assign(X->getNonconstMultiVector().ptr(), *DxDp0);
227 if (DxdotDp0 == Teuchos::null)
228 assign(Xdot->getNonconstMultiVector().ptr(), zero);
230 assign(Xdot->getNonconstMultiVector().ptr(), *DxdotDp0);
233 if (DxdotDp0 == Teuchos::null)
234 assign(Xdotdot->getNonconstMultiVector().ptr(), zero);
236 assign(Xdotdot->getNonconstMultiVector().ptr(), *DxdotdotDp0);
238 state_integrator_->initializeSolutionHistory(t0, x0, xdot0, xdotdot0);
239 sens_integrator_->initializeSolutionHistory(t0, X, Xdot, Xdotdot);
242 template<
class Scalar>
243 Teuchos::RCP<const Thyra::VectorBase<Scalar> >
248 using Teuchos::rcp_dynamic_cast;
249 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
252 rcp_dynamic_cast<
const DMVPV>(solutionHistory_->getCurrentState()->getX());
253 return X->getMultiVector()->col(0);
256 template<
class Scalar>
257 Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
262 using Teuchos::rcp_dynamic_cast;
263 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
266 rcp_dynamic_cast<
const DMVPV>(solutionHistory_->getCurrentState()->getX());
267 const int num_param = X->getMultiVector()->domain()->dim()-1;
268 const Teuchos::Range1D rng(1,num_param);
269 return X->getMultiVector()->subView(rng);
272 template<
class Scalar>
273 Teuchos::RCP<const Thyra::VectorBase<Scalar> >
278 using Teuchos::rcp_dynamic_cast;
279 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
281 RCP<const DMVPV> Xdot =
282 rcp_dynamic_cast<
const DMVPV>(solutionHistory_->getCurrentState()->getXDot());
283 return Xdot->getMultiVector()->col(0);
286 template<
class Scalar>
287 Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
292 using Teuchos::rcp_dynamic_cast;
293 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
295 RCP<const DMVPV> Xdot =
296 rcp_dynamic_cast<
const DMVPV>(solutionHistory_->getCurrentState()->getXDot());
297 const int num_param = Xdot->getMultiVector()->domain()->dim()-1;
298 const Teuchos::Range1D rng(1,num_param);
299 return Xdot->getMultiVector()->subView(rng);
302 template<
class Scalar>
303 Teuchos::RCP<const Thyra::VectorBase<Scalar> >
308 using Teuchos::rcp_dynamic_cast;
309 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
311 RCP<const DMVPV> Xdotdot =
312 rcp_dynamic_cast<
const DMVPV>(solutionHistory_->getCurrentState()->getXDotDot());
313 return Xdotdot->getMultiVector()->col(0);
316 template<
class Scalar>
317 Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
322 using Teuchos::rcp_dynamic_cast;
323 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
325 RCP<const DMVPV> Xdotdot =
326 rcp_dynamic_cast<
const DMVPV>(solutionHistory_->getCurrentState()->getXDotDot());
327 const int num_param = Xdotdot->getMultiVector()->domain()->dim()-1;
328 const Teuchos::Range1D rng(1,num_param);
329 return Xdotdot->getMultiVector()->subView(rng);
332 template<
class Scalar>
337 std::string name =
"Tempus::IntegratorPseudoTransientForwardSensitivity";
341 template<
class Scalar>
345 Teuchos::FancyOStream &out,
346 const Teuchos::EVerbosityLevel verbLevel)
const
348 out << description() <<
"::describe" << std::endl;
349 state_integrator_->describe(out, verbLevel);
350 sens_integrator_->describe(out, verbLevel);
353 template<
class Scalar>
358 state_integrator_->setParameterList(inputPL);
359 sens_integrator_->setParameterList(inputPL);
361 inputPL->sublist(
"Sensitivities").get(
"Reuse State Linear Solver",
false);
363 inputPL->sublist(
"Sensitivities").get(
"Force W Update",
false);
366 template<
class Scalar>
367 Teuchos::RCP<Teuchos::ParameterList>
371 state_integrator_->unsetParameterList();
372 return sens_integrator_->unsetParameterList();
375 template<
class Scalar>
376 Teuchos::RCP<const Teuchos::ParameterList>
380 Teuchos::RCP<Teuchos::ParameterList> pl =
381 Teuchos::rcp(
new Teuchos::ParameterList);
382 Teuchos::RCP<const Teuchos::ParameterList> integrator_pl =
383 state_integrator_->getValidParameters();
384 Teuchos::RCP<const Teuchos::ParameterList> sensitivity_pl =
386 pl->setParameters(*integrator_pl);
387 pl->sublist(
"Sensitivities").setParameters(*sensitivity_pl);
388 pl->sublist(
"Sensitivities").set(
"Reuse State Linear Solver",
false);
389 pl->sublist(
"Sensitivities").set(
"Force W Update",
false);
394 template<
class Scalar>
395 Teuchos::RCP<Teuchos::ParameterList>
399 return state_integrator_->getNonconstParameterList();
402 template <
class Scalar>
403 Teuchos::RCP<SensitivityModelEvaluatorBase<Scalar> >
406 const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model,
407 const Teuchos::RCP<Teuchos::ParameterList>& inputPL)
411 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
412 if (inputPL != Teuchos::null) {
413 *pl = inputPL->sublist(
"Sensitivities");
415 reuse_solver_ = pl->get(
"Reuse State Linear Solver",
false);
416 force_W_update_ = pl->get(
"Force W Update",
true);
417 pl->remove(
"Reuse State Linear Solver");
418 pl->remove(
"Force W Update");
422 template<
class Scalar>
429 using Teuchos::rcp_dynamic_cast;
430 using Teuchos::ParameterList;
431 using Thyra::VectorBase;
432 using Thyra::MultiVectorBase;
433 using Thyra::VectorSpaceBase;
434 using Thyra::createMembers;
435 using Thyra::multiVectorProductVector;
437 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
441 RCP<ParameterList> shPL =
442 Teuchos::sublist(state_integrator_->getIntegratorParameterList(),
443 "Solution History",
true);
446 const int num_param =
447 rcp_dynamic_cast<
const DMVPV>(sens_integrator_->getX())->getMultiVector()->domain()->dim();
448 RCP<const VectorSpaceBase<Scalar> > x_space = model_->get_x_space();
449 RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> > prod_space =
450 Thyra::multiVectorProductVectorSpace(x_space, num_param+1);
451 const Teuchos::Range1D rng(1,num_param);
452 const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
454 RCP<const SolutionHistory<Scalar> > state_solution_history =
455 state_integrator_->getSolutionHistory();
456 int num_states = state_solution_history->getNumStates();
457 for (
int i=0; i<num_states; ++i) {
458 RCP<const SolutionState<Scalar> > state = (*state_solution_history)[i];
461 RCP< MultiVectorBase<Scalar> > x_mv =
462 createMembers(x_space, num_param+1);
463 assign(x_mv->col(0).ptr(), *(state->getX()));
464 assign(x_mv->subView(rng).ptr(), zero);
465 RCP<VectorBase<Scalar> > x = multiVectorProductVector(prod_space, x_mv);
468 RCP<VectorBase<Scalar> > x_dot;
469 if (state->getXDot() != Teuchos::null) {
470 RCP< MultiVectorBase<Scalar> > x_dot_mv =
471 createMembers(x_space, num_param+1);
472 assign(x_dot_mv->col(0).ptr(), *(state->getXDot()));
473 assign(x_dot_mv->subView(rng).ptr(), zero);
474 x_dot = multiVectorProductVector(prod_space, x_dot_mv);
478 RCP<VectorBase<Scalar> > x_dot_dot;
479 if (state->getXDotDot() != Teuchos::null) {
480 RCP< MultiVectorBase<Scalar> > x_dot_dot_mv =
481 createMembers(x_space, num_param+1);
482 assign(x_dot_dot_mv->col(0).ptr(), *(state->getXDotDot()));
483 assign(x_dot_dot_mv->subView(rng).ptr(), zero);
484 x_dot_dot = multiVectorProductVector(prod_space, x_dot_dot_mv);
487 RCP<SolutionState<Scalar> > prod_state = state->clone();
489 prod_state->setXDot(x_dot);
490 prod_state->setXDotDot(x_dot_dot);
491 solutionHistory_->addState(prod_state);
494 RCP<const VectorBase<Scalar> > frozen_x =
495 state_solution_history->getCurrentState()->getX();
496 RCP<const VectorBase<Scalar> > frozen_x_dot =
497 state_solution_history->getCurrentState()->getXDot();
498 RCP<const VectorBase<Scalar> > frozen_x_dot_dot =
499 state_solution_history->getCurrentState()->getXDotDot();
500 RCP<const SolutionHistory<Scalar> > sens_solution_history =
501 sens_integrator_->getSolutionHistory();
502 num_states = sens_solution_history->getNumStates();
503 for (
int i=0; i<num_states; ++i) {
504 RCP<const SolutionState<Scalar> > state = (*sens_solution_history)[i];
507 RCP< MultiVectorBase<Scalar> > x_mv =
508 createMembers(x_space, num_param+1);
509 RCP<const MultiVectorBase<Scalar> > dxdp =
510 rcp_dynamic_cast<
const DMVPV>(state->getX())->getMultiVector();
511 assign(x_mv->col(0).ptr(), *(frozen_x));
512 assign(x_mv->subView(rng).ptr(), *dxdp);
513 RCP<VectorBase<Scalar> > x = multiVectorProductVector(prod_space, x_mv);
516 RCP<VectorBase<Scalar> > x_dot;
517 if (state->getXDot() != Teuchos::null) {
518 RCP< MultiVectorBase<Scalar> > x_dot_mv =
519 createMembers(x_space, num_param+1);
520 RCP<const MultiVectorBase<Scalar> > dxdotdp =
521 rcp_dynamic_cast<
const DMVPV>(state->getXDot())->getMultiVector();
522 assign(x_dot_mv->col(0).ptr(), *(frozen_x_dot));
523 assign(x_dot_mv->subView(rng).ptr(), *dxdotdp);
524 x_dot = multiVectorProductVector(prod_space, x_dot_mv);
528 RCP<VectorBase<Scalar> > x_dot_dot;
529 if (state->getXDotDot() != Teuchos::null) {
530 RCP< MultiVectorBase<Scalar> > x_dot_dot_mv =
531 createMembers(x_space, num_param+1);
532 RCP<const MultiVectorBase<Scalar> > dxdotdotdp =
533 rcp_dynamic_cast<
const DMVPV>(state->getXDotDot())->getMultiVector();
534 assign(x_dot_dot_mv->col(0).ptr(), *(frozen_x_dot_dot));
535 assign(x_dot_dot_mv->subView(rng).ptr(), *dxdotdotdp);
536 x_dot_dot = multiVectorProductVector(prod_space, x_dot_dot_mv);
539 RCP<SolutionState<Scalar> > prod_state = state->clone();
541 prod_state->setXDot(x_dot);
542 prod_state->setXDotDot(x_dot_dot);
543 solutionHistory_->addState(prod_state);
548 template<
class Scalar>
549 Teuchos::RCP<Tempus::IntegratorPseudoTransientForwardSensitivity<Scalar> >
551 Teuchos::RCP<Teuchos::ParameterList> pList,
552 const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model)
554 Teuchos::RCP<Tempus::IntegratorPseudoTransientForwardSensitivity<Scalar> > integrator =
560 template<
class Scalar>
561 Teuchos::RCP<Tempus::IntegratorPseudoTransientForwardSensitivity<Scalar> >
563 const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model,
564 std::string stepperType)
566 Teuchos::RCP<Tempus::IntegratorPseudoTransientForwardSensitivity<Scalar> > integrator =
572 template<
class Scalar>
573 Teuchos::RCP<Tempus::IntegratorPseudoTransientForwardSensitivity<Scalar> >
576 Teuchos::RCP<Tempus::IntegratorPseudoTransientForwardSensitivity<Scalar> > integrator =
582 #endif // Tempus_IntegratorPseudoTransientForwardSensitivity_impl_hpp
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getDxdotDp() const
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
virtual Teuchos::RCP< const SolutionHistory< Scalar > > getSolutionHistory() const override
Get the SolutionHistory.
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getDxDp() const
virtual int getIndex() const override
Get current index.
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getDxdotdotDp() const
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)
Non-member constructor.
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.
Status
Status for the Integrator, the Stepper and the SolutionState.
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const override
std::string description() const override
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_
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getXdot() const
Get current the time derivative of the solution, xdot.
virtual Scalar getTime() const override
Get current time.
Teuchos::RCP< SensitivityModelEvaluatorBase< Scalar > > createSensitivityModel(const Teuchos::RCP< Thyra::ModelEvaluator< Scalar > > &model, const Teuchos::RCP< Teuchos::ParameterList > &inputPL)
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getXdotdot() const
Get current the second time derivative of the solution, xdotdot.
IntegratorPseudoTransientForwardSensitivity()
Destructor.
Time integrator suitable for pseudotransient forward sensitivity analysis.
virtual void setTempusParameterList(Teuchos::RCP< Teuchos::ParameterList > pl) override
void buildSolutionHistory()
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList() override
virtual Teuchos::RCP< const TimeStepControl< Scalar > > getTimeStepControl() const override
Get the TimeStepControl.
Teuchos::RCP< Thyra::ModelEvaluator< Scalar > > model_
Teuchos::RCP< IntegratorBasic< Scalar > > sens_integrator_
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getX() const
Get current the solution, x.