Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Tempus_IntegratorPseudoTransientAdjointSensitivity_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_IntegratorPseudoTransientAdjointSensitivity_impl_hpp
10 #define Tempus_IntegratorPseudoTransientAdjointSensitivity_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 "NOX_Thyra.H"
17 
18 namespace Tempus {
19 
20 template<class Scalar>
23  Teuchos::RCP<Teuchos::ParameterList> inputPL,
24  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model)
25 {
26  model_ = model;
27  state_integrator_ = integratorBasic<Scalar>(inputPL, model_);
28  sens_model_ = createSensitivityModel(model_, inputPL);
29  sens_integrator_ = integratorBasic<Scalar>(inputPL, sens_model_);
30 }
31 
32 template<class Scalar>
35  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model,
36  std::string stepperType)
37 {
38  model_ = model;
39  state_integrator_ = integratorBasic<Scalar>(model_, stepperType);
40  sens_model_ = createSensitivityModel(model, Teuchos::null);
41  sens_integrator_ = integratorBasic<Scalar>(sens_model_, stepperType);
42 }
43 
44 template<class Scalar>
47 {
48  state_integrator_ = integratorBasic<Scalar>();
49  sens_integrator_ = integratorBasic<Scalar>();
50 }
51 
52 template<class Scalar>
53 bool
56 {
57  const Scalar tfinal =
58  state_integrator_->getTimeStepControl()->getFinalTime();
59  return advanceTime(tfinal);
60 }
61 
62 template<class Scalar>
63 bool
65 advanceTime(const Scalar timeFinal)
66 {
67  using Teuchos::RCP;
68  using Thyra::VectorBase;
69  typedef Thyra::ModelEvaluatorBase MEB;
70 
71  // Run state integrator and get solution
72  bool state_status = state_integrator_->advanceTime(timeFinal);
73 
74  // Set solution in sensitivity ME
75  sens_model_->setForwardSolutionHistory(
76  state_integrator_->getSolutionHistory());
77 
78  // Run sensitivity integrator
79  bool sens_status = sens_integrator_->advanceTime(timeFinal);
80 
81  // Compute final dg/dp which is computed by response 0 of the adjoint
82  // model evaluator
83  MEB::InArgs<Scalar> inargs = sens_model_->getNominalValues();
84  MEB::OutArgs<Scalar> outargs = sens_model_->createOutArgs();
85  inargs.set_t(sens_integrator_->getTime());
86  inargs.set_x(sens_integrator_->getX());
87  if (inargs.supports(MEB::IN_ARG_x_dot))
88  inargs.set_x_dot(sens_integrator_->getXdot());
89  if (inargs.supports(MEB::IN_ARG_x_dot_dot))
90  inargs.set_x_dot_dot(sens_integrator_->getXdotdot());
91  RCP<VectorBase<Scalar> > G = dgdp_;
92  if (G == Teuchos::null) {
93  G = Thyra::createMember(sens_model_->get_g_space(0));
94  dgdp_ = Teuchos::rcp_dynamic_cast<DMVPV>(G);
95  }
96  outargs.set_g(0, G);
97  sens_model_->evalModel(inargs, outargs);
98 
99  buildSolutionHistory();
100 
101  return state_status && sens_status;
102 }
103 
104 template<class Scalar>
105 Scalar
107 getTime() const
108 {
109  return solutionHistory_->getCurrentTime();
110 }
111 
112 template<class Scalar>
113 int
115 getIndex() const
116 {
117  return solutionHistory_->getCurrentIndex();
118 }
119 
120 template<class Scalar>
121 Status
123 getStatus() const
124 {
125  Status state_status = state_integrator_->getStatus();
126  Status sens_status = sens_integrator_->getStatus();
127  if (state_status == FAILED || sens_status == FAILED)
128  return FAILED;
129  if (state_status == WORKING || sens_status == WORKING)
130  return WORKING;
131  return PASSED;
132 }
133 
134 template<class Scalar>
135 Teuchos::RCP<Stepper<Scalar> >
137 getStepper() const
138 {
139  return state_integrator_->getStepper();
140 }
141 
142 template<class Scalar>
143 Teuchos::RCP<Teuchos::ParameterList>
146 {
147  return state_integrator_->getTempusParameterList();
148 }
149 
150 template<class Scalar>
151 void
153 setTempusParameterList(Teuchos::RCP<Teuchos::ParameterList> pl)
154 {
155  state_integrator_->setTempusParameterList(pl);
156  sens_integrator_->setTempusParameterList(pl);
157 }
158 
159 template<class Scalar>
160 Teuchos::RCP<const SolutionHistory<Scalar> >
163 {
164  return solutionHistory_;
165 }
166 
167 template<class Scalar>
168 Teuchos::RCP<const TimeStepControl<Scalar> >
171 {
172  return state_integrator_->getTimeStepControl();
173 }
174 
175 template<class Scalar>
176 Teuchos::RCP<TimeStepControl<Scalar> >
179 {
180  return state_integrator_->getNonConstTimeStepControl();
181 }
182 
183 template<class Scalar>
186  Teuchos::RCP<const Thyra::VectorBase<Scalar> > x0,
187  Teuchos::RCP<const Thyra::VectorBase<Scalar> > xdot0,
188  Teuchos::RCP<const Thyra::VectorBase<Scalar> > xdotdot0,
189  Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > y0,
190  Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > ydot0,
191  Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > ydotdot0)
192 {
193  using Teuchos::RCP;
194  using Teuchos::rcp_dynamic_cast;
195  using Thyra::VectorSpaceBase;
196  using Thyra::assign;
197  using Thyra::createMember;
198 
199  //
200  // Create and initialize product X, Xdot, Xdotdot
201 
202  RCP< const VectorSpaceBase<Scalar> > space = sens_model_->get_x_space();
203  RCP<DMVPV> Y = rcp_dynamic_cast<DMVPV>(createMember(space));
204  RCP<DMVPV> Ydot = rcp_dynamic_cast<DMVPV>(createMember(space));
205  RCP<DMVPV> Ydotdot = rcp_dynamic_cast<DMVPV>(createMember(space));
206  const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
207 
208  // y
209  if (y0 == Teuchos::null)
210  assign(Y->getNonconstMultiVector().ptr(), zero);
211  else
212  assign(Y->getNonconstMultiVector().ptr(), *y0);
213 
214  // ydot
215  if (ydot0 == Teuchos::null)
216  assign(Ydot->getNonconstMultiVector().ptr(), zero);
217  else
218  assign(Ydot->getNonconstMultiVector().ptr(), *ydot0);
219 
220  // ydotdot
221  if (ydotdot0 == Teuchos::null)
222  assign(Ydotdot->getNonconstMultiVector().ptr(), zero);
223  else
224  assign(Ydotdot->getNonconstMultiVector().ptr(), *ydotdot0);
225 
226  state_integrator_->initializeSolutionHistory(t0, x0, xdot0, xdotdot0);
227  sens_integrator_->initializeSolutionHistory(t0, Y, Ydot, Ydotdot);
228 }
229 
230 template<class Scalar>
231 Teuchos::RCP<const Thyra::VectorBase<Scalar> >
233 getX() const
234 {
235  return state_integrator_->getX();
236 }
237 
238 template<class Scalar>
239 Teuchos::RCP<const Thyra::VectorBase<Scalar> >
241 getXdot() const
242 {
243  return state_integrator_->getXdot();
244 }
245 
246 template<class Scalar>
247 Teuchos::RCP<const Thyra::VectorBase<Scalar> >
249 getXdotdot() const
250 {
251  return state_integrator_->getXdotdot();
252 }
253 
254 template<class Scalar>
255 Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
257 getDgDp() const
258 {
259  return dgdp_->getMultiVector();
260 }
261 
262 template<class Scalar>
263 std::string
265 description() const
266 {
267  std::string name = "Tempus::IntegratorPseudoTransientAdjointSensitivity";
268  return(name);
269 }
270 
271 template<class Scalar>
272 void
275  Teuchos::FancyOStream &out,
276  const Teuchos::EVerbosityLevel verbLevel) const
277 {
278  out << description() << "::describe" << std::endl;
279  state_integrator_->describe(out, verbLevel);
280  sens_integrator_->describe(out, verbLevel);
281 }
282 
283 template<class Scalar>
284 void
286 setParameterList(const Teuchos::RCP<Teuchos::ParameterList> & inputPL)
287 {
288  state_integrator_->setParameterList(inputPL);
289  sens_integrator_->setParameterList(inputPL);
290 }
291 
292 template<class Scalar>
293 Teuchos::RCP<Teuchos::ParameterList>
296 {
297  state_integrator_->unsetParameterList();
298  return sens_integrator_->unsetParameterList();
299 }
300 
301 template<class Scalar>
302 Teuchos::RCP<const Teuchos::ParameterList>
305 {
306  Teuchos::RCP<Teuchos::ParameterList> pl =
307  Teuchos::rcp(new Teuchos::ParameterList);
308  Teuchos::RCP<const Teuchos::ParameterList> integrator_pl =
309  state_integrator_->getValidParameters();
310  Teuchos::RCP<const Teuchos::ParameterList> sensitivity_pl =
312  pl->setParameters(*integrator_pl);
313  pl->sublist("Sensitivities").setParameters(*sensitivity_pl);
314 
315  return pl;
316 }
317 
318 template<class Scalar>
319 Teuchos::RCP<Teuchos::ParameterList>
322 {
323  return state_integrator_->getNonconstParameterList();
324 }
325 
326 template <class Scalar>
327 Teuchos::RCP<Tempus::AdjointSensitivityModelEvaluator<Scalar> >
330  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model,
331  const Teuchos::RCP<Teuchos::ParameterList>& inputPL)
332 {
333  using Teuchos::rcp;
334 
335  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
336  if (inputPL != Teuchos::null) {
337  *pl = inputPL->sublist("Sensitivities");
338  }
339  const Scalar tfinal = state_integrator_->getTimeStepControl()->getFinalTime();
341  model, tfinal, true, pl));
342 }
343 
344 template<class Scalar>
345 void
348 {
349  using Teuchos::RCP;
350  using Teuchos::rcp;
351  using Teuchos::rcp_dynamic_cast;
352  using Teuchos::ParameterList;
353  using Thyra::VectorBase;
354  using Thyra::VectorSpaceBase;
355  using Thyra::assign;
356  using Thyra::defaultProductVector;
357  typedef Thyra::DefaultProductVector<Scalar> DPV;
358  typedef Thyra::DefaultProductVectorSpace<Scalar> DPVS;
359 
360  // Create combined solution histories, first for the states with zero
361  // adjoint and then for the adjoint with frozen states
362  RCP<ParameterList> shPL =
363  Teuchos::sublist(state_integrator_->getIntegratorParameterList(),
364  "Solution History", true);
365  solutionHistory_ = rcp(new SolutionHistory<Scalar>(shPL));
366 
367  RCP<const VectorSpaceBase<Scalar> > x_space =
368  model_->get_x_space();
369  RCP<const VectorSpaceBase<Scalar> > adjoint_space =
370  sens_model_->get_x_space();
371  Teuchos::Array< RCP<const VectorSpaceBase<Scalar> > > spaces(2);
372  spaces[0] = x_space;
373  spaces[1] = adjoint_space;
374  RCP<const DPVS > prod_space = Thyra::productVectorSpace(spaces());
375  const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
376 
377  RCP<const SolutionHistory<Scalar> > state_solution_history =
378  state_integrator_->getSolutionHistory();
379  int num_states = state_solution_history->getNumStates();
380  for (int i=0; i<num_states; ++i) {
381  RCP<const SolutionState<Scalar> > state = (*state_solution_history)[i];
382 
383  // X
384  RCP<DPV> x = defaultProductVector(prod_space);
385  assign(x->getNonconstVectorBlock(0).ptr(), *(state->getX()));
386  assign(x->getNonconstVectorBlock(1).ptr(), zero);
387  RCP<VectorBase<Scalar> > x_b = x;
388 
389  // X-Dot
390  RCP<VectorBase<Scalar> > x_dot_b;
391  if (state->getXDot() != Teuchos::null) {
392  RCP<DPV> x_dot = defaultProductVector(prod_space);
393  assign(x_dot->getNonconstVectorBlock(0).ptr(), *(state->getXDot()));
394  assign(x_dot->getNonconstVectorBlock(1).ptr(), zero);
395  x_dot_b = x_dot;
396  }
397 
398  // X-Dot-Dot
399  RCP<VectorBase<Scalar> > x_dot_dot_b;
400  if (state->getXDotDot() != Teuchos::null) {
401  RCP<DPV> x_dot_dot = defaultProductVector(prod_space);
402  assign(x_dot_dot->getNonconstVectorBlock(0).ptr(),*(state->getXDotDot()));
403  assign(x_dot_dot->getNonconstVectorBlock(1).ptr(), zero);
404  x_dot_dot_b = x_dot_dot;
405  }
406 
407  RCP<SolutionState<Scalar> > prod_state =
408  rcp(new SolutionState<Scalar>(state->getMetaData()->clone(),
409  x_b, x_dot_b, x_dot_dot_b,
410  state->getStepperState()->clone()));
411  solutionHistory_->addState(prod_state);
412  }
413 
414  RCP<const VectorBase<Scalar> > frozen_x =
415  state_solution_history->getCurrentState()->getX();
416  RCP<const VectorBase<Scalar> > frozen_x_dot =
417  state_solution_history->getCurrentState()->getXDot();
418  RCP<const VectorBase<Scalar> > frozen_x_dot_dot =
419  state_solution_history->getCurrentState()->getXDotDot();
420  RCP<const SolutionHistory<Scalar> > sens_solution_history =
421  sens_integrator_->getSolutionHistory();
422  num_states = sens_solution_history->getNumStates();
423  for (int i=0; i<num_states; ++i) {
424  RCP<const SolutionState<Scalar> > state = (*sens_solution_history)[i];
425 
426  // X
427  RCP<DPV> x = defaultProductVector(prod_space);
428  assign(x->getNonconstVectorBlock(0).ptr(), *frozen_x);
429  assign(x->getNonconstVectorBlock(1).ptr(), *(state->getX()));
430  RCP<VectorBase<Scalar> > x_b = x;
431 
432  // X-Dot
433  RCP<VectorBase<Scalar> > x_dot_b;
434  if (state->getXDot() != Teuchos::null) {
435  RCP<DPV> x_dot = defaultProductVector(prod_space);
436  assign(x_dot->getNonconstVectorBlock(0).ptr(), *frozen_x_dot);
437  assign(x_dot->getNonconstVectorBlock(1).ptr(), *(state->getXDot()));
438  x_dot_b = x_dot;
439  }
440 
441  // X-Dot-Dot
442  RCP<VectorBase<Scalar> > x_dot_dot_b;
443  if (state->getXDotDot() != Teuchos::null) {
444  RCP<DPV> x_dot_dot = defaultProductVector(prod_space);
445  assign(x_dot_dot->getNonconstVectorBlock(0).ptr(), *frozen_x_dot_dot);
446  assign(x_dot_dot->getNonconstVectorBlock(1).ptr(),*(state->getXDotDot()));
447  x_dot_dot_b = x_dot_dot;
448  }
449 
450  RCP<SolutionState<Scalar> > prod_state =
451  rcp(new SolutionState<Scalar>(state->getMetaData()->clone(),
452  x_b, x_dot_b, x_dot_dot_b,
453  state->getStepperState()->clone()));
454  solutionHistory_->addState(prod_state);
455  }
456 }
457 
458 /// Non-member constructor
459 template<class Scalar>
460 Teuchos::RCP<Tempus::IntegratorPseudoTransientAdjointSensitivity<Scalar> >
462  Teuchos::RCP<Teuchos::ParameterList> pList,
463  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model)
464 {
465  Teuchos::RCP<Tempus::IntegratorPseudoTransientAdjointSensitivity<Scalar> > integrator =
466  Teuchos::rcp(new Tempus::IntegratorPseudoTransientAdjointSensitivity<Scalar>(pList, model));
467  return(integrator);
468 }
469 
470 /// Non-member constructor
471 template<class Scalar>
472 Teuchos::RCP<Tempus::IntegratorPseudoTransientAdjointSensitivity<Scalar> >
474  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model,
475  std::string stepperType)
476 {
477  Teuchos::RCP<Tempus::IntegratorPseudoTransientAdjointSensitivity<Scalar> > integrator =
478  Teuchos::rcp(new Tempus::IntegratorPseudoTransientAdjointSensitivity<Scalar>(model, stepperType));
479  return(integrator);
480 }
481 
482 /// Non-member constructor
483 template<class Scalar>
484 Teuchos::RCP<Tempus::IntegratorPseudoTransientAdjointSensitivity<Scalar> >
486 {
487  Teuchos::RCP<Tempus::IntegratorPseudoTransientAdjointSensitivity<Scalar> > integrator =
489  return(integrator);
490 }
491 
492 } // namespace Tempus
493 #endif // Tempus_IntegratorPseudoTransientAdjointSensitivity_impl_hpp
static Teuchos::RCP< const Teuchos::ParameterList > getValidParameters()
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &pl) override
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const override
virtual void setTempusParameterList(Teuchos::RCP< Teuchos::ParameterList > pl) override
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const override
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getDgDp() const
Return adjoint sensitivity stored in gradient format.
Time integrator suitable for pseudotransient adjoint sensitivity analysis.
Teuchos::RCP< Tempus::IntegratorPseudoTransientAdjointSensitivity< Scalar > > integratorPseudoTransientAdjointSensitivity(Teuchos::RCP< Teuchos::ParameterList > pList, const Teuchos::RCP< Thyra::ModelEvaluator< Scalar > > &model)
Non-member constructor.
virtual Teuchos::RCP< TimeStepControl< Scalar > > getNonConstTimeStepControl() override
Status
Status for the Integrator, the Stepper and the SolutionState.
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getXdot() const
Get current the time derivative of the solution, xdot.
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 > > y0=Teuchos::null, Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > ydot0=Teuchos::null, Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > ydotdot0=Teuchos::null)
Set the initial state from Thyra::VectorBase(s)
Teuchos::RCP< AdjointSensitivityModelEvaluator< Scalar > > createSensitivityModel(const Teuchos::RCP< Thyra::ModelEvaluator< Scalar > > &model, const Teuchos::RCP< Teuchos::ParameterList > &inputPL)
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getXdotdot() const
Get current the second time derivative of the solution, xdotdot.
virtual Teuchos::RCP< Teuchos::ParameterList > getTempusParameterList() override
Return a copy of the Tempus ParameterList.
virtual Teuchos::RCP< Stepper< Scalar > > getStepper() const override
Get the Stepper.
virtual Teuchos::RCP< const SolutionHistory< Scalar > > getSolutionHistory() const override
Get the SolutionHistory.
virtual Teuchos::RCP< const TimeStepControl< Scalar > > getTimeStepControl() const override
Get the TimeStepControl.
ModelEvaluator for forming adjoint sensitivity equations.
virtual bool advanceTime()
Advance the solution to timeMax, and return true if successful.
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getX() const
Get current the solution, x.
Solution state for integrators and steppers. SolutionState contains the metadata for solutions and th...