Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
20 namespace Tempus {
21 
22 template <class Scalar>
25  const Teuchos::RCP<IntegratorBasic<Scalar>> &integrator,
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 
38 template<class Scalar>
41 {
42  integrator_ = createIntegratorBasic<Scalar>();
43  integrator_->initialize();
44 }
45 
46 template<class Scalar>
50 {
51  if (use_combined_method_)
52  integrator_->setModel(sens_model_);
53  else
54  integrator_->setStepper(sens_stepper_);
55 }
56 
57 template<class Scalar>
62  Teuchos::RCP<const Thyra::VectorBase<Scalar> > xdotdot0,
65  Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > DxdotdotDp0)
66 {
67  using Teuchos::RCP;
68  using Teuchos::rcp_dynamic_cast;
70  using Thyra::assign;
71  using Thyra::createMember;
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 
120 template<class Scalar>
123 getX() const
124 {
125  using Teuchos::RCP;
126  using Teuchos::rcp_dynamic_cast;
128 
129  RCP<const DMVPV> X = rcp_dynamic_cast<const DMVPV>(integrator_->getX());
130  return X->getMultiVector()->col(0);
131 }
132 
133 template<class Scalar>
136 getDxDp() const
137 {
138  using Teuchos::RCP;
139  using Teuchos::rcp_dynamic_cast;
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 
148 template<class Scalar>
151 getXDot() const
152 {
153  using Teuchos::RCP;
154  using Teuchos::rcp_dynamic_cast;
156 
157  RCP<const DMVPV> Xdot = rcp_dynamic_cast<const DMVPV>(integrator_->getXDot());
158  return Xdot->getMultiVector()->col(0);
159 }
160 
161 template<class Scalar>
164 getDXDotDp() const
165 {
166  using Teuchos::RCP;
167  using Teuchos::rcp_dynamic_cast;
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 
176 template<class Scalar>
179 getXDotDot() const
180 {
181  using Teuchos::RCP;
182  using Teuchos::rcp_dynamic_cast;
184 
185  RCP<const DMVPV> Xdotdot =
186  rcp_dynamic_cast<const DMVPV>(integrator_->getXDotDot());
187  return Xdotdot->getMultiVector()->col(0);
188 }
189 
190 template<class Scalar>
194 {
195  using Teuchos::RCP;
196  using Teuchos::rcp_dynamic_cast;
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 
206 template<class Scalar>
209 getG() const
210 {
211  typedef Thyra::ModelEvaluatorBase MEB;
212 
213  // Compute g which is computed by response 1 of the
214  // sensitivity model evaluator
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 
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 
237 template<class Scalar>
240 getDgDp() const
241 {
242  typedef Thyra::ModelEvaluatorBase MEB;
244 
245  // Compute final dg/dp which is computed by response 0 of the
246  // sensitivity model evaluator
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 
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 
270 template<class Scalar>
271 std::string
273 description() const
274 {
275  std::string name = "Tempus::IntegratorForwardSensitivity";
276  return(name);
277 }
278 
279 template<class Scalar>
280 void
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 
294 template <class Scalar>
297 getStepMode() const
298 {
299  if (use_combined_method_)
301  return sens_stepper_->getStepMode(); // Staggered case
302 }
303 
305 template<class Scalar>
310  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& sens_residual_model,
311  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& sens_solve_model)
312 {
313 
314 
317 
318  //1. create integrator
319  auto fwd_integrator = createIntegratorBasic<Scalar>(pList, model);
320 
321  //2. set parameter list
323  {
324  Teuchos::RCP<const Teuchos::ParameterList> integrator_pl = fwd_integrator->getValidParameters();
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");
350  sens_model = wrapCombinedFSAModelEvaluator(
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
365  Teuchos::rcp(new IntegratorForwardSensitivity<Scalar>(model, fwd_integrator, sens_model, sens_stepper, use_combined_method));
366  return (integrator);
367 }
368 
370 template<class Scalar>
373 {
374 
377  return(integrator);
378 }
379 
380 } // namespace Tempus
381 #endif // Tempus_IntegratorForwardSensitivity_impl_hpp
Teuchos::RCP< IntegratorBasic< Scalar > > integrator_
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...
virtual void initializeSolutionHistory(Teuchos::RCP< SolutionState< Scalar > > state=Teuchos::null)
Set the initial state which has the initial conditions.
SensitivityStepMode getStepMode() const
What mode the current time integration step is in.
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getDXDotDp() const
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 > > getDXDotDotDp() const
A stepper implementing staggered forward sensitivity analysis.
virtual void setStepper(Teuchos::RCP< Thyra::ModelEvaluator< Scalar > > model)
Set the Stepper.
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getDgDp() const
Return forward sensitivity stored in Jacobian format.
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getDxDp() const
Get the forward sensitivities .
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
RCP< const MultiVectorBase< Scalar > > getMultiVector() const
Teuchos::RCP< IntegratorForwardSensitivity< Scalar > > createIntegratorForwardSensitivity(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.
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const override
ParameterList & setParameters(const ParameterList &source)
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getG() const
Return response function g.
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getX() const
Get the current solution, x, only. If looking for the solution vector and the sensitivities, use SolutionState-&gt;getX() which will return a Block MultiVector with the first block containing the current solution, x, and the remaining blocks are the forward sensitivities .
RCP< std::basic_ostream< char_type, traits_type > > getOStream()
ParameterList & setParametersNotAlreadySet(const ParameterList &source)
RCP< MultiVectorBase< Scalar > > getNonconstMultiVector()
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
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...
A ModelEvaluator decorator for sensitivity analysis.
Time integrator implementing forward sensitivity analysis.
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)