Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Tempus_IntegratorPseudoTransientAdjointSensitivity_decl.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_decl_hpp
10 #define Tempus_IntegratorPseudoTransientAdjointSensitivity_decl_hpp
11 
12 // Tempus
13 #include "Tempus_IntegratorBasic.hpp"
14 #include "Tempus_AdjointSensitivityModelEvaluator.hpp"
15 
16 namespace Tempus {
17 
18 
19 /** \brief Time integrator suitable for pseudotransient adjoint sensitivity
20  * analysis */
21 /**
22  * For some problems, time integrators are used to compute steady-state
23  * solutions (also known as pseudo-transient solvers). When computing
24  * sensitivities, it is not necessary in these cases to propagate sensitivities
25  * all the way through the forward time integration. Instead the steady-state
26  * is first computed as usual, and then the sensitivities are computed using
27  * a similar pseudo-transient time integration applied to the adjoint
28  * sensitivity equations with the state frozen to the computed steady-state.
29  * This integrator specializes the transient sensitivity methods implemented by
30  * Tempus::IntegratorAdjointSensitivity to this case.
31  *
32  * Consider an implicit ODE f(x_dot,x,p) = 0 with a stable steady-state solution
33  * x = x^s, x_dot = 0 where f(0,x^s,p) = 0 and all of the eigenvalues of
34  * df/dx(0,x^s,p) are in the right half-plane (for an explicit ODE, the
35  * eigenvalues must be in the left half-plane). In the pseudo-transient method
36  * a time-integrator is applied to f(x_dot,x,p) = 0 until x_dot is sufficiently
37  * small. Now consider the adjoint sensitivity equations for some response
38  * function g(x,p):
39  * df/dx_dot^T*y_dot + df/dx^T*y - dg/dx^T = 0
40  * after the transformation tau = T - t has been applied, where T is the final
41  * time. For pseudo-transient adjoint sensitivities, the above is integrated
42  * from y(0) = 0 until y_dot is sufficiently small, in which case
43  * y^s = (df/dx)^{-T}*(dg/dx)^T.
44  * Then the final sensitivity of g is
45  * dg/dp^T - df/dp^T*y^s.
46  * One can see that y^s is the only steady-state solution of the adjoint
47  * equations, since df/dx and dg/dx are constant, and must be linearly stable
48  * (since the eigenvalues of df/dx^T are the same as df/dx).
49  */
50 template<class Scalar>
52  virtual public Tempus::Integrator<Scalar>
53 {
54 public:
55 
56  /** \brief Constructor with ParameterList and model, and will be fully
57  * initialized. */
58  /*!
59  * In addition to all of the regular integrator options, the supplied
60  * parameter list supports the following options contained within a sublist
61  * "Sensitivities" from the top-level parameter list:
62  * <ul>
63  * <li> "Sensitivity Parameter Index", (default: 0) The model evaluator
64  * parameter index for which sensitivities will be computed.
65  * <li> "Response Function Index", (default: 0) The model evaluator
66  * response index for which sensitivities will be computed.
67  * <li> "Mass Matrix Is Constant" (default: true) Whether the mass matrix
68  * df/dx_dot is a constant matrix. As describe above, this is
69  * currently required to be true.
70  * <li> "Mass Matrix Is Identity" (default: false) Whether the mass matrix
71  * is the identity matrix, in which some computations can be skipped.
72  * </ul>
73  */
75  Teuchos::RCP<Teuchos::ParameterList> pList,
76  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model);
77 
78  /** \brief Constructor with model and "Stepper Type" and is fully initialized with default settings. */
80  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model,
81  std::string stepperType);
82 
83  /// Destructor
84  /** \brief Constructor that requires a subsequent setParameterList, setStepper, and initialize calls. */
86 
87  /// Destructor
89 
90  /// \name Basic integrator methods
91  //@{
92 
93  /// Advance the solution to timeMax, and return true if successful.
94  virtual bool advanceTime();
95  /// Advance the solution to timeFinal, and return true if successful.
96  virtual bool advanceTime(const Scalar timeFinal) override;
97  /// Get current time
98  virtual Scalar getTime() const override;
99  /// Get current index
100  virtual Scalar getIndex() const override;
101  /// Get Status
102  virtual Status getStatus() const override;
103  /// Get the Stepper
104  virtual Teuchos::RCP<Stepper<Scalar> > getStepper() const override;
105  /// Return a copy of the Tempus ParameterList
106  virtual Teuchos::RCP<Teuchos::ParameterList> getTempusParameterList() override;
107  virtual void setTempusParameterList(Teuchos::RCP<Teuchos::ParameterList> pl) override;
108  /// Get the SolutionHistory
109  virtual Teuchos::RCP<const SolutionHistory<Scalar> > getSolutionHistory() const override;
110  /// Get the TimeStepControl
111  virtual Teuchos::RCP<const TimeStepControl<Scalar> > getTimeStepControl() const override;
112  /// Returns the IntegratorTimer_ for this Integrator
113  virtual Teuchos::RCP<Teuchos::Time> getIntegratorTimer() const override
114  {return state_integrator_->getIntegratorTimer();}
115  virtual Teuchos::RCP<Teuchos::Time> getStepperTimer() const override
116  {return state_integrator_->getStepperTimer();}
117 
118  //@}
119 
120  /// Set the initial state from Thyra::VectorBase(s)
121  virtual void initializeSolutionHistory(
122  Scalar t0,
123  Teuchos::RCP<const Thyra::VectorBase<Scalar> > x0,
124  Teuchos::RCP<const Thyra::VectorBase<Scalar> > xdot0 = Teuchos::null,
125  Teuchos::RCP<const Thyra::VectorBase<Scalar> > xdotdot0 = Teuchos::null,
126  Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > y0 = Teuchos::null,
127  Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > ydot0 = Teuchos::null,
128  Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > ydotdot0 = Teuchos::null);
129 
130  /// Get current the solution, x
131  virtual Teuchos::RCP<const Thyra::VectorBase<Scalar> > getX() const;
132  /// Get current the time derivative of the solution, xdot
133  virtual Teuchos::RCP<const Thyra::VectorBase<Scalar> > getXdot() const;
134  /// Get current the second time derivative of the solution, xdotdot
135  virtual Teuchos::RCP<const Thyra::VectorBase<Scalar> > getXdotdot() const;
136 
137  /// Return adjoint sensitivity stored in gradient format
138  virtual Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > getDgDp() const;
139 
140  /// \name Overridden from Teuchos::ParameterListAcceptor
141  //@{
142  void setParameterList(const Teuchos::RCP<Teuchos::ParameterList> & pl)
143  override;
144  Teuchos::RCP<Teuchos::ParameterList> getNonconstParameterList() override;
145  Teuchos::RCP<Teuchos::ParameterList> unsetParameterList() override;
146 
147  Teuchos::RCP<const Teuchos::ParameterList> getValidParameters()
148  const override;
149  //@}
150 
151  /// \name Overridden from Teuchos::Describable
152  //@{
153  std::string description() const override;
154  void describe(Teuchos::FancyOStream & out,
155  const Teuchos::EVerbosityLevel verbLevel) const override;
156  //@}
157 
158 protected:
159  typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
160 
161  // Create sensitivity model evaluator from application model
162  Teuchos::RCP<AdjointSensitivityModelEvaluator<Scalar> >
164  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model,
165  const Teuchos::RCP<Teuchos::ParameterList>& inputPL);
166 
167  void buildSolutionHistory();
168 
169  Teuchos::RCP<Thyra::ModelEvaluator<Scalar> > model_;
170  Teuchos::RCP<AdjointSensitivityModelEvaluator<Scalar> > sens_model_;
171  Teuchos::RCP<IntegratorBasic<Scalar> > state_integrator_;
172  Teuchos::RCP<IntegratorBasic<Scalar> > sens_integrator_;
173  Teuchos::RCP<SolutionHistory<Scalar> > solutionHistory_;
174  Teuchos::RCP<DMVPV> dgdp_;
175 };
176 
177 /// Non-member constructor
178 template<class Scalar>
179 Teuchos::RCP<Tempus::IntegratorPseudoTransientAdjointSensitivity<Scalar> >
181  Teuchos::RCP<Teuchos::ParameterList> pList,
182  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model);
183 
184 /// Non-member constructor
185 template<class Scalar>
186 Teuchos::RCP<Tempus::IntegratorPseudoTransientAdjointSensitivity<Scalar> >
188  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model,
189  std::string stepperType);
190 
191 /// Non-member constructor
192 template<class Scalar>
193 Teuchos::RCP<Tempus::IntegratorPseudoTransientAdjointSensitivity<Scalar> >
195 
196 } // namespace Tempus
197 
198 #endif // Tempus_IntegratorPseudoTransientAdjointSensitivity_decl_hpp
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &pl) override
virtual Teuchos::RCP< Teuchos::Time > getIntegratorTimer() const override
Returns the IntegratorTimer_ for this Integrator.
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.
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)
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.
Thyra Base interface for time integrators. Time integrators are designed to advance the solution from...
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.
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.