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