Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros 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 
12 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
13 #include "Thyra_DefaultMultiVectorProductVector.hpp"
14 #include "Thyra_VectorStdOps.hpp"
15 #include "Thyra_MultiVectorStdOps.hpp"
16 #include "Tempus_CombinedForwardSensitivityModelEvaluator.hpp"
18 
19 namespace Tempus {
20 
21 template<class Scalar>
24  Teuchos::RCP<Teuchos::ParameterList> inputPL,
25  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model)
26 {
27  model_ = model;
28  integrator_ = integratorBasic<Scalar>();
29  this->setParameterList(inputPL);
30  createSensitivityModelAndStepper(model);
31  if (use_combined_method_)
32  integrator_ = integratorBasic<Scalar>(tempus_pl_, sens_model_);
33  else {
34  integrator_ = integratorBasic<Scalar>();
35  integrator_->setTempusParameterList(tempus_pl_);
36  integrator_->setStepperWStepper(sens_stepper_);
37  integrator_->initialize();
38  }
39 }
40 
41 template<class Scalar>
44  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model,
45  std::string stepperType)
46 {
47  model_ = model;
48  integrator_ = integratorBasic<Scalar>();
49  this->setParameterList(Teuchos::null);
50  createSensitivityModelAndStepper(model);
51  if (use_combined_method_)
52  integrator_ = integratorBasic<Scalar>(sens_model_, stepperType);
53  else {
54  integrator_ = integratorBasic<Scalar>();
55  integrator_->setParameterList(tempus_pl_);
56  integrator_->setStepperWStepper(sens_stepper_);
57  integrator_->initialize();
58  }
59 
60 }
61 
62 template<class Scalar>
65 {
66  integrator_ = integratorBasic<Scalar>();
67  this->setParameterList(Teuchos::null);
68 }
69 
70 template<class Scalar>
73  Teuchos::RCP<Thyra::ModelEvaluator<Scalar> > model)
74 {
75  createSensitivityModelAndStepper(model);
76  if (use_combined_method_)
77  integrator_->setStepper(sens_model_);
78  else
79  integrator_->setStepperWStepper(sens_stepper_);
80 }
81 
82 template<class Scalar>
85  Teuchos::RCP<const Thyra::VectorBase<Scalar> > x0,
86  Teuchos::RCP<const Thyra::VectorBase<Scalar> > xdot0,
87  Teuchos::RCP<const Thyra::VectorBase<Scalar> > xdotdot0,
88  Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > DxDp0,
89  Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > DxdotDp0,
90  Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > DxdotdotDp0)
91 {
92  using Teuchos::RCP;
93  using Teuchos::rcp_dynamic_cast;
94  using Thyra::VectorSpaceBase;
95  using Thyra::assign;
96  using Thyra::createMember;
97  typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
98 
99  //
100  // Create and initialize product X, Xdot, Xdotdot
101 
102  RCP< const VectorSpaceBase<Scalar> > space;
103  if (use_combined_method_)
104  space = sens_model_->get_x_space();
105  else
106  space = sens_stepper_->get_x_space();
107  RCP<DMVPV> X = rcp_dynamic_cast<DMVPV>(createMember(space));
108  RCP<DMVPV> Xdot = rcp_dynamic_cast<DMVPV>(createMember(space));
109  RCP<DMVPV> Xdotdot = rcp_dynamic_cast<DMVPV>(createMember(space));
110 
111  const int num_param = X->getNonconstMultiVector()->domain()->dim()-1;
112  const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
113  const Teuchos::Range1D rng(1,num_param);
114 
115  // x
116  assign(X->getNonconstMultiVector()->col(0).ptr(), *x0);
117  if (DxDp0 == Teuchos::null)
118  assign(X->getNonconstMultiVector()->subView(rng).ptr(), zero);
119  else
120  assign(X->getNonconstMultiVector()->subView(rng).ptr(), *DxDp0);
121 
122  // xdot
123  if (xdot0 == Teuchos::null)
124  assign(Xdot->getNonconstMultiVector()->col(0).ptr(), zero);
125  else
126  assign(Xdot->getNonconstMultiVector()->col(0).ptr(), *xdot0);
127  if (DxdotDp0 == Teuchos::null)
128  assign(Xdot->getNonconstMultiVector()->subView(rng).ptr(), zero);
129  else
130  assign(Xdot->getNonconstMultiVector()->subView(rng).ptr(), *DxdotDp0);
131 
132  // xdotdot
133  if (xdotdot0 == Teuchos::null)
134  assign(Xdotdot->getNonconstMultiVector()->col(0).ptr(), zero);
135  else
136  assign(Xdotdot->getNonconstMultiVector()->col(0).ptr(), *xdotdot0);
137  if (DxdotDp0 == Teuchos::null)
138  assign(Xdotdot->getNonconstMultiVector()->subView(rng).ptr(), zero);
139  else
140  assign(Xdotdot->getNonconstMultiVector()->subView(rng).ptr(), *DxdotdotDp0);
141 
142  integrator_->initializeSolutionHistory(t0, X, Xdot, Xdotdot);
143 }
144 
145 template<class Scalar>
146 Teuchos::RCP<const Thyra::VectorBase<Scalar> >
148 getX() const
149 {
150  using Teuchos::RCP;
151  using Teuchos::rcp_dynamic_cast;
152  typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
153 
154  RCP<const DMVPV> X = rcp_dynamic_cast<const DMVPV>(integrator_->getX());
155  return X->getMultiVector()->col(0);
156 }
157 
158 template<class Scalar>
159 Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
161 getDxDp() const
162 {
163  using Teuchos::RCP;
164  using Teuchos::rcp_dynamic_cast;
165  typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
166 
167  RCP<const DMVPV> X = rcp_dynamic_cast<const DMVPV>(integrator_->getX());
168  const int num_param = X->getMultiVector()->domain()->dim()-1;
169  const Teuchos::Range1D rng(1,num_param);
170  return X->getMultiVector()->subView(rng);
171 }
172 
173 template<class Scalar>
174 Teuchos::RCP<const Thyra::VectorBase<Scalar> >
176 getXdot() const
177 {
178  using Teuchos::RCP;
179  using Teuchos::rcp_dynamic_cast;
180  typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
181 
182  RCP<const DMVPV> Xdot = rcp_dynamic_cast<const DMVPV>(integrator_->getXdot());
183  return Xdot->getMultiVector()->col(0);
184 }
185 
186 template<class Scalar>
187 Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
189 getDxdotDp() const
190 {
191  using Teuchos::RCP;
192  using Teuchos::rcp_dynamic_cast;
193  typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
194 
195  RCP<const DMVPV> Xdot = rcp_dynamic_cast<const DMVPV>(integrator_->getXdot());
196  const int num_param = Xdot->getMultiVector()->domain()->dim()-1;
197  const Teuchos::Range1D rng(1,num_param);
198  return Xdot->getMultiVector()->subView(rng);
199 }
200 
201 template<class Scalar>
202 Teuchos::RCP<const Thyra::VectorBase<Scalar> >
204 getXdotdot() const
205 {
206  using Teuchos::RCP;
207  using Teuchos::rcp_dynamic_cast;
208  typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
209 
210  RCP<const DMVPV> Xdotdot =
211  rcp_dynamic_cast<const DMVPV>(integrator_->getXdotdot());
212  return Xdotdot->getMultiVector()->col(0);
213 }
214 
215 template<class Scalar>
216 Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
219 {
220  using Teuchos::RCP;
221  using Teuchos::rcp_dynamic_cast;
222  typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
223 
224  RCP<const DMVPV> Xdotdot =
225  rcp_dynamic_cast<const DMVPV>(integrator_->getXdotdot());
226  const int num_param = Xdotdot->getMultiVector()->domain()->dim()-1;
227  const Teuchos::Range1D rng(1,num_param);
228  return Xdotdot->getMultiVector()->subView(rng);
229 }
230 
231 template<class Scalar>
232 std::string
234 description() const
235 {
236  std::string name = "Tempus::IntegratorForwardSensitivity";
237  return(name);
238 }
239 
240 template<class Scalar>
241 void
244  Teuchos::FancyOStream &out,
245  const Teuchos::EVerbosityLevel verbLevel) const
246 {
247  out << description() << "::describe" << std::endl;
248  integrator_->describe(out, verbLevel);
249 }
250 
251 template<class Scalar>
252 void
254 setParameterList(const Teuchos::RCP<Teuchos::ParameterList> & inputPL)
255 {
256  tempus_pl_ = Teuchos::parameterList();
257  if (inputPL != Teuchos::null)
258  *tempus_pl_ = *inputPL;
259  tempus_pl_->setParametersNotAlreadySet(*this->getValidParameters());
260  sens_pl_ = Teuchos::sublist(tempus_pl_, "Sensitivities", false);
261  std::string integratorName =
262  tempus_pl_->get<std::string>("Integrator Name", "Default Integrator");
263  std::string stepperName =
264  tempus_pl_->sublist(integratorName).get<std::string>("Stepper Name");
265  stepper_pl_ = Teuchos::sublist(tempus_pl_, stepperName, true);
266  use_combined_method_ =
267  sens_pl_->get<std::string>("Sensitivity Method") == "Combined";
268 }
269 
270 template<class Scalar>
271 Teuchos::RCP<Teuchos::ParameterList>
274 {
275  Teuchos::RCP<Teuchos::ParameterList> temp_param_list = tempus_pl_;
276  tempus_pl_ = Teuchos::null;
277  sens_pl_ = Teuchos::null;
278  stepper_pl_ = Teuchos::null;
279  integrator_->unsetParameterList();
280  return temp_param_list;
281 }
282 
283 template<class Scalar>
284 Teuchos::RCP<const Teuchos::ParameterList>
287 {
288  Teuchos::RCP<Teuchos::ParameterList> pl =
289  Teuchos::rcp(new Teuchos::ParameterList);
290  Teuchos::RCP<const Teuchos::ParameterList> integrator_pl =
291  integrator_->getValidParameters();
292  Teuchos::RCP<const Teuchos::ParameterList> sensitivity_pl =
294  pl->setParameters(*integrator_pl);
295  Teuchos::ParameterList& spl = pl->sublist("Sensitivities");
296  spl.setParameters(*sensitivity_pl);
297  spl.set("Sensitivity Method", "Combined");
298  spl.set("Reuse State Linear Solver", false);
299 
300  return pl;
301 }
302 
303 template <class Scalar>
304 void
307  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& /* model */)
308 {
309  using Teuchos::rcp;
310 
311  Teuchos::RCP<Teuchos::ParameterList> spl = Teuchos::parameterList();
312  *spl = *sens_pl_;
313  spl->remove("Sensitivity Method");
314 
315  if (use_combined_method_) {
316  spl->remove("Reuse State Linear Solver");
317  sens_model_ =
318  wrapCombinedFSAModelEvaluator(model_, spl);
319  }
320  else {
321  sens_stepper_ =
323  model_, stepper_pl_, spl));
324  }
325 }
326 
327 /// Non-member constructor
328 template<class Scalar>
329 Teuchos::RCP<Tempus::IntegratorForwardSensitivity<Scalar> >
331  Teuchos::RCP<Teuchos::ParameterList> pList,
332  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model)
333 {
334  Teuchos::RCP<Tempus::IntegratorForwardSensitivity<Scalar> > integrator =
335  Teuchos::rcp(new Tempus::IntegratorForwardSensitivity<Scalar>(pList, model));
336  return(integrator);
337 }
338 
339 /// Non-member constructor
340 template<class Scalar>
341 Teuchos::RCP<Tempus::IntegratorForwardSensitivity<Scalar> >
343  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model,
344  std::string stepperType)
345 {
346  Teuchos::RCP<Tempus::IntegratorForwardSensitivity<Scalar> > integrator =
347  Teuchos::rcp(new Tempus::IntegratorForwardSensitivity<Scalar>(model, stepperType));
348  return(integrator);
349 }
350 
351 /// Non-member constructor
352 template<class Scalar>
353 Teuchos::RCP<Tempus::IntegratorForwardSensitivity<Scalar> >
355 {
356  Teuchos::RCP<Tempus::IntegratorForwardSensitivity<Scalar> > integrator =
358  return(integrator);
359 }
360 
361 } // namespace Tempus
362 #endif // Tempus_IntegratorForwardSensitivity_impl_hpp
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getXdot() const
Get current the time derivative of the solution, xdot.
virtual void initializeSolutionHistory(Teuchos::RCP< SolutionState< Scalar > > state=Teuchos::null)
Set the initial state which has the initial conditions.
A stepper implementing staggered forward sensitivity analysis.
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getDxdotDp() const
virtual void setStepper(Teuchos::RCP< Thyra::ModelEvaluator< Scalar > > model)
Set the Stepper.
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList() override
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getDxdotdotDp() const
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getXdotdot() const
Get current the second time derivative of the solution, xdotdot.
void createSensitivityModelAndStepper(const Teuchos::RCP< Thyra::ModelEvaluator< Scalar > > &model)
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getDxDp() const
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &pl) override
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const override
Teuchos::RCP< Tempus::IntegratorForwardSensitivity< Scalar > > integratorForwardSensitivity(Teuchos::RCP< Teuchos::ParameterList > pList, const Teuchos::RCP< Thyra::ModelEvaluator< Scalar > > &model)
Non-member constructor.
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getX() const
Get current the solution, x.
Teuchos::RCP< SensitivityModelEvaluatorBase< Scalar > > wrapCombinedFSAModelEvaluator(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, const Teuchos::RCP< const Teuchos::ParameterList > &pList=Teuchos::null)
Time integrator implementing forward sensitivity analysis.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const override