Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tempus_IntegratorPseudoTransientForwardSensitivity_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_IntegratorPseudoTransientForwardSensitivity_impl_hpp
10 #define Tempus_IntegratorPseudoTransientForwardSensitivity_impl_hpp
11 
12 #include "Thyra_DefaultMultiVectorProductVector.hpp"
13 #include "Thyra_VectorStdOps.hpp"
14 #include "Thyra_MultiVectorStdOps.hpp"
15 
17 
18 
19 namespace Tempus {
20 
21 template<class Scalar>
26  reuse_solver_(false)
27 {
28  model_ = model;
30  state_integrator_ = integratorBasic<Scalar>(inputPL, model_);
31  sens_integrator_ = integratorBasic<Scalar>(inputPL, sens_model_);
32 }
33 
34 template<class Scalar>
38  std::string stepperType) :
39  reuse_solver_(false),
40  force_W_update_(false)
41 {
42  model_ = model;
43  sens_model_ = createSensitivityModel(model, Teuchos::null);
44  state_integrator_ = integratorBasic<Scalar>(model_, stepperType);
45  sens_integrator_ = integratorBasic<Scalar>(sens_model_, stepperType);
46 }
47 
48 template<class Scalar>
51  reuse_solver_(false),
52  force_W_update_(false)
53 {
54  state_integrator_ = integratorBasic<Scalar>();
55  sens_integrator_ = integratorBasic<Scalar>();
56 }
57 
58 template<class Scalar>
59 bool
62 {
63  using Teuchos::RCP;
64  using Thyra::VectorBase;
65 
66  // Run state integrator and get solution
67  bool state_status = state_integrator_->advanceTime();
68 
69  // Set solution in sensitivity ME
70  sens_model_->setForwardSolutionState(state_integrator_->getCurrentState());
71 
72  // Reuse state solver if requested
73  if (reuse_solver_ &&
74  state_integrator_->getStepper()->getSolver() != Teuchos::null) {
75  sens_model_->setSolver(state_integrator_->getStepper()->getSolver(),
76  force_W_update_);
77  }
78 
79  // Run sensitivity integrator
80  bool sens_status = sens_integrator_->advanceTime();
81 
82  buildSolutionHistory();
83 
84  return state_status && sens_status;
85 }
86 
87 template<class Scalar>
88 bool
90 advanceTime(const Scalar timeFinal)
91 {
92  using Teuchos::RCP;
93  using Thyra::VectorBase;
94 
95  // Run state integrator and get solution
96  bool state_status = state_integrator_->advanceTime(timeFinal);
97 
98  // Set solution in sensitivity ME
99  sens_model_->setForwardSolutionState(state_integrator_->getCurrentState());
100 
101  // Reuse state solver if requested
102  if (reuse_solver_ &&
103  state_integrator_->getStepper()->getSolver() != Teuchos::null) {
104  sens_model_->setSolver(state_integrator_->getStepper()->getSolver(),
105  force_W_update_);
106  }
107 
108  // Run sensitivity integrator
109  bool sens_status = sens_integrator_->advanceTime(timeFinal);
110 
111  buildSolutionHistory();
112 
113  return state_status && sens_status;
114 }
115 
116 template<class Scalar>
117 Scalar
119 getTime() const
120 {
121  return solutionHistory_->getCurrentTime();
122 }
123 
124 template<class Scalar>
125 int
127 getIndex() const
128 {
129  return solutionHistory_->getCurrentIndex();
130 }
131 
132 template<class Scalar>
133 Status
135 getStatus() const
136 {
137  Status state_status = state_integrator_->getStatus();
138  Status sens_status = sens_integrator_->getStatus();
139  if (state_status == FAILED || sens_status == FAILED)
140  return FAILED;
141  if (state_status == WORKING || sens_status == WORKING)
142  return WORKING;
143  return PASSED;
144 }
145 
146 template<class Scalar>
149 getStepper() const
150 {
151  return state_integrator_->getStepper();
152 }
153 
154 template<class Scalar>
158 {
159  return state_integrator_->getTempusParameterList();
160 }
161 
162 template<class Scalar>
163 void
166 {
167  state_integrator_->setTempusParameterList(pl);
168  sens_integrator_->setTempusParameterList(pl);
169 }
170 
171 template<class Scalar>
175 {
176  return solutionHistory_;
177 }
178 
179 template<class Scalar>
183 {
184  return state_integrator_->getTimeStepControl();
185 }
186 
187 template<class Scalar>
191 {
192  return state_integrator_->getNonConstTimeStepControl();
193 }
194 
195 template<class Scalar>
200  Teuchos::RCP<const Thyra::VectorBase<Scalar> > xdotdot0,
203  Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > DxdotdotDp0)
204 {
205  using Teuchos::RCP;
206  using Teuchos::rcp_dynamic_cast;
208  using Thyra::assign;
209  using Thyra::createMember;
211 
212  //
213  // Create and initialize product X, Xdot, Xdotdot
214 
215  RCP< const VectorSpaceBase<Scalar> > space = sens_model_->get_x_space();
216  RCP<DMVPV> X = rcp_dynamic_cast<DMVPV>(createMember(space));
217  RCP<DMVPV> Xdot = rcp_dynamic_cast<DMVPV>(createMember(space));
218  RCP<DMVPV> Xdotdot = rcp_dynamic_cast<DMVPV>(createMember(space));
219  const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
220 
221  // x
222  if (DxDp0 == Teuchos::null)
223  assign(X->getNonconstMultiVector().ptr(), zero);
224  else
225  assign(X->getNonconstMultiVector().ptr(), *DxDp0);
226 
227  // xdot
228  if (DxdotDp0 == Teuchos::null)
229  assign(Xdot->getNonconstMultiVector().ptr(), zero);
230  else
231  assign(Xdot->getNonconstMultiVector().ptr(), *DxdotDp0);
232 
233  // xdotdot
234  if (DxdotDp0 == Teuchos::null)
235  assign(Xdotdot->getNonconstMultiVector().ptr(), zero);
236  else
237  assign(Xdotdot->getNonconstMultiVector().ptr(), *DxdotdotDp0);
238 
239  state_integrator_->initializeSolutionHistory(t0, x0, xdot0, xdotdot0);
240  sens_integrator_->initializeSolutionHistory(t0, X, Xdot, Xdotdot);
241 }
242 
243 template<class Scalar>
246 getX() const
247 {
248  using Teuchos::RCP;
249  using Teuchos::rcp_dynamic_cast;
251 
252  RCP<const DMVPV> X =
253  rcp_dynamic_cast<const DMVPV>(solutionHistory_->getCurrentState()->getX());
254  return X->getMultiVector()->col(0);
255 }
256 
257 template<class Scalar>
260 getDxDp() const
261 {
262  using Teuchos::RCP;
263  using Teuchos::rcp_dynamic_cast;
265 
266  RCP<const DMVPV> X =
267  rcp_dynamic_cast<const DMVPV>(solutionHistory_->getCurrentState()->getX());
268  const int num_param = X->getMultiVector()->domain()->dim()-1;
269  const Teuchos::Range1D rng(1,num_param);
270  return X->getMultiVector()->subView(rng);
271 }
272 
273 template<class Scalar>
276 getXDot() const
277 {
278  using Teuchos::RCP;
279  using Teuchos::rcp_dynamic_cast;
281 
282  RCP<const DMVPV> Xdot =
283  rcp_dynamic_cast<const DMVPV>(solutionHistory_->getCurrentState()->getXDot());
284  return Xdot->getMultiVector()->col(0);
285 }
286 
287 template<class Scalar>
290 getDXDotDp() const
291 {
292  using Teuchos::RCP;
293  using Teuchos::rcp_dynamic_cast;
295 
296  RCP<const DMVPV> Xdot =
297  rcp_dynamic_cast<const DMVPV>(solutionHistory_->getCurrentState()->getXDot());
298  const int num_param = Xdot->getMultiVector()->domain()->dim()-1;
299  const Teuchos::Range1D rng(1,num_param);
300  return Xdot->getMultiVector()->subView(rng);
301 }
302 
303 template<class Scalar>
306 getXDotDot() const
307 {
308  using Teuchos::RCP;
309  using Teuchos::rcp_dynamic_cast;
311 
312  RCP<const DMVPV> Xdotdot =
313  rcp_dynamic_cast<const DMVPV>(solutionHistory_->getCurrentState()->getXDotDot());
314  return Xdotdot->getMultiVector()->col(0);
315 }
316 
317 template<class Scalar>
321 {
322  using Teuchos::RCP;
323  using Teuchos::rcp_dynamic_cast;
325 
326  RCP<const DMVPV> Xdotdot =
327  rcp_dynamic_cast<const DMVPV>(solutionHistory_->getCurrentState()->getXDotDot());
328  const int num_param = Xdotdot->getMultiVector()->domain()->dim()-1;
329  const Teuchos::Range1D rng(1,num_param);
330  return Xdotdot->getMultiVector()->subView(rng);
331 }
332 
333 template<class Scalar>
334 std::string
336 description() const
337 {
338  std::string name = "Tempus::IntegratorPseudoTransientForwardSensitivity";
339  return(name);
340 }
341 
342 template<class Scalar>
343 void
346  Teuchos::FancyOStream &in_out,
347  const Teuchos::EVerbosityLevel verbLevel) const
348 {
349  auto out = Teuchos::fancyOStream( in_out.getOStream() );
350  out->setOutputToRootOnly(0);
351  *out << description() << "::describe" << std::endl;
352  state_integrator_->describe(in_out, verbLevel);
353  sens_integrator_->describe(in_out, verbLevel);
354 }
355 
356 template<class Scalar>
357 void
360 {
361  state_integrator_->setParameterList(inputPL);
362  sens_integrator_->setParameterList(inputPL);
363  reuse_solver_ =
364  inputPL->sublist("Sensitivities").get("Reuse State Linear Solver", false);
365  force_W_update_ =
366  inputPL->sublist("Sensitivities").get("Force W Update", false);
367 }
368 
369 template<class Scalar>
373 {
374  state_integrator_->unsetParameterList();
375  return sens_integrator_->unsetParameterList();
376 }
377 
378 template<class Scalar>
382 {
386  state_integrator_->getValidParameters();
389  pl->setParameters(*integrator_pl);
390  pl->sublist("Sensitivities").setParameters(*sensitivity_pl);
391  pl->sublist("Sensitivities").set("Reuse State Linear Solver", false);
392  pl->sublist("Sensitivities").set("Force W Update", false);
393 
394  return pl;
395 }
396 
397 template<class Scalar>
401 {
402  return state_integrator_->getNonconstParameterList();
403 }
404 
405 template <class Scalar>
411 {
412  using Teuchos::rcp;
413 
414  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
415  if (inputPL != Teuchos::null) {
416  *pl = inputPL->sublist("Sensitivities");
417  }
418  reuse_solver_ = pl->get("Reuse State Linear Solver", false);
419  force_W_update_ = pl->get("Force W Update", true);
420  pl->remove("Reuse State Linear Solver");
421  pl->remove("Force W Update");
422  return wrapStaggeredFSAModelEvaluator(model, pl);
423 }
424 
425 template<class Scalar>
426 void
429 {
430  using Teuchos::RCP;
431  using Teuchos::rcp;
432  using Teuchos::rcp_dynamic_cast;
434  using Thyra::VectorBase;
437  using Thyra::createMembers;
438  using Thyra::multiVectorProductVector;
439  using Thyra::assign;
441 
442  // Create combined solution histories, first for the states with zero
443  // sensitivities and then for the sensitivities with frozen states
444  RCP<ParameterList> shPL =
445  Teuchos::sublist(state_integrator_->getIntegratorParameterList(),
446  "Solution History", true);
447  solutionHistory_ = createSolutionHistoryPL<Scalar>(shPL);
448 
449  const int num_param =
450  rcp_dynamic_cast<const DMVPV>(sens_integrator_->getX())->getMultiVector()->domain()->dim();
451  RCP<const VectorSpaceBase<Scalar> > x_space = model_->get_x_space();
452  RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> > prod_space =
453  Thyra::multiVectorProductVectorSpace(x_space, num_param+1);
454  const Teuchos::Range1D rng(1,num_param);
455  const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
456 
457  RCP<const SolutionHistory<Scalar> > state_solution_history =
458  state_integrator_->getSolutionHistory();
459  int num_states = state_solution_history->getNumStates();
460  for (int i=0; i<num_states; ++i) {
461  RCP<const SolutionState<Scalar> > state = (*state_solution_history)[i];
462 
463  // X
464  RCP< MultiVectorBase<Scalar> > x_mv =
465  createMembers(x_space, num_param+1);
466  assign(x_mv->col(0).ptr(), *(state->getX()));
467  assign(x_mv->subView(rng).ptr(), zero);
468  RCP<VectorBase<Scalar> > x = multiVectorProductVector(prod_space, x_mv);
469 
470  // X-Dot
471  RCP<VectorBase<Scalar> > x_dot;
472  if (state->getXDot() != Teuchos::null) {
473  RCP< MultiVectorBase<Scalar> > x_dot_mv =
474  createMembers(x_space, num_param+1);
475  assign(x_dot_mv->col(0).ptr(), *(state->getXDot()));
476  assign(x_dot_mv->subView(rng).ptr(), zero);
477  x_dot = multiVectorProductVector(prod_space, x_dot_mv);
478  }
479 
480  // X-Dot-Dot
481  RCP<VectorBase<Scalar> > x_dot_dot;
482  if (state->getXDotDot() != Teuchos::null) {
483  RCP< MultiVectorBase<Scalar> > x_dot_dot_mv =
484  createMembers(x_space, num_param+1);
485  assign(x_dot_dot_mv->col(0).ptr(), *(state->getXDotDot()));
486  assign(x_dot_dot_mv->subView(rng).ptr(), zero);
487  x_dot_dot = multiVectorProductVector(prod_space, x_dot_dot_mv);
488  }
489 
490  RCP<SolutionState<Scalar> > prod_state = state->clone();
491  prod_state->setX(x);
492  prod_state->setXDot(x_dot);
493  prod_state->setXDotDot(x_dot_dot);
494  solutionHistory_->addState(prod_state);
495  }
496 
497  RCP<const VectorBase<Scalar> > frozen_x =
498  state_solution_history->getCurrentState()->getX();
499  RCP<const VectorBase<Scalar> > frozen_x_dot =
500  state_solution_history->getCurrentState()->getXDot();
501  RCP<const VectorBase<Scalar> > frozen_x_dot_dot =
502  state_solution_history->getCurrentState()->getXDotDot();
503  RCP<const SolutionHistory<Scalar> > sens_solution_history =
504  sens_integrator_->getSolutionHistory();
505  num_states = sens_solution_history->getNumStates();
506  for (int i=0; i<num_states; ++i) {
507  RCP<const SolutionState<Scalar> > state = (*sens_solution_history)[i];
508 
509  // X
510  RCP< MultiVectorBase<Scalar> > x_mv =
511  createMembers(x_space, num_param+1);
512  RCP<const MultiVectorBase<Scalar> > dxdp =
513  rcp_dynamic_cast<const DMVPV>(state->getX())->getMultiVector();
514  assign(x_mv->col(0).ptr(), *(frozen_x));
515  assign(x_mv->subView(rng).ptr(), *dxdp);
516  RCP<VectorBase<Scalar> > x = multiVectorProductVector(prod_space, x_mv);
517 
518  // X-Dot
519  RCP<VectorBase<Scalar> > x_dot;
520  if (state->getXDot() != Teuchos::null) {
521  RCP< MultiVectorBase<Scalar> > x_dot_mv =
522  createMembers(x_space, num_param+1);
523  RCP<const MultiVectorBase<Scalar> > dxdotdp =
524  rcp_dynamic_cast<const DMVPV>(state->getXDot())->getMultiVector();
525  assign(x_dot_mv->col(0).ptr(), *(frozen_x_dot));
526  assign(x_dot_mv->subView(rng).ptr(), *dxdotdp);
527  x_dot = multiVectorProductVector(prod_space, x_dot_mv);
528  }
529 
530  // X-Dot-Dot
531  RCP<VectorBase<Scalar> > x_dot_dot;
532  if (state->getXDotDot() != Teuchos::null) {
533  RCP< MultiVectorBase<Scalar> > x_dot_dot_mv =
534  createMembers(x_space, num_param+1);
535  RCP<const MultiVectorBase<Scalar> > dxdotdotdp =
536  rcp_dynamic_cast<const DMVPV>(state->getXDotDot())->getMultiVector();
537  assign(x_dot_dot_mv->col(0).ptr(), *(frozen_x_dot_dot));
538  assign(x_dot_dot_mv->subView(rng).ptr(), *dxdotdotdp);
539  x_dot_dot = multiVectorProductVector(prod_space, x_dot_dot_mv);
540  }
541 
542  RCP<SolutionState<Scalar> > prod_state = state->clone();
543  prod_state->setX(x);
544  prod_state->setXDot(x_dot);
545  prod_state->setXDotDot(x_dot_dot);
546  solutionHistory_->addState(prod_state, false);
547  }
548 }
549 
551 template<class Scalar>
556 {
559  return(integrator);
560 }
561 
563 template<class Scalar>
567  std::string stepperType)
568 {
571  return(integrator);
572 }
573 
575 template<class Scalar>
578 {
581  return(integrator);
582 }
583 
584 } // namespace Tempus
585 #endif // Tempus_IntegratorPseudoTransientForwardSensitivity_impl_hpp
virtual Teuchos::RCP< Stepper< Scalar > > getStepper() const override
Get the Stepper.
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &pl) override
T & get(const std::string &name, T def_value)
virtual Teuchos::RCP< const SolutionHistory< Scalar > > getSolutionHistory() const override
Get the SolutionHistory.
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 > > getDxDp() const
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
Teuchos::RCP< Tempus::IntegratorPseudoTransientForwardSensitivity< Scalar > > integratorPseudoTransientForwardSensitivity(Teuchos::RCP< Teuchos::ParameterList > pList, const Teuchos::RCP< Thyra::ModelEvaluator< Scalar > > &model)
Nonmember constructor.
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.
bool remove(std::string const &name, bool throwIfNotExists=true)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Status
Status for the Integrator, the Stepper and the SolutionState.
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getXDotDot() const
Get current the second time derivative of the solution, xdotdot.
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const override
RCP< const MultiVectorBase< Scalar > > getMultiVector() const
virtual Teuchos::RCP< Teuchos::ParameterList > getTempusParameterList() override
Return a copy of the Tempus ParameterList.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const override
ParameterList & setParameters(const ParameterList &source)
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getDXDotDp() const
Teuchos::RCP< SensitivityModelEvaluatorBase< Scalar > > createSensitivityModel(const Teuchos::RCP< Thyra::ModelEvaluator< Scalar > > &model, const Teuchos::RCP< Teuchos::ParameterList > &inputPL)
RCP< std::basic_ostream< char_type, traits_type > > getOStream()
Time integrator suitable for pseudotransient forward sensitivity analysis.
virtual RCP< const VectorSpaceBase< Scalar > > domain() const
virtual void setTempusParameterList(Teuchos::RCP< Teuchos::ParameterList > pl) override
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
virtual Teuchos::RCP< const TimeStepControl< Scalar > > getTimeStepControl() const override
Get the TimeStepControl.
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getDXDotDotDp() const
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getXDot() const
Get current the time derivative of the solution, xdot.
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getX() const
Get current the solution, x.