Flow123d  release_2.2.0-23-g01927c6
time_governor.hh
Go to the documentation of this file.
1 /*!
2  *
3  * Copyright (C) 2015 Technical University of Liberec. All rights reserved.
4  *
5  * This program is free software; you can redistribute it and/or modify it under
6  * the terms of the GNU General Public License version 3 as published by the
7  * Free Software Foundation. (http://www.gnu.org/licenses/gpl-3.0.en.html)
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
12  *
13  *
14  * @file time_governor.hh
15  * @brief Basic time management class.
16  * @author Jan Brezina
17  */
18 
19 #ifndef TIME_HH_
20 #define TIME_HH_
21 
22 #include <limits>
23 #include <cmath>
24 #include <algorithm>
25 #include <boost/circular_buffer.hpp>
26 
27 
28 #include "system/global_defs.h"
29 #include "system/system.hh"
31 #include "input/input_exception.hh"
32 #include "tools/time_marks.hh"
33 
34 namespace Input {
35  class Record;
36  namespace Type {
37  class Record;
38  }
39 }
40 
41 
42 
43 
44 /**
45  * @brief Representation of one time step.\
46  *
47  * Time step consists of the time step @p length() and from time step @end() time.
48  * More over we store the index of the time step within it time governor.
49  *
50  * The reason to store both the end time and the length of the time step is to allow
51  * safe comparisons of the time with safety margin small relative to the
52  * time step length.
53  */
54 class TimeStep {
55 public:
56  /**
57  * Constructor of the zero time step.
58  */
59  TimeStep(double init_time);
60 
61  /**
62  * Default constructor.
63  * Creates undefined time step.
64  */
65  TimeStep();
66 
67  /**
68  * Copy constructor.
69  */
70  TimeStep(const TimeStep &other);
71 
72  /**
73  * Create subsequent time step.
74  */
75  TimeStep make_next(double new_length) const;
76 
77  /**
78  * Create subsequent time step, with the @end_time
79  * explicitly specified. This allow slight discrepancy to
80  * overcome rounding errors in the case of fixed time step.
81  * Otherwise using small fixed time step, we may miss long term fixed
82  * goal time.
83  *
84  */
85  TimeStep make_next(double new_lenght, double end_time) const;
86 
87  /**
88  * Getters.
89  */
90  unsigned int index() const {return index_;}
91  double length() const { return length_;}
92  double end() const { return end_;}
93  /**
94  * Performs rounding safe comparison time > other_time, i.e. time is strictly greater than given parameter
95  * other_time with precision relative to the magnitude of the numbers time step.
96  * TODO: introduce type TimeDouble with overloaded comparison operators, use it consistently in TimeMarks.
97  */
98  inline bool gt(double other_time) const
99  { return ! safe_compare(other_time, end());}
100 
101  /**
102  * Performs rounding safe comparison time >= other_time See @fn gt
103  */
104  inline bool ge(double other_time) const
105  { return safe_compare(end(), other_time); }
106 
107  /**
108  * Performs rounding safe comparison time < other_time. See @fn gt
109  */
110  inline bool lt(double other_time) const
111  { return ! safe_compare(end(), other_time); }
112 
113  /**
114  * Performs rounding safe comparison time <= other_time. See @fn gt
115  */
116  inline bool le(double other_time) const
117  { return safe_compare(other_time, end()); }
118 
119  /**
120  * Performs rounding safe comparison time (step) == other_time. See @fn gt
121  */
122  inline bool eq(double other_time) const
123  { return this->le(other_time) && this->ge(other_time); }
124 
125  inline bool contains(double other_time) const
126  { return this->ge(other_time) && this->lt(other_time + length_); }
127 
128  /**
129  * Returns true if two time steps are exactly the same.
130  */
131  bool operator==(const TimeStep & other)
132  { return (index_ == other.index_)
133  && (length_ == other.length_)
134  && (end_ == other.end_);
135  }
136 private:
137 
138  /* Returns true if t1-t0 > delta. Where delta is choosen
139  * related to the current time step and magnitude of t1, t0.
140  */
141  bool safe_compare(double t1, double t0) const;
142 
143  /// Index of the step is index if the end time. Zero time step is artificial.
144  unsigned int index_;
145  /// Length of the time step. Theoretically @p end minus end of the previous time step.
146  /// However may be slightly different due to rounding errors.
147  double length_;
148  /// End time point of the time step.
149  double end_;
150 };
151 
152 std::ostream& operator<<(std::ostream& out, const TimeStep& t_step);
153 
154 
155 
156 /**
157  * @brief
158  * Basic time management functionality for unsteady (and steady) solvers (class Equation).
159  *
160  * <h2> Common features and unsteady time governor (TG) </h2>
161  *
162  * This class provides algorithm for selecting next time step, and information about current time step frame.
163  * Step estimating is constrained by several bounds (permanent maximal and minimal time step, upper
164  * and lower constraint of time step). The permanent constraints are set in the constructor from the input
165  * record so that user can set the time step constraints for the whole simulation.
166  * Function set_permanent_constraint() should be used only in very specific cases and possibly right after
167  * the constructor before using other functions of TG.
168  *
169  * Choice of the very next time step can be constrained using functions set_upper_constraint()
170  * and set_lower_constraint().
171  * Lower and upper constraints are set equal to permanent ones in the constructor and can only
172  * become stricter. If one tries to set these constraints outside the interval of the previous constraints,
173  * nothing is changed and a specified value is returned. Upper and lower constraints are reset in function
174  * next_time() to the permanent constraints.
175  *
176  * The later one can be called multiple times with various constraint values and we use the minimum of them.
177  * Function next_time() choose the next time step in such a way that it meets actual constraints and
178  * a uniform discrete time grid with this step hits the nearest fixed time in lowest possible number of steps.
179  *
180  * The fixed times are time marks of TimeMarks object passed at construction time with particular mask.
181  *
182  * There is just one set of time marks for the whole problem. Therefore TimeMarks object is static and is shared umong
183  * all the equations and time governors. Each equation creates its own specific time mark type.
184  *
185  * Information provided by TG includes:
186  * - actual time, last time, end time
187  * - actual time step
188  * - number of the time level
189  * - end of interval with fixed time step
190  * - time comparison
191  * - static pointer to time marks
192  *
193  * <h2> Steady time governor</h2>
194  *
195  * Steady TG can be constructed by default constructor (initial time is zero) or by
196  * constructor with initial time as parameter. End time and time step are set to infinity.
197  * One can check if the time governor is steady by calling is_steady().
198  * Calling estimate_dt() will return infinity.
199  *
200  * Setting constraints have no consequences. Calling fix_dt_until_mark() will only return zero
201  * and will not do anything.
202  *
203  * The steady TG works in two states. At first the time is set to initial and time level
204  * is equal zero. To use steady TG properly one should call next_time() after the computation
205  * of steady problem is done. Current time is then set to infinity, time level is set to 1 and
206  * calling estimate_dt() will return zero.
207  *
208  * Note: For example class TransportNothing (which computes really nothing) uses also steady TG but
209  * it calls next_time() immediately after TG's construction. This means that the 'computation'of transport
210  * is done.
211  *
212  *
213  */
214 
216 {
217 public:
218 
219  DECLARE_INPUT_EXCEPTION(ExcTimeGovernorMessage, << EI_Message::val);
220  TYPEDEF_ERR_INFO( EI_Index, int);
221  TYPEDEF_ERR_INFO( EI_BackIndex, unsigned int);
222  TYPEDEF_ERR_INFO( EI_HistorySize, unsigned int);
223  DECLARE_EXCEPTION(ExcMissingTimeStep,
224  << "Time step index: " << EI_Index::val
225  << ", history index: " << EI_BackIndex::val
226  << " out of history of size: " << EI_HistorySize::val);
227 
228  static const Input::Type::Record & get_input_type();
229 
230  /**
231  * Getter for time marks.
232  */
233  static inline TimeMarks &marks()
234  {return time_marks_;}
235 
236  /**
237  * @brief Constructor for unsteady solvers.
238  *
239  * @param input accessor to input data
240  * @param fixed_time_mask TimeMark mask used to select fixed time marks from all the time marks.
241  * This value is bitwise added to the default one defined in TimeMarks::type_fixed_time().
242  *
243  */
244  TimeGovernor(const Input::Record &input,
245  TimeMark::Type fixed_time_mask = TimeMark::none_type);
246 
247  /**
248  * @brief Default constructor - steady time governor.
249  *
250  * OBSOLETE.
251  *
252  * We can have "zero step" steady problem (no computation, e.g. EquationNothing) and one step steady problem
253  * (e.g. steady water flow).
254  *
255  * Time is set to zero, time step and end time to infinity.
256  *
257  * First call of next_time() pushes the actual time to infinity.
258  *
259  * However, you have to use full constructor for the "steady problem" that has time-variable input data.
260  *
261  * Has a private pointer to static TimeMarks and can access them by marks().
262  */
263  explicit TimeGovernor(double init_time=0.0,
264  TimeMark::Type fixed_time_mask = TimeMark::none_type);
265 
266  /**
267  * The aim of this constuctor is simple way to make a time governor without Input interface.
268  *
269  * TODO: Partially tested as part of field test. Needs its own unit test.
270  */
271  TimeGovernor(double init_time, double dt);
272 
273  /**
274  * Returns true if the time governor was set from default values
275  */
276  bool is_default() {
277  return (end_time_ == max_end_time)
278  && (max_time_step_ == end_time_ - init_time_);
279  }
280 
281  /**
282  * @brief Sets permanent constraints for time step.
283  *
284  * This function should not be normally used. These values are to be set in constructor
285  * from the input record or by default.
286  * @param min_dt is the minimal value allowed for time step
287  * @param max_dt is the maximal value allowed for time step
288  */
289  void set_permanent_constraint( double min_dt, double max_dt);
290 
291  /**
292  * @brief Sets upper constraint for the next time step estimating.
293  *
294  * This function can only make the constraint stricter. Upper constraint is reset to @p max_dt after the next_time() call.
295  * The return value mimics result of the comparison: current constraint compared to @param upper.
296  * @param message describes the origin of the constraint
297  * In particular the return values is:
298  * - -1: Current upper constraint is less then the @param upper. No change happen.
299  * - 0: Current constraint interval contains the @param upper. The upper constraint is set.
300  * - +1: Current lower constraint is greater then the @param upper. No change happen.
301  */
302  int set_upper_constraint(double upper, std::string message);
303 
304  /**
305  * @brief Sets lower constraint for the next time step estimating.
306  *
307  * This function can only make the constraint stricter. Lower constraint is reset to @p min_dt after the next_time() call.
308  * The return value mimics result of the comparison: current constraint compared to @param upper.
309  * In particular the return values is:
310  * - -1: Current upper constraint is less then the @param lower. No change happen.
311  * - 0: Current constraint interval contains the @param lower. The lower constraint is set.
312  * - +1: Current upper constraint is greater then the @param lower. No change happen.
313  */
314  /**
315  * @brief Sets lower constraint for the next time step estimating.
316  * @param lower is the lower constraint for time step
317  * @param message describes the origin of the constraint
318  * @return -1, 0 or 1 according to the success.
319  * @see set_upper_constrain().
320  */
321  int set_lower_constraint(double lower, std::string message);
322 
323  /**
324  * @brief Fixing time step until fixed time mark.
325  *
326  * Fix time step until first fixed time mark. When called inside an already fixed interval,
327  * it overwrites previous setting.
328  * @return actual end of fixed time step.
329  */
330  double fix_dt_until_mark();
331 
332  /**
333  * @brief Proceed to the next time according to current estimated time step.
334  *
335  * The timestep constraints are relaxed to the permanent constraints.
336  */
337  void next_time();
338 
339  /**
340  * @brief Force timestep reduction in particular in the case of failure of the non-linear solver.
341  *
342  * Calling this method also force immediate end of the fixed timestep interval.
343  * Returns true reduce factor used. It is larger then given factor if we hit the lower timestep constraint.
344  *
345  * TODO: How to keep constraints active for the last next_time call.
346  */
347  double reduce_timestep(double factor);
348 
349  /**
350  * Returns reference to required time step in the recent history.
351  * Without parameter the actual time step is returned.
352  * Use negative indices to get recent time steps: step(-1) the actual step, step(-2) the last one.
353  * Use positive index to get time step by its index: step(0) the first time step.
354  * However only limited number of last time steps is stored.
355  * If the time step is not accessible any more, we throw an exception ExcMissingTimeStep.
356  */
357  const TimeStep &step(int index=-1) const;
358 
359  /**
360  * Specific time mark of the equation owning the time governor.
361  */
363  { return eq_mark_type_;}
364 
365  /**
366  * Specific time mark of the fixed times of the equation owning the time governor.
367  */
369  { return eq_mark_type_ | marks().type_fixed_time(); }
370 
371  /**
372  * Add sequence of time marks starting from the initial time up to the end time with given @p step.
373  * Time marks type combines given mark_type (none by default) and native mark type of the time governor.
374  */
375  void add_time_marks_grid(double step, TimeMark::Type mark_type= TimeMark::none_type) const;
376 
377  /**
378  * Simpler interface to TimeMarks::is_current().
379  */
380  bool is_current(const TimeMark::Type &mask) const;
381 
382 
383  /**
384  * Simpler interface to TimeMarks::next().
385  */
386  inline TimeMarks::iterator next(const TimeMark::Type &mask) const
387  {return time_marks_.next(*this, mask);}
388 
389  /**
390  * Simpler interface to TimeMarks::last().
391  */
392  inline TimeMarks::iterator last(const TimeMark::Type &mask) const
393  {return time_marks_.last(*this, mask);}
394 
395  /**
396  * Getter for upper constrain.
397  */
398  inline double upper_constraint() const
399  {return upper_constraint_;}
400 
401  /**
402  * Returns lower constraint.
403  */
404  inline double lower_constraint() const
405  {return lower_constraint_;}
406 
407  /**
408  * End of interval with currently fixed time step. Can be changed by next call of method fix_dt_until_mark.
409  */
410  inline double end_of_fixed_dt() const
411  {return end_of_fixed_dt_interval_;}
412 
413  /**
414  * Getter for dt_changed. Returns whether the time step has been changed.
415  */
416  inline bool is_changed_dt() const
417  {return time_step_changed_;}
418 
419 
420  /**
421  * Initial time getter.
422  */
423  inline double init_time() const
424  {return this->init_time_;}
425 
426  /**
427  * End of actual time interval; i.e. where the solution is computed.
428  */
429  inline double t() const
430  {return step().end();}
431 
432  /**
433  * Previous time step.
434  */
435  inline double last_dt() const
436  {if (step().index() >0) return step(-2).length();
437  else return inf_time;
438  }
439 
440  /**
441  * Previous time.
442  */
443  inline double last_t() const
444  {if (step().index() >0) return step(-2).end();
445  else return step().end() - step().length();
446  }
447 
448 
449  /**
450  * Length of actual time interval; i.e. the actual time step.
451  */
452  inline double dt() const
453  {return step().length();}
454 
455  /**
456  * @brief Estimate choice of next time step according to actual setting of constraints.
457  *
458  * Precedence of constraints:
459  *
460  * -# meet next fixed time (always satisfied)
461  * -# permanent upper constraint (always satisfied)
462  * -# upper constraint (always satisfied)
463  * -# lower constraint (satisfied if not in conflict with 1.)
464  * -# permanent lower constraint (satisfied if 4.)
465  * -# else writes the difference between lower constraint and estimated time step
466  * -# If there are more then one step to the next fixed time, try to
467  * use the last time step if it is nearly the same.
468  */
469  double estimate_dt() const;
470 
471  /**
472  * Estimate next time.
473  */
474  inline double estimate_time() const
475  {return t()+estimate_dt();}
476 
477  /// End time.
478  inline double end_time() const
479  { return end_time_; }
480 
481  /// Returns true if the actual time is greater than or equal to the end time.
482  inline bool is_end() const
483  { return (this->step().ge(end_time_) || t() == inf_time); }
484 
485  /// Returns true if the time governor is used for steady problem.
486  inline bool is_steady() const
487  { return steady_; }
488 
489 
490 
491  /**
492  * Returns the time level.
493  */
494  inline int tlevel() const
495  {return step().index();}
496 
497  /**
498  * Prints output of TimeGovernor.
499  * @param name is the name of time governor that you want to show up in output (just for your convenience)
500  *
501  */
502  void view(const char *name="") const;
503 
504  // Maximal tiem of simulation. More then age of the universe in seconds.
505  static const double max_end_time;
506 
507  /// Infinity time used for steady case.
508  static const double inf_time;
509 
510  /**
511  * Rounding precision for computing time_step.
512  * Used as technical lower bound for the time step.
513  */
514  static const double time_step_precision;
515 
516 private:
517 
518 
519  /**
520  * \brief Common part of the constructors. Set most important parameters, check they are valid and set default values to other.
521  *
522  * Set main parameters to given values.
523  * Check they are correct.
524  * Set soft and permanent constrains to the same, the least restricting values.
525  * Set time marks for the start time and end time (if finite).
526  */
527  void init_common(double init_time, double end_time, TimeMark::Type type);
528 
529  /**
530  * Size of the time step buffer, i.e. recent_time_steps_.
531  */
532  static const unsigned int size_of_recent_steps_ = 3;
533 
534  /// Circular buffer of recent time steps. Implicit size is 3.
535  boost::circular_buffer<TimeStep> recent_steps_;
536  /// Initial time.
537  double init_time_;
538  /// End of interval if fixed time step.
540  /// End time of the simulation.
541  double end_time_;
542 
543  /// Next fixed time step.
545  /// Flag that is set when the fixed step is set (lasts only one time step).
547  /// Flag is set if the time step has been changed (lasts only one time step).
549 
550  /// Description of the upper constraint.
552  /// Description of the upper constraint.
554  /// Upper constraint for the choice of the next time step.
556  /// Lower constraint for the choice of the next time step.
558  /// Permanent upper limit for the time step.
560  /// Permanent lower limit for the time step.
562 
563  /// Upper constraint used for choice of current time.
565  /// Lower constraint used for choice of current time.
567 
568  /**
569  * When the next time is chosen we need only the lowest fix time. Therefore we use
570  * minimum priority queue of doubles based on the vector container.
571  * This is one global set of time marks for the whole problem and is shared among all equations.
572  * Therefore this object is static constant pointer.
573  */
575 
576  /// TimeMark type of the equation.
578 
579  /// True if the time governor is used for steady problem.
580  bool steady_;
581 
582  friend TimeMarks;
583 };
584 
585 /**
586  * \brief Stream output operator for TimeGovernor.
587  *
588  * Currently for debugging purposes.
589  * In the future it should be customized for use in combination with
590  * streams for various log targets.
591  *
592  */
593 std::ostream& operator<<(std::ostream& out, const TimeGovernor& tg);
594 
595 
596 
597 #endif /* TIME_HH_ */
bool operator==(const TimeStep &other)
double init_time_
Initial time.
double lower_constraint_
Lower constraint for the choice of the next time step.
Iterator over TimeMark objects in TimeMarks object (database of TimeMark objects).
Definition: time_marks.hh:342
double end_time() const
End time.
double estimate_time() const
bool is_time_step_fixed_
Flag that is set when the fixed step is set (lasts only one time step).
boost::circular_buffer< TimeStep > recent_steps_
Circular buffer of recent time steps. Implicit size is 3.
int tlevel() const
static TimeMarks time_marks_
bool lt(double other_time) const
TYPEDEF_ERR_INFO(EI_KeyName, const string)
double upper_constraint() const
double end_time_
End time of the simulation.
bool contains(double other_time) const
bool is_end() const
Returns true if the actual time is greater than or equal to the end time.
TimeMarks::iterator last(const TimeMark::Type &mask) const
TimeMarks::iterator next(const TimeMark::Type &mask) const
double end_of_fixed_dt() const
static const double max_end_time
double t() const
Basic time management functionality for unsteady (and steady) solvers (class Equation).
static const Type none_type
Mark Type with all bits unset.
Definition: time_marks.hh:86
static TimeMarks & marks()
bool time_step_changed_
Flag is set if the time step has been changed (lasts only one time step).
double last_t() const
bool le(double other_time) const
double upper_constraint_
Upper constraint for the choice of the next time step.
double init_time() const
TimeMark::Type eq_mark_type_
TimeMark type of the equation.
Global macros to enhance readability and debugging, general constants.
bool ge(double other_time) const
static const double time_step_precision
TimeMark::Type equation_fixed_mark_type() const
Accessor to the data with type Type::Record.
Definition: accessors.hh:292
TimeMark::Type equation_mark_type() const
bool is_changed_dt() const
double length() const
double last_lower_constraint_
Lower constraint used for choice of current time.
double lower_constraint() const
unsigned int index_
Index of the step is index if the end time. Zero time step is artificial.
double end() const
This class is a collection of time marks to manage various events occurring during simulation time...
Definition: time_marks.hh:197
double end_of_fixed_dt_interval_
End of interval if fixed time step.
bool steady_
True if the time governor is used for steady problem.
bool is_steady() const
Returns true if the time governor is used for steady problem.
double max_time_step_
Permanent upper limit for the time step.
double dt() const
std::string upper_constraint_message_
Description of the upper constraint.
bool gt(double other_time) const
DECLARE_EXCEPTION(ExcWrongDefaultJSON,<< "Consistency Error: Not valid JSON of Default value "<< EI_DefaultStr::qval<< " of type "<< EI_TypeName::qval<< ";\n"<< "During declaration of the key: "<< EI_KeyName::qval)
std::ostream & operator<<(std::ostream &stream, const TypeBase &type)
For convenience we provide also redirection operator for output documentation of Input:Type classes...
Definition: type_base.cc:240
unsigned int index() const
Record type proxy class.
Definition: type_record.hh:182
double length_
double fixed_time_step_
Next fixed time step.
DECLARE_INPUT_EXCEPTION(ExcInputMessage,<< EI_Message::val)
Simple input exception that accepts just string message.
double last_upper_constraint_
Upper constraint used for choice of current time.
double last_dt() const
static const double inf_time
Infinity time used for steady case.
Representation of one time step..
double end_
End time point of the time step.
std::string lower_constraint_message_
Description of the upper constraint.
bool eq(double other_time) const
double min_time_step_
Permanent lower limit for the time step.