Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Tempus_SolutionHistory_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_SolutionHistory_decl_hpp
10 #define Tempus_SolutionHistory_decl_hpp
11 
12 
13 #include "Teuchos_VerboseObject.hpp"
14 #include "Teuchos_Describable.hpp"
15 #include "Teuchos_ParameterListAcceptorDefaultBase.hpp"
16 
17 #include "Thyra_VectorStdOps.hpp"
18 
19 #include "Tempus_config.hpp"
20 #include "Tempus_SolutionState.hpp"
21 #include "Tempus_Interpolator.hpp"
22 
23 
24 namespace Tempus {
25 
27  STORAGE_TYPE_INVALID = 0, ///< Invalid storage type
28  STORAGE_TYPE_KEEP_NEWEST = 1, ///< Keep the single newest state
29  STORAGE_TYPE_UNDO = 2, ///< Keep the 2 newest states for undo
30  STORAGE_TYPE_STATIC = 3, ///< Keep a fix number of states
31  STORAGE_TYPE_UNLIMITED = 4, ///< Grow the history as needed
32 };
33 
34 
35 /** \brief SolutionHistory is basically a container of SolutionStates.
36  * SolutionHistory maintains a collection of SolutionStates for later
37  * retrival and reuse, such as checkpointing, restart, and undo
38  * operations.
39  *
40  * The actual storage of the SolutionStates may take several forms:
41  * - in memory
42  * - on disk
43  * - combination
44  * - use ATDM DataWareHouse
45  * but the interface should be unchanged and very similar to other
46  * containers.
47  *
48  * SolutionHistory can fill in SolutionStates between other SolutionStates
49  * by either
50  * - Integrating from one SolutionState to the desired time
51  * - This might include methods like Griewank's algorithm.
52  * - Interpolating between SolutionStates
53  * - Interpolated SolutionStates may not be suitable for adjoint
54  * solutions, restart, or undo operations (see SolutionState).
55  *
56  * There are two sets of indices associated with SolutionHistory:
57  * timestep index and the SolutionHistory index. Users are interested
58  * in the timestep index as it indicates the solution increment, e.g.,
59  * \f$[x^0, ... , x^{n-2}, x^{n-1}, x^{n}]\f$, where \f$n\f$ is the timestep
60  * index. While developers can also be interested in the SolutionHistory
61  * index to access the correct SolutionState in the SolutionHistory, e.g.,
62  * \f$[S^0, ... , S^{M-1}]\f$, where \f$M\f$ is the number of SolutionStates
63  * in the SolutionHistory (1 \f$\le M \le\f$ StorageLimit_).
64  *
65  * SolutionHistory will hold upto the StorageLimit_ number of states.
66  * Rules of thumb for the <b>minimum</b> storage limit:
67  * - Explicit one-step methods can update the solution state in-place
68  * --> StorageLimit_ = 1.
69  * - Implicit one-step methods need two solution states --> StorageLimit_ = 2.
70  * - MultiStep methods require k past solution states --> StorageLimit_ = k+1.
71  *
72  * The states contained in the SolutionHistory will often be the last
73  * \f$M\f$ states for forward time integration, e.g.,
74  * \f$[S^0(x^{n-M+1}), ... , S^{M-3}(x^{n-2}), S^{M-2}(x^{n-1}),
75  * S^{M-1}(x^{n})]\f$, but this is not guaranteed.
76  * For example, during transient adjoint sensitivity calculations, not
77  * all the forward states can be stored and therefore will be recalculated.
78  * Thus intermediate states are kept to reduce recomputation costs, i.e.,
79  * so one does not recalculate from the initial conditions to the desired time.
80  * This means the states in the SolutionHistory may not have consecutive
81  * timestep indices, e.g., \f$[..., S^{M-4}(x^{n-200}), S^{M-3}(x^{n-100}),
82  * S^{M-2}(x^{n-1}), S^{M-1}(x^{n})]\f$.
83  *
84  * The SolutionHistory is kept in chronological order, starting with the
85  * oldest state and ending with the latest.
86  *
87  * The "current" state is the latest solution that has been successfully
88  * solved. The solution from the last successful time step or the initial
89  * conditions.
90  *
91  * The "working" state is the state which is being worked on. It is valid
92  * from initialization (e.g., copied from the "current" state), during the
93  * the timestep, and until it is accepted and promoted to the new "current"
94  * state. Between the promotion and the initialization, the workingState_
95  * is invalid (i.e., Teuchos::null). If the timestep fails, the
96  * workingState_ is maintained, and the timestep is retried until the
97  * Integrator declares a successful timestep or the time integration is a
98  * failure. The SolutionHistory indices associated with the currentState_
99  * and the workingState_ vary during the time step loop, and are
100  * <table>
101  * <tr><th> Loop Portion <th> currentState_ <th> workingState_
102  * <tr><td> Before initializing working state
103  * <td> PASS -> \f$M\f$-1 <br>
104  * FAIL -> \f$M\f$-2
105  * <td> PASS -> Invalid <br>
106  * FAIL -> \f$M\f$-1
107  * <tr><td> After initializing working state <br>
108  * During timestep <br>
109  * Before accepting timestep
110  * <td> \f$M\f$-2
111  * <td> \f$M\f$-1
112  * <tr><td> After accepting timestep
113  * <td> PASS -> \f$M\f$-1 <br>
114  * FAIL -> \f$M\f$-2
115  * <td> PASS -> Invalid <br>
116  * FAIL -> \f$M\f$-1
117  * </table>
118  * Initial conditions are considered PASSing, which sets up the loop.
119  * The difference between the currentState_ timestep index and the
120  * workingState_ timestep index is guaranteed to be one (except for when
121  * StorageLimit_ = 1, i.e., explicit one-step methods that can update
122  * the solution in-place).
123  */
124 template<class Scalar>
125 class SolutionHistory
126  : virtual public Teuchos::Describable,
127  virtual public Teuchos::VerboseObject<SolutionHistory<Scalar> >,
128  virtual public Teuchos::ParameterListAcceptor
129 //, virtual public InterpolatorAcceptingObjectBase<Scalar>
130 {
131 public:
132 
133  /// Default Contructor
134  SolutionHistory();
135 
136  /// Contructor
138  std::string name,
139  Teuchos::RCP<std::vector<Teuchos::RCP<SolutionState<Scalar> > > > history,
140  Teuchos::RCP<Interpolator<Scalar> > interpolator,
141  StorageType storageType,
142  int storageLimit);
143 
144  /// Contructor (Slated for deprecation. Use createSolutionHistory(shPL);)
145  SolutionHistory(Teuchos::RCP<Teuchos::ParameterList> shPL);
146 
147  /// Destructor
149 
150  /// \name Basic SolutionHistory Methods
151  //@{
152  /// Add solution state to history
153  void addState(const Teuchos::RCP<SolutionState<Scalar> >& state);
154 
155  /// Add a working solution state to history
156  void addWorkingState(const Teuchos::RCP<SolutionState<Scalar> >& state,
157  const bool updateTime = true);
158 
159  /// Remove solution state
160  void removeState(const Teuchos::RCP<SolutionState<Scalar> >& state);
161 
162  /// Remove solution state based on time
163  void removeState(const Scalar time);
164 
165  /// Find solution state at requested time (no interpolation)
166  Teuchos::RCP<SolutionState<Scalar> > findState(const Scalar time) const;
167 
168  /// Generate and interpolate a new solution state at requested time
169  Teuchos::RCP<SolutionState<Scalar> > interpolateState(const Scalar time) const;
170 
171  /// Interpolate solution state at requested time and store in supplied state
172  void interpolateState(const Scalar time,
173  SolutionState<Scalar>* state_out) const;
174 
175  /// Initialize the working state
176  void initWorkingState();
177 
178  /// Promote the working state to current state
179  void promoteWorkingState();
180 
181  void clear() {history_->clear();}
182  //@}
183 
184  /// \name Accessor methods
185  //@{
186  /// Get this SolutionHistory's name
187  std::string getName() const {return name_;}
188 
189  /// Set this SolutionHistory's name
190  void setName(std::string name) {name_ = name;}
191 
192  /// Get underlining history
193  Teuchos::RCP<std::vector<Teuchos::RCP<SolutionState<Scalar> > > >
194  getHistory() const {return history_;}
195 
196  /// Subscript operator
197  Teuchos::RCP<SolutionState<Scalar> > operator[](const int i) {
198  TEUCHOS_TEST_FOR_EXCEPTION(
199  !((0 <= i) and (i < (int)history_->size())), std::out_of_range,
200  "Error - SolutionHistory index is out of range.\n"
201  << " [Min, Max] = [ 0, " << history_->size()<< "]\n"
202  << " index = " << i << "\n");
203  return (*history_)[i];
204  }
205 
206  /// Subscript operator (const version)
207  Teuchos::RCP<const SolutionState<Scalar> > operator[](const int i) const {
208  TEUCHOS_TEST_FOR_EXCEPTION(
209  !((0 <= i) and (i < (int)history_->size())), std::out_of_range,
210  "Error - SolutionHistory index is out of range.\n"
211  << " [Min, Max] = [ 0, " << history_->size()<< "]\n"
212  << " index = " << i << "\n");
213  return (*history_)[i];
214  }
215 
216  /// Return the current state, i.e., the last accepted state
217  Teuchos::RCP<SolutionState<Scalar> > getCurrentState() const
218  {
219  const int m = history_->size();
220  Teuchos::RCP<SolutionState<Scalar> > state;
221  if (m == 0) state=Teuchos::null;
222  else if (m == 1 or workingState_ == Teuchos::null) state=(*history_)[m-1];
223  else if (m > 1) state=(*history_)[m-2];
224  return state;
225  }
226 
227  /// Return the working state
228  Teuchos::RCP<SolutionState<Scalar> > getWorkingState(bool warn=true) const
229  {
230  if (workingState_ == Teuchos::null && warn) {
231  Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
232  Teuchos::OSTab ostab(out,1,"SolutionHistory::getWorkingState()");
233  *out << "Warning - WorkingState is null and likely has been promoted. "
234  << "You might want to call getCurrentState() instead.\n"
235  << std::endl;
236  }
237  return workingState_;
238  }
239 
240  /// Get the number of states
241  int getNumStates() const {return history_->size();}
242 
243  /// Get the current time
244  Scalar getCurrentTime() const {return getCurrentState()->getTime();}
245 
246  /// Get the current timestep index
247  int getCurrentIndex() const {return getCurrentState()->getIndex();}
248 
249  /// Set the maximum storage of this history
250  void setStorageLimit(int storage_limit);
251 
252  /// Get the maximum storage of this history
253  int getStorageLimit() const {return storageLimit_;}
254 
257 
258  /// Return the current minimum time of the SolutionStates
259  Scalar minTime() const {return (history_->front())->getTime();}
260 
261  /// Return the current maximum time of the SolutionStates
262  Scalar maxTime() const {return (history_->back())->getTime();}
263 
264  /// Get the state with timestep index equal to n
265  Teuchos::RCP<SolutionState<Scalar> > getStateTimeIndexN() const;
266 
267  /// Get the state with timestep index equal to n-1
268  Teuchos::RCP<SolutionState<Scalar> > getStateTimeIndexNM1() const;
269 
270  /// Get the state with timestep index equal to n-2
271  Teuchos::RCP<SolutionState<Scalar> > getStateTimeIndexNM2() const;
272 
273  /// Get the state with timestep index equal to "index"
274  Teuchos::RCP<SolutionState<Scalar> > getStateTimeIndex(int index) const;
275  //@}
276 
277  /// \name Overridden from Teuchos::ParameterListAcceptor
278  //@{
279  void setParameterList(const Teuchos::RCP<Teuchos::ParameterList> & pl);
280  Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
281  Teuchos::RCP<Teuchos::ParameterList> getNonconstParameterList();
282  Teuchos::RCP<Teuchos::ParameterList> unsetParameterList();
283  //@}
284 
285  /// \name Overridden from Teuchos::Describable
286  //@{
287  virtual std::string description() const;
288  virtual void describe(Teuchos::FancyOStream & out,
289  const Teuchos::EVerbosityLevel verbLevel) const;
290  //@}
291 
292  /// \name Interpolation Methods
293  //@{
294  /// Set the interpolator for this history
295  void setInterpolator(const Teuchos::RCP<Interpolator<Scalar> >& interpolator);
296  Teuchos::RCP<Interpolator<Scalar> > getNonconstInterpolator();
297  Teuchos::RCP<const Interpolator<Scalar> > getInterpolator() const;
298  /// Unset the interpolator for this history
299  Teuchos::RCP<Interpolator<Scalar> > unSetInterpolator();
300  //@}
301 
302  void printHistory(std::string verb="low") const
303  {
304  Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
305  Teuchos::OSTab ostab(out,1,"SolutionHistory::printHistory");
306  *out << name_ << " (size=" << history_->size() << ")"
307  << " (w - working; c - current; i - interpolated)" << std::endl;
308  for (int i=0; i<(int)history_->size() ; ++i) {
309  auto state = (*history_)[i];
310  *out << " ";
311  if (state == getWorkingState()) *out << "w - ";
312  else if (state == getCurrentState()) *out << "c - ";
313  else if (state->getIsInterpolated() == true) *out<<"i - ";
314  else *out << " ";
315  *out << "[" << i << "] = " << state << std::endl;
316  if (verb == "medium" or verb == "high") {
317  if (state != Teuchos::null) {
318  auto x = state->getX();
319  *out << " x = " << x << std::endl
320  << " norm(x) = " << Thyra::norm(*x) << std::endl;
321  }
322  }
323  if (verb == "high") {
324  (*history_)[i]->describe(*out,this->getVerbLevel());
325  }
326  }
327  }
328 
329 protected:
330 
331  std::string name_;
332  Teuchos::RCP<Teuchos::ParameterList> shPL_;
333  Teuchos::RCP<std::vector<Teuchos::RCP<SolutionState<Scalar> > > > history_;
334  Teuchos::RCP<Interpolator<Scalar> > interpolator_;
337 
338  Teuchos::RCP<SolutionState<Scalar> > workingState_; ///< The state being worked on
339 };
340 
341 
342 /// Nonmember constructor from a ParameterList
343 template<class Scalar>
344 Teuchos::RCP<SolutionHistory<Scalar> >
345 createSolutionHistoryPL(Teuchos::RCP<Teuchos::ParameterList> pList);
346 
347 /// Nonmember contructor from a SolutionState.
348 template<class Scalar>
349 Teuchos::RCP<SolutionHistory<Scalar> >
350 createSolutionHistoryState(const Teuchos::RCP<SolutionState<Scalar> >& state);
351 
352 /// Nonmember contructor from a Thyra ModelEvaluator.
353 template<class Scalar>
354 Teuchos::RCP<SolutionHistory<Scalar> >
356  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model);
357 
358 
359 } // namespace Tempus
360 
361 #endif // Tempus_SolutionHistory_decl_hpp
Teuchos::RCP< SolutionState< Scalar > > getStateTimeIndexNM2() const
Get the state with timestep index equal to n-2.
Teuchos::RCP< Interpolator< Scalar > > unSetInterpolator()
Unset the interpolator for this history.
Scalar minTime() const
Return the current minimum time of the SolutionStates.
std::string getName() const
Get this SolutionHistory&#39;s name.
Keep the 2 newest states for undo.
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList()
int getCurrentIndex() const
Get the current timestep index.
Teuchos::RCP< Interpolator< Scalar > > interpolator_
void initWorkingState()
Initialize the working state.
void printHistory(std::string verb="low") const
Teuchos::RCP< SolutionState< Scalar > > interpolateState(const Scalar time) const
Generate and interpolate a new solution state at requested time.
int getNumStates() const
Get the number of states.
Teuchos::RCP< const SolutionState< Scalar > > operator[](const int i) const
Subscript operator (const version)
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList()
Teuchos::RCP< Interpolator< Scalar > > getNonconstInterpolator()
Scalar getCurrentTime() const
Get the current time.
void setStorageLimit(int storage_limit)
Set the maximum storage of this history.
Teuchos::RCP< SolutionHistory< Scalar > > createSolutionHistoryPL(Teuchos::RCP< Teuchos::ParameterList > pList)
Nonmember constructor from a ParameterList.
Teuchos::RCP< SolutionHistory< Scalar > > createSolutionHistoryME(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model)
Nonmember contructor from a Thyra ModelEvaluator.
Teuchos::RCP< SolutionState< Scalar > > getWorkingState(bool warn=true) const
Return the working state.
Teuchos::RCP< SolutionState< Scalar > > workingState_
The state being worked on.
Teuchos::RCP< SolutionState< Scalar > > getStateTimeIndexNM1() const
Get the state with timestep index equal to n-1.
void promoteWorkingState()
Promote the working state to current state.
void removeState(const Teuchos::RCP< SolutionState< Scalar > > &state)
Remove solution state.
Teuchos::RCP< SolutionState< Scalar > > getCurrentState() const
Return the current state, i.e., the last accepted state.
Teuchos::RCP< std::vector< Teuchos::RCP< SolutionState< Scalar > > > > getHistory() const
Get underlining history.
Base strategy class for interpolation functionality.
virtual std::string description() const
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &pl)
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
void addState(const Teuchos::RCP< SolutionState< Scalar > > &state)
Add solution state to history.
Teuchos::RCP< Teuchos::ParameterList > shPL_
void setName(std::string name)
Set this SolutionHistory&#39;s name.
void addWorkingState(const Teuchos::RCP< SolutionState< Scalar > > &state, const bool updateTime=true)
Add a working solution state to history.
Keep a fix number of states.
Teuchos::RCP< SolutionHistory< Scalar > > createSolutionHistoryState(const Teuchos::RCP< SolutionState< Scalar > > &state)
Nonmember contructor from a SolutionState.
Teuchos::RCP< SolutionState< Scalar > > getStateTimeIndex(int index) const
Get the state with timestep index equal to &quot;index&quot;.
Teuchos::RCP< SolutionState< Scalar > > getStateTimeIndexN() const
Get the state with timestep index equal to n.
Teuchos::RCP< const Interpolator< Scalar > > getInterpolator() const
Scalar maxTime() const
Return the current maximum time of the SolutionStates.
void setInterpolator(const Teuchos::RCP< Interpolator< Scalar > > &interpolator)
Set the interpolator for this history.
int getStorageLimit() const
Get the maximum storage of this history.
Teuchos::RCP< SolutionState< Scalar > > findState(const Scalar time) const
Find solution state at requested time (no interpolation)
Teuchos::RCP< SolutionState< Scalar > > operator[](const int i)
Subscript operator.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Solution state for integrators and steppers. SolutionState contains the metadata for solutions and th...
Teuchos::RCP< std::vector< Teuchos::RCP< SolutionState< Scalar > > > > history_