Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups 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 "Thyra_DefaultMultiVectorProductVector.hpp"
13 #include "Thyra_VectorStdOps.hpp"
14 #include "Thyra_MultiVectorStdOps.hpp"
15 #include "Thyra_ImplicitAdjointModelEvaluator.hpp"
16 
17 
18 namespace Tempus {
19 
20 template<class Scalar>
21 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
22 IntegratorPseudoTransientAdjointSensitivity(
23  Teuchos::RCP<Teuchos::ParameterList> inputPL,
24  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model,
25  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& adjoint_residual_model,
26  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& adjoint_solve_model)
27 {
28  model_ = model;
29  adjoint_residual_model_ = adjoint_residual_model;
30  adjoint_solve_model_ = adjoint_solve_model;
31  state_integrator_ = createIntegratorBasic<Scalar>(inputPL, model_);
32  sens_model_ = createSensitivityModel(model_, adjoint_residual_model_,
33  adjoint_solve_model_, inputPL);
34  sens_integrator_ = createIntegratorBasic<Scalar>(inputPL, sens_model_);
35  stepMode_ = SensitivityStepMode::Forward;
36  do_forward_integration_ = true;
37  do_adjoint_integration_ = true;
38 }
39 
40 template<class Scalar>
41 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
42 IntegratorPseudoTransientAdjointSensitivity(
43  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model,
44  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& adjoint_residual_model,
45  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& adjoint_solve_model,
46  std::string stepperType)
47 {
48  model_ = model;
49  adjoint_residual_model_ = adjoint_residual_model;
50  adjoint_solve_model_ = adjoint_solve_model;
51  state_integrator_ = createIntegratorBasic<Scalar>(model_, stepperType);
52  sens_model_ = createSensitivityModel(model_, adjoint_residual_model_,
53  adjoint_solve_model_, Teuchos::null);
54  sens_integrator_ = createIntegratorBasic<Scalar>(sens_model_, stepperType);
55  stepMode_ = SensitivityStepMode::Forward;
56  do_forward_integration_ = true;
57  do_adjoint_integration_ = true;
58 }
59 
60 template<class Scalar>
61 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
62 IntegratorPseudoTransientAdjointSensitivity(
63  Teuchos::RCP<Teuchos::ParameterList> inputPL,
64  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model,
65  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& adjoint_model) :
66  IntegratorPseudoTransientAdjointSensitivity(
67  inputPL, model, adjoint_model, adjoint_model)
68 {
69 }
70 
71 template<class Scalar>
72 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
73 IntegratorPseudoTransientAdjointSensitivity(
74  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model,
75  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& adjoint_model,
76  std::string stepperType) :
77  IntegratorPseudoTransientAdjointSensitivity(
78  model, adjoint_model, adjoint_model, stepperType)
79 {
80 }
81 
82 template<class Scalar>
83 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
84 IntegratorPseudoTransientAdjointSensitivity(
85  Teuchos::RCP<Teuchos::ParameterList> inputPL,
86  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model) :
87  IntegratorPseudoTransientAdjointSensitivity(
88  inputPL, model, Thyra::implicitAdjointModelEvaluator(model))
89 {
90 }
91 
92 template<class Scalar>
93 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
94 IntegratorPseudoTransientAdjointSensitivity(
95  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model,
96  std::string stepperType) :
97  IntegratorPseudoTransientAdjointSensitivity(
98  model, Thyra::implicitAdjointModelEvaluator(model), stepperType)
99 {
100 }
101 
102 template<class Scalar>
103 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
104 IntegratorPseudoTransientAdjointSensitivity()
105 {
106  state_integrator_ = createIntegratorBasic<Scalar>();
107  sens_integrator_ = createIntegratorBasic<Scalar>();
108  stepMode_ = SensitivityStepMode::Forward;
109 }
110 
111 template<class Scalar>
112 bool
113 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
114 advanceTime()
115 {
116  const Scalar tfinal =
117  state_integrator_->getTimeStepControl()->getFinalTime();
118  return advanceTime(tfinal);
119 }
120 
121 template<class Scalar>
122 bool
123 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
124 advanceTime(const Scalar timeFinal)
125 {
126  TEMPUS_FUNC_TIME_MONITOR_DIFF("Tempus::IntegratorPseudoTransientAdjointSensitivity::advanceTime()", TEMPUS_PTAS_AT);
127 
128  using Teuchos::RCP;
129  using Thyra::VectorBase;
130  typedef Thyra::ModelEvaluatorBase MEB;
131 
132  bool state_status = true;
133  if (do_forward_integration_) {
134  TEMPUS_FUNC_TIME_MONITOR_DIFF("Tempus::IntegratorPseudoTransientAdjointSensitivity::advanceTime::state", TEMPUS_PTAS_AT_FWD);
135 
136  // Run state integrator and get solution
137  stepMode_ = SensitivityStepMode::Forward;
138  state_status = state_integrator_->advanceTime(timeFinal);
139  }
140 
141  bool sens_status = true;
142  if (do_adjoint_integration_) {
143  TEMPUS_FUNC_TIME_MONITOR_DIFF("Tempus::IntegratorPseudoTransientAdjointSensitivity::advanceTime::adjoint", TEMPUS_PTAS_AT_ADJ);
144 
145  // For at least some time-stepping methods, the time of the last time step
146  // may not be timeFinal (e.g., it may be greater by at most delta_t).
147  // But since the adjoint model requires timeFinal in its formulation, reset
148  // it to the achieved final time.
149  sens_model_->setFinalTime(state_integrator_->getTime());
150 
151  // Set solution in sensitivity ME
152  sens_model_->setForwardSolutionHistory(
153  state_integrator_->getSolutionHistory());
154 
155  // Run sensitivity integrator
156  stepMode_ = SensitivityStepMode::Adjoint;
157  sens_status = sens_integrator_->advanceTime(timeFinal);
158 
159  // Compute final dg/dp, g which is computed by response 0, 1 of the adjoint
160  // model evaluator
161  MEB::InArgs<Scalar> inargs = sens_model_->getNominalValues();
162  MEB::OutArgs<Scalar> outargs = sens_model_->createOutArgs();
163  inargs.set_t(sens_integrator_->getTime());
164  inargs.set_x(sens_integrator_->getX());
165  if (inargs.supports(MEB::IN_ARG_x_dot))
166  inargs.set_x_dot(sens_integrator_->getXDot());
167  if (inargs.supports(MEB::IN_ARG_x_dot_dot))
168  inargs.set_x_dot_dot(sens_integrator_->getXDotDot());
169  RCP<VectorBase<Scalar> > G = dgdp_;
170  if (G == Teuchos::null) {
171  G = Thyra::createMember(sens_model_->get_g_space(0));
172  dgdp_ = Teuchos::rcp_dynamic_cast<DMVPV>(G);
173  }
174  if (g_ == Teuchos::null)
175  g_ = Thyra::createMember(sens_model_->get_g_space(1));
176  outargs.set_g(0, G);
177  outargs.set_g(1, g_);
178  sens_model_->evalModel(inargs, outargs);
179 
180  buildSolutionHistory();
181  }
182 
183  return state_status && sens_status;
184 }
185 
186 template<class Scalar>
187 Scalar
188 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
189 getTime() const
190 {
191  return solutionHistory_->getCurrentTime();
192 }
193 
194 template<class Scalar>
195 int
196 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
197 getIndex() const
198 {
199  return solutionHistory_->getCurrentIndex();
200 }
201 
202 template<class Scalar>
203 Status
204 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
205 getStatus() const
206 {
207  Status state_status = state_integrator_->getStatus();
208  Status sens_status = sens_integrator_->getStatus();
209  if (state_status == FAILED || sens_status == FAILED)
210  return FAILED;
211  if (state_status == WORKING || sens_status == WORKING)
212  return WORKING;
213  return PASSED;
214 }
215 
216 template <class Scalar>
217 void IntegratorPseudoTransientAdjointSensitivity<Scalar>::setStatus(
218  const Status st) {
219  state_integrator_->setStatus(st);
220  sens_integrator_->setStatus(st);
221 }
222 
223 template<class Scalar>
224 Teuchos::RCP<Stepper<Scalar> >
225 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
226 getStepper() const
227 {
228  return state_integrator_->getStepper();
229 }
230 
231 template<class Scalar>
232 Teuchos::RCP<Stepper<Scalar> >
233 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
234 getStateStepper() const
235 {
236  return state_integrator_->getStepper();
237 }
238 
239 template<class Scalar>
240 Teuchos::RCP<Stepper<Scalar> >
241 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
242 getSensStepper() const
243 {
244  return sens_integrator_->getStepper();
245 }
246 
247 template<class Scalar>
248 Teuchos::RCP<const SolutionHistory<Scalar> >
249 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
250 getSolutionHistory() const
251 {
252  return solutionHistory_;
253 }
254 
255 template<class Scalar>
256 Teuchos::RCP<const SolutionHistory<Scalar> >
257 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
258 getStateSolutionHistory() const
259 {
260  return state_integrator_->getSolutionHistory();
261 }
262 
263 template<class Scalar>
264 Teuchos::RCP<const SolutionHistory<Scalar> >
265 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
266 getSensSolutionHistory() const
267 {
268  return sens_integrator_->getSolutionHistory();
269 }
270 
271 template<class Scalar>
272 Teuchos::RCP<SolutionHistory<Scalar> >
273 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
274 getNonConstSolutionHistory()
275 {
276  return solutionHistory_;
277 }
278 
279 template<class Scalar>
280 Teuchos::RCP<const TimeStepControl<Scalar> >
281 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
282 getTimeStepControl() const
283 {
284  return state_integrator_->getTimeStepControl();
285 }
286 
287 template<class Scalar>
288 Teuchos::RCP<TimeStepControl<Scalar> >
289 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
290 getNonConstTimeStepControl()
291 {
292  return state_integrator_->getNonConstTimeStepControl();
293 }
294 
295 template<class Scalar>
296 Teuchos::RCP<TimeStepControl<Scalar> >
297 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
298 getStateNonConstTimeStepControl()
299 {
300  return state_integrator_->getNonConstTimeStepControl();
301 }
302 
303 template<class Scalar>
304 Teuchos::RCP<TimeStepControl<Scalar> >
305 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
306 getSensNonConstTimeStepControl()
307 {
308  return sens_integrator_->getNonConstTimeStepControl();
309 }
310 
311 template<class Scalar>
312 Teuchos::RCP<IntegratorObserver<Scalar> >
313 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
314 getObserver()
315 {
316  return state_integrator_->getObserver();
317 }
318 
319 template<class Scalar>
320 void
321 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
322 setObserver(Teuchos::RCP<IntegratorObserver<Scalar> > obs)
323 {
324  state_integrator_->setObserver(obs);
325  sens_integrator_->setObserver(obs);
326 }
327 
328 template<class Scalar>
329 void IntegratorPseudoTransientAdjointSensitivity<Scalar>::
330 initializeSolutionHistory(Scalar t0,
331  Teuchos::RCP<const Thyra::VectorBase<Scalar> > x0,
332  Teuchos::RCP<const Thyra::VectorBase<Scalar> > xdot0,
333  Teuchos::RCP<const Thyra::VectorBase<Scalar> > xdotdot0,
334  Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > y0,
335  Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > ydot0,
336  Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > ydotdot0)
337 {
338  using Teuchos::RCP;
339  using Teuchos::rcp_dynamic_cast;
340  using Thyra::VectorSpaceBase;
341  using Thyra::assign;
342  using Thyra::createMember;
343 
344  //
345  // Create and initialize product X, Xdot, Xdotdot
346 
347  RCP< const VectorSpaceBase<Scalar> > space = sens_model_->get_x_space();
348  RCP<DMVPV> Y = rcp_dynamic_cast<DMVPV>(createMember(space));
349  RCP<DMVPV> Ydot = rcp_dynamic_cast<DMVPV>(createMember(space));
350  RCP<DMVPV> Ydotdot = rcp_dynamic_cast<DMVPV>(createMember(space));
351  const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
352 
353  // y
354  if (y0 == Teuchos::null)
355  assign(Y->getNonconstMultiVector().ptr(), zero);
356  else
357  assign(Y->getNonconstMultiVector().ptr(), *y0);
358 
359  // ydot
360  if (ydot0 == Teuchos::null)
361  assign(Ydot->getNonconstMultiVector().ptr(), zero);
362  else
363  assign(Ydot->getNonconstMultiVector().ptr(), *ydot0);
364 
365  // ydotdot
366  if (ydotdot0 == Teuchos::null)
367  assign(Ydotdot->getNonconstMultiVector().ptr(), zero);
368  else
369  assign(Ydotdot->getNonconstMultiVector().ptr(), *ydotdot0);
370 
371  state_integrator_->initializeSolutionHistory(t0, x0, xdot0, xdotdot0);
372  sens_integrator_->initializeSolutionHistory(t0, Y, Ydot, Ydotdot);
373 }
374 
375 template<class Scalar>
376 Teuchos::RCP<const Thyra::VectorBase<Scalar> >
377 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
378 getX() const
379 {
380  return state_integrator_->getX();
381 }
382 
383 template<class Scalar>
384 Teuchos::RCP<const Thyra::VectorBase<Scalar> >
385 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
386 getXDot() const
387 {
388  return state_integrator_->getXDot();
389 }
390 
391 template<class Scalar>
392 Teuchos::RCP<const Thyra::VectorBase<Scalar> >
393 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
394 getXDotDot() const
395 {
396  return state_integrator_->getXDotDot();
397 }
398 
399 template<class Scalar>
400 Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
401 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
402 getY() const
403 {
404  using Teuchos::RCP;
405  using Teuchos::rcp_dynamic_cast;
406  RCP<const DMVPV> mvpv =
407  rcp_dynamic_cast<const DMVPV>(sens_integrator_->getX());
408  return mvpv->getMultiVector();
409 }
410 
411 template<class Scalar>
412 Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
413 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
414 getYDot() const
415 {
416  using Teuchos::RCP;
417  using Teuchos::rcp_dynamic_cast;
418  RCP<const DMVPV> mvpv =
419  rcp_dynamic_cast<const DMVPV>(sens_integrator_->getXDot());
420  return mvpv->getMultiVector();
421 }
422 
423 template<class Scalar>
424 Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
425 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
426 getYDotDot() const
427 {
428  using Teuchos::RCP;
429  using Teuchos::rcp_dynamic_cast;
430  RCP<const DMVPV> mvpv =
431  rcp_dynamic_cast<const DMVPV>(sens_integrator_->getXDotDot());
432  return mvpv->getMultiVector();
433 }
434 
435 template<class Scalar>
436 Teuchos::RCP<const Thyra::VectorBase<Scalar> >
437 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
438 getG() const
439 {
440  return g_;
441 }
442 
443 template<class Scalar>
444 Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
445 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
446 getDgDp() const
447 {
448  return dgdp_->getMultiVector();
449 }
450 
451 template<class Scalar>
452 std::string
453 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
454 description() const
455 {
456  std::string name = "Tempus::IntegratorPseudoTransientAdjointSensitivity";
457  return(name);
458 }
459 
460 template<class Scalar>
461 void
462 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
463 describe(
464  Teuchos::FancyOStream &out,
465  const Teuchos::EVerbosityLevel verbLevel) const
466 {
467  auto l_out = Teuchos::fancyOStream( out.getOStream() );
468  Teuchos::OSTab ostab(*l_out, 2, this->description());
469  l_out->setOutputToRootOnly(0);
470 
471  *l_out << description() << "::describe" << std::endl;
472  state_integrator_->describe(*l_out, verbLevel);
473  sens_integrator_->describe(*l_out, verbLevel);
474 }
475 
476 template<class Scalar>
477 void
478 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
479 setParameterList(const Teuchos::RCP<Teuchos::ParameterList> & inputPL)
480 {
481  // IntegratorBasic is no longer a Teuchos::ParameterListAcceptor.
482  // Since setting the ParameterList is essentially a complete reset,
483  // we will rebuild from scratch and reuse the ModelEvaluator.
484  auto model = Teuchos::rcp_const_cast<Thyra::ModelEvaluator<Scalar>> (
485  state_integrator_->getStepper()->getModel());
486  auto tmp_state_integrator = createIntegratorBasic<Scalar>(inputPL, model);
487  state_integrator_->copy(tmp_state_integrator);
488 
489  model = Teuchos::rcp_const_cast<Thyra::ModelEvaluator<Scalar>> (
490  sens_integrator_->getStepper()->getModel());
491  auto tmp_sens_integrator = createIntegratorBasic<Scalar>(inputPL, model);
492  sens_integrator_->copy(tmp_sens_integrator);
493 }
494 
495 template<class Scalar>
496 Teuchos::RCP<Teuchos::ParameterList>
497 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
498 unsetParameterList()
499 {
500  // IntegratorBasic is no longer a Teuchos::ParameterListAcceptor.
501  // We will treat unsetting the ParameterList as a reset to default
502  // settings, and reuse the ModelEvaluator.
503  auto tmp_state_integrator = createIntegratorBasic<Scalar>();
504  auto model = state_integrator_->getStepper()->getModel();
505  tmp_state_integrator->setModel(model);
506  state_integrator_->copy(tmp_state_integrator);
507 
508  auto tmp_sens_integrator = createIntegratorBasic<Scalar>();
509  model = sens_integrator_->getStepper()->getModel();
510  tmp_sens_integrator->setModel(model);
511  sens_integrator_->copy(tmp_sens_integrator);
512 
513  auto pl = Teuchos::rcp_const_cast<Teuchos::ParameterList> (
514  sens_integrator_->getValidParameters());
515  return pl;
516 }
517 
518 template<class Scalar>
519 Teuchos::RCP<const Teuchos::ParameterList>
520 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
521 getValidParameters() const
522 {
523  Teuchos::RCP<Teuchos::ParameterList> pl =
524  Teuchos::rcp(new Teuchos::ParameterList);
525  Teuchos::RCP<const Teuchos::ParameterList> integrator_pl =
526  state_integrator_->getValidParameters();
527  Teuchos::RCP<const Teuchos::ParameterList> sensitivity_pl =
528  AdjointSensitivityModelEvaluator<Scalar>::getValidParameters();
529  pl->setParameters(*integrator_pl);
530  pl->sublist("Sensitivities").setParameters(*sensitivity_pl);
531 
532  return pl;
533 }
534 
535 template<class Scalar>
536 Teuchos::RCP<Teuchos::ParameterList>
537 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
538 getNonconstParameterList()
539 {
540  auto pl = Teuchos::rcp_const_cast<Teuchos::ParameterList> (
541  state_integrator_->getValidParameters());
542  return pl;
543 }
544 
545 template<class Scalar>
546 SensitivityStepMode
547 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
548 getStepMode() const
549 {
550  return stepMode_;
551 }
552 
553 template <class Scalar>
554 Teuchos::RCP<AdjointSensitivityModelEvaluator<Scalar> >
555 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
556 createSensitivityModel(
557  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model,
558  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& adjoint_residual_model,
559  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& adjoint_solve_model,
560  const Teuchos::RCP<Teuchos::ParameterList>& inputPL)
561 {
562  using Teuchos::rcp;
563 
564  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
565  if (inputPL != Teuchos::null) {
566  *pl = inputPL->sublist("Sensitivities");
567  }
568  const Scalar tinit = state_integrator_->getTimeStepControl()->getInitTime();
569  const Scalar tfinal = state_integrator_->getTimeStepControl()->getFinalTime();
570  return rcp(new AdjointSensitivityModelEvaluator<Scalar>(
571  model, adjoint_residual_model, adjoint_solve_model_,
572  tinit, tfinal, true, pl));
573 }
574 
575 template<class Scalar>
576 void
577 IntegratorPseudoTransientAdjointSensitivity<Scalar>::
578 buildSolutionHistory()
579 {
580  using Teuchos::RCP;
581  using Teuchos::rcp;
582  using Teuchos::rcp_dynamic_cast;
583  using Teuchos::ParameterList;
584  using Thyra::VectorBase;
585  using Thyra::VectorSpaceBase;
586  using Thyra::assign;
587  using Thyra::defaultProductVector;
588  typedef Thyra::DefaultProductVector<Scalar> DPV;
589  typedef Thyra::DefaultProductVectorSpace<Scalar> DPVS;
590 
591  // Create combined solution histories, first for the states with zero
592  // adjoint and then for the adjoint with frozen states
593  auto shPL = Teuchos::rcp_const_cast<Teuchos::ParameterList> (
594  state_integrator_->getSolutionHistory()->getValidParameters());
595  solutionHistory_ = createSolutionHistoryPL<Scalar>(shPL);
596 
597  RCP<const VectorSpaceBase<Scalar> > x_space =
598  model_->get_x_space();
599  RCP<const VectorSpaceBase<Scalar> > adjoint_space =
600  sens_model_->get_x_space();
601  Teuchos::Array< RCP<const VectorSpaceBase<Scalar> > > spaces(2);
602  spaces[0] = x_space;
603  spaces[1] = adjoint_space;
604  RCP<const DPVS > prod_space = Thyra::productVectorSpace(spaces());
605  const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
606 
607  RCP<const SolutionHistory<Scalar> > state_solution_history =
608  state_integrator_->getSolutionHistory();
609  int num_states = state_solution_history->getNumStates();
610  for (int i=0; i<num_states; ++i) {
611  RCP<const SolutionState<Scalar> > state = (*state_solution_history)[i];
612 
613  // X
614  RCP<DPV> x = defaultProductVector(prod_space);
615  assign(x->getNonconstVectorBlock(0).ptr(), *(state->getX()));
616  assign(x->getNonconstVectorBlock(1).ptr(), zero);
617  RCP<VectorBase<Scalar> > x_b = x;
618 
619  // X-Dot
620  RCP<VectorBase<Scalar> > x_dot_b;
621  if (state->getXDot() != Teuchos::null) {
622  RCP<DPV> x_dot = defaultProductVector(prod_space);
623  assign(x_dot->getNonconstVectorBlock(0).ptr(), *(state->getXDot()));
624  assign(x_dot->getNonconstVectorBlock(1).ptr(), zero);
625  x_dot_b = x_dot;
626  }
627 
628  // X-Dot-Dot
629  RCP<VectorBase<Scalar> > x_dot_dot_b;
630  if (state->getXDotDot() != Teuchos::null) {
631  RCP<DPV> x_dot_dot = defaultProductVector(prod_space);
632  assign(x_dot_dot->getNonconstVectorBlock(0).ptr(),*(state->getXDotDot()));
633  assign(x_dot_dot->getNonconstVectorBlock(1).ptr(), zero);
634  x_dot_dot_b = x_dot_dot;
635  }
636 
637  RCP<SolutionState<Scalar> > prod_state = state->clone();
638  prod_state->setX(x_b);
639  prod_state->setXDot(x_dot_b);
640  prod_state->setXDotDot(x_dot_dot_b);
641  solutionHistory_->addState(prod_state);
642  }
643 
644  RCP<const VectorBase<Scalar> > frozen_x =
645  state_solution_history->getCurrentState()->getX();
646  RCP<const VectorBase<Scalar> > frozen_x_dot =
647  state_solution_history->getCurrentState()->getXDot();
648  RCP<const VectorBase<Scalar> > frozen_x_dot_dot =
649  state_solution_history->getCurrentState()->getXDotDot();
650  RCP<const SolutionHistory<Scalar> > sens_solution_history =
651  sens_integrator_->getSolutionHistory();
652  num_states = sens_solution_history->getNumStates();
653  for (int i=0; i<num_states; ++i) {
654  RCP<const SolutionState<Scalar> > state = (*sens_solution_history)[i];
655 
656  // X
657  RCP<DPV> x = defaultProductVector(prod_space);
658  assign(x->getNonconstVectorBlock(0).ptr(), *frozen_x);
659  assign(x->getNonconstVectorBlock(1).ptr(), *(state->getX()));
660  RCP<VectorBase<Scalar> > x_b = x;
661 
662  // X-Dot
663  RCP<VectorBase<Scalar> > x_dot_b;
664  if (state->getXDot() != Teuchos::null) {
665  RCP<DPV> x_dot = defaultProductVector(prod_space);
666  assign(x_dot->getNonconstVectorBlock(0).ptr(), *frozen_x_dot);
667  assign(x_dot->getNonconstVectorBlock(1).ptr(), *(state->getXDot()));
668  x_dot_b = x_dot;
669  }
670 
671  // X-Dot-Dot
672  RCP<VectorBase<Scalar> > x_dot_dot_b;
673  if (state->getXDotDot() != Teuchos::null) {
674  RCP<DPV> x_dot_dot = defaultProductVector(prod_space);
675  assign(x_dot_dot->getNonconstVectorBlock(0).ptr(), *frozen_x_dot_dot);
676  assign(x_dot_dot->getNonconstVectorBlock(1).ptr(),*(state->getXDotDot()));
677  x_dot_dot_b = x_dot_dot;
678  }
679 
680  RCP<SolutionState<Scalar> > prod_state = state->clone();
681  prod_state->setX(x_b);
682  prod_state->setXDot(x_dot_b);
683  prod_state->setXDotDot(x_dot_dot_b);
684  solutionHistory_->addState(prod_state, false);
685  }
686 }
687 
689 template<class Scalar>
690 Teuchos::RCP<IntegratorPseudoTransientAdjointSensitivity<Scalar> >
691 integratorPseudoTransientAdjointSensitivity(
692  Teuchos::RCP<Teuchos::ParameterList> pList,
693  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model)
694 {
695  Teuchos::RCP<IntegratorPseudoTransientAdjointSensitivity<Scalar> > integrator =
696  Teuchos::rcp(new IntegratorPseudoTransientAdjointSensitivity<Scalar>(pList, model));
697  return(integrator);
698 }
699 
701 template<class Scalar>
702 Teuchos::RCP<IntegratorPseudoTransientAdjointSensitivity<Scalar> >
703 integratorPseudoTransientAdjointSensitivity(
704  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model,
705  std::string stepperType)
706 {
707  Teuchos::RCP<IntegratorPseudoTransientAdjointSensitivity<Scalar> > integrator =
708  Teuchos::rcp(new IntegratorPseudoTransientAdjointSensitivity<Scalar>(model, stepperType));
709  return(integrator);
710 }
711 
713 template<class Scalar>
714 Teuchos::RCP<IntegratorPseudoTransientAdjointSensitivity<Scalar> >
715 integratorPseudoTransientAdjointSensitivity(
716  Teuchos::RCP<Teuchos::ParameterList> pList,
717  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model,
718  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& adjoint_model)
719 {
720  Teuchos::RCP<IntegratorPseudoTransientAdjointSensitivity<Scalar> > integrator =
721  Teuchos::rcp(new IntegratorPseudoTransientAdjointSensitivity<Scalar>(pList, model, adjoint_model));
722  return(integrator);
723 }
724 
726 template<class Scalar>
727 Teuchos::RCP<IntegratorPseudoTransientAdjointSensitivity<Scalar> >
728 integratorPseudoTransientAdjointSensitivity(
729  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model,
730  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& adjoint_model,
731  std::string stepperType)
732 {
733  Teuchos::RCP<IntegratorPseudoTransientAdjointSensitivity<Scalar> > integrator =
734  Teuchos::rcp(new IntegratorPseudoTransientAdjointSensitivity<Scalar>(model, adjoint_model, stepperType));
735  return(integrator);
736 }
737 
739 template<class Scalar>
740 Teuchos::RCP<IntegratorPseudoTransientAdjointSensitivity<Scalar> >
741 integratorPseudoTransientAdjointSensitivity(
742  Teuchos::RCP<Teuchos::ParameterList> pList,
743  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model,
744  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& adjoint_residual_model,
745  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& adjoint_solve_model)
746 {
747  Teuchos::RCP<IntegratorPseudoTransientAdjointSensitivity<Scalar> > integrator =
748  Teuchos::rcp(new IntegratorPseudoTransientAdjointSensitivity<Scalar>(pList, model, adjoint_residual_model, adjoint_solve_model));
749  return(integrator);
750 }
751 
753 template<class Scalar>
754 Teuchos::RCP<IntegratorPseudoTransientAdjointSensitivity<Scalar> >
755 integratorPseudoTransientAdjointSensitivity(
756  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model,
757  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& adjoint_residual_model,
758  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& adjoint_solve_model,
759  std::string stepperType)
760 {
761  Teuchos::RCP<IntegratorPseudoTransientAdjointSensitivity<Scalar> > integrator =
762  Teuchos::rcp(new IntegratorPseudoTransientAdjointSensitivity<Scalar>(model, adjoint_residual_model, adjoint_solve_model, stepperType));
763  return(integrator);
764 }
765 
767 template<class Scalar>
768 Teuchos::RCP<IntegratorPseudoTransientAdjointSensitivity<Scalar> >
769 integratorPseudoTransientAdjointSensitivity()
770 {
771  Teuchos::RCP<IntegratorPseudoTransientAdjointSensitivity<Scalar> > integrator =
772  Teuchos::rcp(new IntegratorPseudoTransientAdjointSensitivity<Scalar>());
773  return(integrator);
774 }
775 
776 } // namespace Tempus
777 #endif // Tempus_IntegratorPseudoTransientAdjointSensitivity_impl_hpp
Tempus::AdjointSensitivityModelEvaluator::getValidParameters
static Teuchos::RCP< const Teuchos::ParameterList > getValidParameters()
Definition: Tempus_AdjointSensitivityModelEvaluator_impl.hpp:555
Tempus::IntegratorPseudoTransientAdjointSensitivity::setParameterList
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &pl) override
Definition: Tempus_IntegratorPseudoTransientAdjointSensitivity_impl.hpp:479
Tempus::IntegratorPseudoTransientAdjointSensitivity::describe
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const override
Definition: Tempus_IntegratorPseudoTransientAdjointSensitivity_impl.hpp:463
Tempus::IntegratorPseudoTransientAdjointSensitivity::getObserver
virtual Teuchos::RCP< IntegratorObserver< Scalar > > getObserver()
Get the Observer.
Definition: Tempus_IntegratorPseudoTransientAdjointSensitivity_impl.hpp:314
Tempus::IntegratorPseudoTransientAdjointSensitivity::getY
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getY() const
Get the current adjoint solution, y.
Definition: Tempus_IntegratorPseudoTransientAdjointSensitivity_impl.hpp:402
Tempus::IntegratorPseudoTransientAdjointSensitivity::getXDot
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getXDot() const
Get the current time derivative of the solution, xdot.
Definition: Tempus_IntegratorPseudoTransientAdjointSensitivity_impl.hpp:386
Tempus::integratorPseudoTransientAdjointSensitivity
Teuchos::RCP< Tempus::IntegratorPseudoTransientAdjointSensitivity< Scalar > > integratorPseudoTransientAdjointSensitivity(Teuchos::RCP< Teuchos::ParameterList > pList, const Teuchos::RCP< Thyra::ModelEvaluator< Scalar > > &model)
Nonmember constructor.
Definition: Tempus_IntegratorPseudoTransientAdjointSensitivity_impl.hpp:691
Tempus::IntegratorPseudoTransientAdjointSensitivity::setObserver
virtual void setObserver(Teuchos::RCP< IntegratorObserver< Scalar > > obs=Teuchos::null)
Set the Observer.
Definition: Tempus_IntegratorPseudoTransientAdjointSensitivity_impl.hpp:322
Tempus::IntegratorPseudoTransientAdjointSensitivity::getNonConstSolutionHistory
virtual Teuchos::RCP< SolutionHistory< Scalar > > getNonConstSolutionHistory() override
Get the SolutionHistory.
Definition: Tempus_IntegratorPseudoTransientAdjointSensitivity_impl.hpp:274
Teuchos::EVerbosityLevel
EVerbosityLevel
Teuchos::basic_OSTab
Tempus::IntegratorPseudoTransientAdjointSensitivity::getValidParameters
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const override
Definition: Tempus_IntegratorPseudoTransientAdjointSensitivity_impl.hpp:521
Tempus::IntegratorPseudoTransientAdjointSensitivity::getDgDp
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getDgDp() const
Return adjoint sensitivity stored in gradient format.
Definition: Tempus_IntegratorPseudoTransientAdjointSensitivity_impl.hpp:446
Tempus::IntegratorPseudoTransientAdjointSensitivity::getYDot
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getYDot() const
Get the current time derivative of the adjoint solution, ydot.
Definition: Tempus_IntegratorPseudoTransientAdjointSensitivity_impl.hpp:414
Tempus::IntegratorPseudoTransientAdjointSensitivity
Time integrator suitable for pseudotransient adjoint sensitivity analysis.
Definition: Tempus_IntegratorPseudoTransientAdjointSensitivity_decl.hpp:60
Tempus::IntegratorPseudoTransientAdjointSensitivity::getStepMode
SensitivityStepMode getStepMode() const
What mode the current time integration step is in.
Definition: Tempus_IntegratorPseudoTransientAdjointSensitivity_impl.hpp:548
Tempus::IntegratorPseudoTransientAdjointSensitivity::getG
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getG() const
Return response function g.
Definition: Tempus_IntegratorPseudoTransientAdjointSensitivity_impl.hpp:438
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Tempus::IntegratorPseudoTransientAdjointSensitivity::getNonConstTimeStepControl
virtual Teuchos::RCP< TimeStepControl< Scalar > > getNonConstTimeStepControl() override
Definition: Tempus_IntegratorPseudoTransientAdjointSensitivity_impl.hpp:290
Tempus::Status
Status
Status for the Integrator, the Stepper and the SolutionState.
Definition: Tempus_Types.hpp:19
Tempus::IntegratorObserver
IntegratorObserver class for time integrators.
Definition: Tempus_IntegratorObserver.hpp:42
Tempus::IntegratorPseudoTransientAdjointSensitivity::getXDotDot
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getXDotDot() const
Get the current second time derivative of the solution, xdotdot.
Definition: Tempus_IntegratorPseudoTransientAdjointSensitivity_impl.hpp:394
Teuchos::basic_FancyOStream
Tempus::IntegratorPseudoTransientAdjointSensitivity::initializeSolutionHistory
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)
Definition: Tempus_IntegratorPseudoTransientAdjointSensitivity_impl.hpp:330
Thyra::DefaultMultiVectorProductVector::getMultiVector
RCP< const MultiVectorBase< Scalar > > getMultiVector() const
Tempus::IntegratorPseudoTransientAdjointSensitivity::getYDotDot
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getYDotDot() const
Get the current second time derivative of the adjoint solution, ydotdot.
Definition: Tempus_IntegratorPseudoTransientAdjointSensitivity_impl.hpp:426
Teuchos::ParameterList::setParameters
ParameterList & setParameters(const ParameterList &source)
Thyra::implicitAdjointModelEvaluator
RCP< ImplicitAdjointModelEvaluator< Scalar > > implicitAdjointModelEvaluator(const RCP< const ModelEvaluator< Scalar > > &model)
Definition: Thyra_ImplicitAdjointModelEvaluator.hpp:104
Thyra_ImplicitAdjointModelEvaluator.hpp
Teuchos::basic_FancyOStream::getOStream
RCP< std::basic_ostream< char_type, traits_type > > getOStream()
Tempus::IntegratorPseudoTransientAdjointSensitivity::getStepper
virtual Teuchos::RCP< Stepper< Scalar > > getStepper() const override
Get the Stepper.
Definition: Tempus_IntegratorPseudoTransientAdjointSensitivity_impl.hpp:226
Tempus::IntegratorPseudoTransientAdjointSensitivity::getSolutionHistory
virtual Teuchos::RCP< const SolutionHistory< Scalar > > getSolutionHistory() const override
Get the SolutionHistory.
Definition: Tempus_IntegratorPseudoTransientAdjointSensitivity_impl.hpp:250
Tempus::IntegratorPseudoTransientAdjointSensitivity::getTimeStepControl
virtual Teuchos::RCP< const TimeStepControl< Scalar > > getTimeStepControl() const override
Get the TimeStepControl.
Definition: Tempus_IntegratorPseudoTransientAdjointSensitivity_impl.hpp:282
Tempus::AdjointSensitivityModelEvaluator
ModelEvaluator for forming adjoint sensitivity equations.
Definition: Tempus_AdjointSensitivityModelEvaluator_decl.hpp:54
Teuchos::ParameterList::sublist
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
Tempus::IntegratorPseudoTransientAdjointSensitivity::advanceTime
virtual bool advanceTime()
Advance the solution to timeMax, and return true if successful.
Definition: Tempus_IntegratorPseudoTransientAdjointSensitivity_impl.hpp:114
Tempus::IntegratorPseudoTransientAdjointSensitivity::createSensitivityModel
Teuchos::RCP< AdjointSensitivityModelEvaluator< Scalar > > createSensitivityModel(const Teuchos::RCP< Thyra::ModelEvaluator< Scalar > > &model, const Teuchos::RCP< Thyra::ModelEvaluator< Scalar > > &adjoint_residual_model, const Teuchos::RCP< Thyra::ModelEvaluator< Scalar > > &adjoint_solve_model, const Teuchos::RCP< Teuchos::ParameterList > &inputPL)
Definition: Tempus_IntegratorPseudoTransientAdjointSensitivity_impl.hpp:556
Tempus::IntegratorPseudoTransientAdjointSensitivity::getX
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getX() const
Get the current solution, x.
Definition: Tempus_IntegratorPseudoTransientAdjointSensitivity_impl.hpp:378
Generated on Mon Sep 25 2023 07:58:56 for Tempus by   doxygen 1.8.5