Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_IntegratorForwardSensitivity_impl.hpp
Go to the documentation of this file.
1// @HEADER
2// ****************************************************************************
3// Tempus: Copyright (2017) Sandia Corporation
4//
5// Distributed under BSD 3-clause license (See accompanying file Copyright.txt)
6// ****************************************************************************
7// @HEADER
8
9#ifndef Tempus_IntegratorForwardSensitivity_impl_hpp
10#define Tempus_IntegratorForwardSensitivity_impl_hpp
11
13#include "Thyra_DefaultMultiVectorProductVector.hpp"
14#include "Thyra_VectorStdOps.hpp"
15#include "Thyra_MultiVectorStdOps.hpp"
16
17#include "Tempus_CombinedForwardSensitivityModelEvaluator.hpp"
19
20namespace Tempus {
21
22template <class Scalar>
24 const Teuchos::RCP<Thyra::ModelEvaluator<Scalar>> &model,
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)
29 : model_(model)
30 , integrator_(integrator)
31 , sens_model_(sens_model)
32 , sens_stepper_(sens_stepper)
33 , use_combined_method_(use_combined_method)
34{
35 integrator_->initialize();
36}
37
38template<class Scalar>
41{
42 integrator_ = createIntegratorBasic<Scalar>();
43 integrator_->initialize();
44}
45
46template<class Scalar>
49 Teuchos::RCP<Thyra::ModelEvaluator<Scalar> > model)
50{
51 if (use_combined_method_)
52 integrator_->setModel(sens_model_);
53 else
54 integrator_->setStepper(sens_stepper_);
55}
56
57template<class Scalar>
60 Teuchos::RCP<const Thyra::VectorBase<Scalar> > x0,
61 Teuchos::RCP<const Thyra::VectorBase<Scalar> > xdot0,
62 Teuchos::RCP<const Thyra::VectorBase<Scalar> > xdotdot0,
63 Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > DxDp0,
64 Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > DxdotDp0,
65 Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > DxdotdotDp0)
66{
67 using Teuchos::RCP;
68 using Teuchos::rcp_dynamic_cast;
69 using Thyra::VectorSpaceBase;
70 using Thyra::assign;
71 using Thyra::createMember;
72 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
73
74 //
75 // Create and initialize product X, Xdot, Xdotdot
76
77 RCP< const VectorSpaceBase<Scalar> > space;
78 if (use_combined_method_)
79 space = sens_model_->get_x_space();
80 else
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));
85
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);
89
90 // x
91 assign(X->getNonconstMultiVector()->col(0).ptr(), *x0);
92 if (DxDp0 == Teuchos::null)
93 assign(X->getNonconstMultiVector()->subView(rng).ptr(), zero);
94 else
95 assign(X->getNonconstMultiVector()->subView(rng).ptr(), *DxDp0);
96
97 // xdot
98 if (xdot0 == Teuchos::null)
99 assign(Xdot->getNonconstMultiVector()->col(0).ptr(), zero);
100 else
101 assign(Xdot->getNonconstMultiVector()->col(0).ptr(), *xdot0);
102 if (DxdotDp0 == Teuchos::null)
103 assign(Xdot->getNonconstMultiVector()->subView(rng).ptr(), zero);
104 else
105 assign(Xdot->getNonconstMultiVector()->subView(rng).ptr(), *DxdotDp0);
106
107 // xdotdot
108 if (xdotdot0 == Teuchos::null)
109 assign(Xdotdot->getNonconstMultiVector()->col(0).ptr(), zero);
110 else
111 assign(Xdotdot->getNonconstMultiVector()->col(0).ptr(), *xdotdot0);
112 if (DxdotDp0 == Teuchos::null)
113 assign(Xdotdot->getNonconstMultiVector()->subView(rng).ptr(), zero);
114 else
115 assign(Xdotdot->getNonconstMultiVector()->subView(rng).ptr(), *DxdotdotDp0);
116
117 integrator_->initializeSolutionHistory(t0, X, Xdot, Xdotdot);
118}
119
120template<class Scalar>
121Teuchos::RCP<const Thyra::VectorBase<Scalar> >
123getX() const
124{
125 using Teuchos::RCP;
126 using Teuchos::rcp_dynamic_cast;
127 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
128
129 RCP<const DMVPV> X = rcp_dynamic_cast<const DMVPV>(integrator_->getX());
130 return X->getMultiVector()->col(0);
131}
132
133template<class Scalar>
134Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
136getDxDp() const
137{
138 using Teuchos::RCP;
139 using Teuchos::rcp_dynamic_cast;
140 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
141
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);
146}
147
148template<class Scalar>
149Teuchos::RCP<const Thyra::VectorBase<Scalar> >
151getXDot() const
152{
153 using Teuchos::RCP;
154 using Teuchos::rcp_dynamic_cast;
155 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
156
157 RCP<const DMVPV> Xdot = rcp_dynamic_cast<const DMVPV>(integrator_->getXDot());
158 return Xdot->getMultiVector()->col(0);
159}
160
161template<class Scalar>
162Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
164getDXDotDp() const
165{
166 using Teuchos::RCP;
167 using Teuchos::rcp_dynamic_cast;
168 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
169
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);
174}
175
176template<class Scalar>
177Teuchos::RCP<const Thyra::VectorBase<Scalar> >
179getXDotDot() const
180{
181 using Teuchos::RCP;
182 using Teuchos::rcp_dynamic_cast;
183 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
184
185 RCP<const DMVPV> Xdotdot =
186 rcp_dynamic_cast<const DMVPV>(integrator_->getXDotDot());
187 return Xdotdot->getMultiVector()->col(0);
188}
189
190template<class Scalar>
191Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
193getDXDotDotDp() const
194{
195 using Teuchos::RCP;
196 using Teuchos::rcp_dynamic_cast;
197 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
198
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);
204}
205
206template<class Scalar>
207Teuchos::RCP<const Thyra::VectorBase<Scalar> >
209getG() const
210{
211 typedef Thyra::ModelEvaluatorBase MEB;
212
213 // Compute g which is computed by response 1 of the
214 // sensitivity model evaluator
215 Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > smodel;
216 if (use_combined_method_)
217 smodel = sens_model_;
218 else
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());
228
229 Teuchos::RCP<Thyra::VectorBase<Scalar> > g =
230 Thyra::createMember(smodel->get_g_space(1));
231 outargs.set_g(1, g);
232
233 smodel->evalModel(inargs, outargs);
234 return g;
235}
236
237template<class Scalar>
238Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
240getDgDp() const
241{
242 typedef Thyra::ModelEvaluatorBase MEB;
243 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
244
245 // Compute final dg/dp which is computed by response 0 of the
246 // sensitivity model evaluator
247 Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > smodel;
248 if (use_combined_method_)
249 smodel = sens_model_;
250 else
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());
260
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);
264 outargs.set_g(0, G);
265
266 smodel->evalModel(inargs, outargs);
267 return dgdp->getMultiVector();
268}
269
270template<class Scalar>
271std::string
273description() const
274{
275 std::string name = "Tempus::IntegratorForwardSensitivity";
276 return(name);
277}
278
279template<class Scalar>
280void
283 Teuchos::FancyOStream &out,
284 const Teuchos::EVerbosityLevel verbLevel) const
285{
286 auto l_out = Teuchos::fancyOStream( out.getOStream() );
287 Teuchos::OSTab ostab(*l_out, 2, this->description());
288 l_out->setOutputToRootOnly(0);
289
290 *l_out << description() << "::describe" << std::endl;
291 integrator_->describe(*l_out, verbLevel);
292}
293
294template <class Scalar>
297getStepMode() const
298{
299 if (use_combined_method_)
301 return sens_stepper_->getStepMode(); // Staggered case
302}
303
305template<class Scalar>
306Teuchos::RCP<IntegratorForwardSensitivity<Scalar> >
308 Teuchos::RCP<Teuchos::ParameterList> pList,
309 const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model,
310 const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& sens_residual_model,
311 const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& sens_solve_model)
312{
313
314
315 Teuchos::RCP<SensitivityModelEvaluatorBase<Scalar> > sens_model;
316 Teuchos::RCP<StepperStaggeredForwardSensitivity<Scalar>> sens_stepper;
317
318 //1. create integrator
319 auto fwd_integrator = createIntegratorBasic<Scalar>(pList, model);
320
321 //2. set parameter list
322 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::rcp(new Teuchos::ParameterList);
323 {
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);
332 }
333 pList->setParametersNotAlreadySet(*pl);
334
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";
341
342 //3. create sensitivity model and stepper
343 // createSensitivityModelAndStepper
344 {
345 sens_pl->remove("Sensitivity Method");
346
347 if (use_combined_method)
348 {
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);
353 }
354 else
355 {
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);
360 }
361 }
362
363 //4. initialize propoer integrator
364 Teuchos::RCP<IntegratorForwardSensitivity<Scalar>> integrator =
365 Teuchos::rcp(new IntegratorForwardSensitivity<Scalar>(model, fwd_integrator, sens_model, sens_stepper, use_combined_method));
366 return (integrator);
367}
368
370template<class Scalar>
371Teuchos::RCP<IntegratorForwardSensitivity<Scalar> >
373{
374
375 Teuchos::RCP<IntegratorForwardSensitivity<Scalar> > integrator =
376 Teuchos::rcp(new IntegratorForwardSensitivity<Scalar>());
377 return(integrator);
378}
379
380} // namespace Tempus
381#endif // Tempus_IntegratorForwardSensitivity_impl_hpp
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.
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...
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.