9#ifndef Tempus_IntegratorForwardSensitivity_impl_hpp
10#define Tempus_IntegratorForwardSensitivity_impl_hpp
13#include "Thyra_DefaultMultiVectorProductVector.hpp"
14#include "Thyra_VectorStdOps.hpp"
15#include "Thyra_MultiVectorStdOps.hpp"
17#include "Tempus_CombinedForwardSensitivityModelEvaluator.hpp"
22template <
class Scalar>
25 const Teuchos::RCP<IntegratorBasic<Scalar>> &integrator,
26 const Teuchos::RCP<SensitivityModelEvaluatorBase<Scalar>> &sens_model,
27 const Teuchos::RCP<StepperStaggeredForwardSensitivity<Scalar>> &sens_stepper,
28 bool use_combined_method)
30 , integrator_(integrator)
31 , sens_model_(sens_model)
32 , sens_stepper_(sens_stepper)
33 , use_combined_method_(use_combined_method)
42 integrator_ = createIntegratorBasic<Scalar>();
43 integrator_->initialize();
51 if (use_combined_method_)
52 integrator_->setModel(sens_model_);
54 integrator_->setStepper(sens_stepper_);
68 using Teuchos::rcp_dynamic_cast;
69 using Thyra::VectorSpaceBase;
71 using Thyra::createMember;
72 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
77 RCP< const VectorSpaceBase<Scalar> > space;
78 if (use_combined_method_)
79 space = sens_model_->get_x_space();
81 space = sens_stepper_->get_x_space();
82 RCP<DMVPV> X = rcp_dynamic_cast<DMVPV>(createMember(space));
83 RCP<DMVPV> Xdot = rcp_dynamic_cast<DMVPV>(createMember(space));
84 RCP<DMVPV> Xdotdot = rcp_dynamic_cast<DMVPV>(createMember(space));
86 const int num_param = X->getNonconstMultiVector()->domain()->dim()-1;
87 const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
88 const Teuchos::Range1D rng(1,num_param);
91 assign(X->getNonconstMultiVector()->col(0).ptr(), *x0);
92 if (DxDp0 == Teuchos::null)
93 assign(X->getNonconstMultiVector()->subView(rng).ptr(), zero);
95 assign(X->getNonconstMultiVector()->subView(rng).ptr(), *DxDp0);
98 if (xdot0 == Teuchos::null)
99 assign(Xdot->getNonconstMultiVector()->col(0).ptr(), zero);
101 assign(Xdot->getNonconstMultiVector()->col(0).ptr(), *xdot0);
102 if (DxdotDp0 == Teuchos::null)
103 assign(Xdot->getNonconstMultiVector()->subView(rng).ptr(), zero);
105 assign(Xdot->getNonconstMultiVector()->subView(rng).ptr(), *DxdotDp0);
108 if (xdotdot0 == Teuchos::null)
109 assign(Xdotdot->getNonconstMultiVector()->col(0).ptr(), zero);
111 assign(Xdotdot->getNonconstMultiVector()->col(0).ptr(), *xdotdot0);
112 if (DxdotDp0 == Teuchos::null)
113 assign(Xdotdot->getNonconstMultiVector()->subView(rng).ptr(), zero);
115 assign(Xdotdot->getNonconstMultiVector()->subView(rng).ptr(), *DxdotdotDp0);
117 integrator_->initializeSolutionHistory(t0, X, Xdot, Xdotdot);
120template<
class Scalar>
121Teuchos::RCP<const Thyra::VectorBase<Scalar> >
126 using Teuchos::rcp_dynamic_cast;
127 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
129 RCP<const DMVPV> X = rcp_dynamic_cast<const DMVPV>(integrator_->getX());
130 return X->getMultiVector()->col(0);
133template<
class Scalar>
134Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
139 using Teuchos::rcp_dynamic_cast;
140 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
142 RCP<const DMVPV> X = rcp_dynamic_cast<const DMVPV>(integrator_->getX());
143 const int num_param = X->getMultiVector()->domain()->dim()-1;
144 const Teuchos::Range1D rng(1,num_param);
145 return X->getMultiVector()->subView(rng);
148template<
class Scalar>
149Teuchos::RCP<const Thyra::VectorBase<Scalar> >
154 using Teuchos::rcp_dynamic_cast;
155 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
157 RCP<const DMVPV> Xdot = rcp_dynamic_cast<const DMVPV>(integrator_->getXDot());
158 return Xdot->getMultiVector()->col(0);
161template<
class Scalar>
162Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
167 using Teuchos::rcp_dynamic_cast;
168 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
170 RCP<const DMVPV> Xdot = rcp_dynamic_cast<const DMVPV>(integrator_->getXDot());
171 const int num_param = Xdot->getMultiVector()->domain()->dim()-1;
172 const Teuchos::Range1D rng(1,num_param);
173 return Xdot->getMultiVector()->subView(rng);
176template<
class Scalar>
177Teuchos::RCP<const Thyra::VectorBase<Scalar> >
182 using Teuchos::rcp_dynamic_cast;
183 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
185 RCP<const DMVPV> Xdotdot =
186 rcp_dynamic_cast<const DMVPV>(integrator_->getXDotDot());
187 return Xdotdot->getMultiVector()->col(0);
190template<
class Scalar>
191Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
196 using Teuchos::rcp_dynamic_cast;
197 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
199 RCP<const DMVPV> Xdotdot =
200 rcp_dynamic_cast<const DMVPV>(integrator_->getXDotDot());
201 const int num_param = Xdotdot->getMultiVector()->domain()->dim()-1;
202 const Teuchos::Range1D rng(1,num_param);
203 return Xdotdot->getMultiVector()->subView(rng);
206template<
class Scalar>
207Teuchos::RCP<const Thyra::VectorBase<Scalar> >
211 typedef Thyra::ModelEvaluatorBase MEB;
215 Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > smodel;
216 if (use_combined_method_)
217 smodel = sens_model_;
219 smodel = sens_stepper_->getModel();
220 MEB::InArgs<Scalar> inargs = smodel->getNominalValues();
221 MEB::OutArgs<Scalar> outargs = smodel->createOutArgs();
222 inargs.set_t(integrator_->getTime());
223 inargs.set_x(integrator_->getX());
224 if (inargs.supports(MEB::IN_ARG_x_dot))
225 inargs.set_x_dot(integrator_->getXDot());
226 if (inargs.supports(MEB::IN_ARG_x_dot_dot))
227 inargs.set_x_dot_dot(integrator_->getXDotDot());
229 Teuchos::RCP<Thyra::VectorBase<Scalar> > g =
230 Thyra::createMember(smodel->get_g_space(1));
233 smodel->evalModel(inargs, outargs);
237template<
class Scalar>
238Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
242 typedef Thyra::ModelEvaluatorBase MEB;
243 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
247 Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > smodel;
248 if (use_combined_method_)
249 smodel = sens_model_;
251 smodel = sens_stepper_->getModel();
252 MEB::InArgs<Scalar> inargs = smodel->getNominalValues();
253 MEB::OutArgs<Scalar> outargs = smodel->createOutArgs();
254 inargs.set_t(integrator_->getTime());
255 inargs.set_x(integrator_->getX());
256 if (inargs.supports(MEB::IN_ARG_x_dot))
257 inargs.set_x_dot(integrator_->getXDot());
258 if (inargs.supports(MEB::IN_ARG_x_dot_dot))
259 inargs.set_x_dot_dot(integrator_->getXDotDot());
261 Teuchos::RCP<Thyra::VectorBase<Scalar> > G =
262 Thyra::createMember(smodel->get_g_space(0));
263 Teuchos::RCP<DMVPV> dgdp = Teuchos::rcp_dynamic_cast<DMVPV>(G);
266 smodel->evalModel(inargs, outargs);
267 return dgdp->getMultiVector();
270template<
class Scalar>
275 std::string name =
"Tempus::IntegratorForwardSensitivity";
279template<
class Scalar>
283 Teuchos::FancyOStream &out,
284 const Teuchos::EVerbosityLevel verbLevel)
const
286 auto l_out = Teuchos::fancyOStream( out.getOStream() );
287 Teuchos::OSTab ostab(*l_out, 2, this->description());
288 l_out->setOutputToRootOnly(0);
290 *l_out << description() <<
"::describe" << std::endl;
291 integrator_->describe(*l_out, verbLevel);
294template <
class Scalar>
299 if (use_combined_method_)
301 return sens_stepper_->getStepMode();
305template<
class Scalar>
306Teuchos::RCP<IntegratorForwardSensitivity<Scalar> >
308 Teuchos::RCP<Teuchos::ParameterList> pList,
315 Teuchos::RCP<SensitivityModelEvaluatorBase<Scalar> > sens_model;
316 Teuchos::RCP<StepperStaggeredForwardSensitivity<Scalar>> sens_stepper;
319 auto fwd_integrator = createIntegratorBasic<Scalar>(pList, model);
322 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::rcp(
new Teuchos::ParameterList);
324 Teuchos::RCP<const Teuchos::ParameterList> integrator_pl = fwd_integrator->getValidParameters();
325 Teuchos::RCP<const Teuchos::ParameterList> sensitivity_pl =
327 pl->setParameters(*integrator_pl);
328 Teuchos::ParameterList &spl = pl->sublist(
"Sensitivities");
329 spl.setParameters(*sensitivity_pl);
330 spl.set(
"Sensitivity Method",
"Combined");
331 spl.set(
"Reuse State Linear Solver",
false);
333 pList->setParametersNotAlreadySet(*pl);
335 auto sens_pl = Teuchos::sublist(pList,
"Sensitivities",
false);
336 std::string integratorName = pList->get<std::string>(
"Integrator Name",
"Default Integrator");
337 std::string stepperName = pList->sublist(integratorName).get<std::string>(
"Stepper Name");
338 auto stepper_pl = Teuchos::sublist(pList, stepperName,
true);
339 std::string sensitivity_method = sens_pl->get<std::string>(
"Sensitivity Method");
340 bool use_combined_method = sensitivity_method ==
"Combined";
345 sens_pl->remove(
"Sensitivity Method");
347 if (use_combined_method)
349 sens_pl->remove(
"Reuse State Linear Solver");
351 model, sens_residual_model, sens_solve_model, sens_pl);
352 fwd_integrator = createIntegratorBasic<Scalar>(pList, sens_model);
356 sens_stepper = Teuchos::rcp(
new StepperStaggeredForwardSensitivity<Scalar>(model, sens_residual_model, sens_solve_model, stepper_pl, sens_pl));
357 auto fsa_staggered_me = Teuchos::rcp_const_cast<Thyra::ModelEvaluator<Scalar>>(sens_stepper->getModel());
358 fwd_integrator = createIntegratorBasic<Scalar>(pList, fsa_staggered_me);
359 fwd_integrator->setStepper(sens_stepper);
364 Teuchos::RCP<IntegratorForwardSensitivity<Scalar>> integrator =
365 Teuchos::rcp(
new IntegratorForwardSensitivity<Scalar>(model, fwd_integrator, sens_model, sens_stepper, use_combined_method));
370template<
class Scalar>
371Teuchos::RCP<IntegratorForwardSensitivity<Scalar> >
375 Teuchos::RCP<IntegratorForwardSensitivity<Scalar> > integrator =
376 Teuchos::rcp(
new IntegratorForwardSensitivity<Scalar>());
static Teuchos::RCP< const Teuchos::ParameterList > getValidParameters()
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getDXDotDp() const
virtual void setStepper(Teuchos::RCP< Thyra::ModelEvaluator< Scalar > > model)
Set the Stepper.
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getXDotDot() const
Get current the second time derivative of the solution, xdotdot, only. This is the first block only a...
SensitivityStepMode getStepMode() const
What mode the current time integration step is in.
virtual void initializeSolutionHistory(Teuchos::RCP< SolutionState< Scalar > > state=Teuchos::null)
Set the initial state which has the initial conditions.
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getDxDp() const
Get the forward sensitivities .
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getX() const
Get the current solution, x, only. If looking for the solution vector and the sensitivities,...
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getDXDotDotDp() const
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getG() const
Return response function g.
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getDgDp() const
Return forward sensitivity stored in Jacobian format.
Teuchos::RCP< IntegratorBasic< Scalar > > integrator_
std::string description() const override
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const override
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getXDot() const
Get current the time derivative of the solution, xdot, only. This is the first block only and not the...
IntegratorForwardSensitivity()
Destructor.
Teuchos::RCP< SensitivityModelEvaluatorBase< Scalar > > wrapCombinedFSAModelEvaluator(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 Teuchos::RCP< const Teuchos::ParameterList > &pList=Teuchos::null)
Teuchos::RCP< IntegratorForwardSensitivity< Scalar > > createIntegratorForwardSensitivity()
Nonmember constructor.