OGS
EvolutionaryPIDcontroller.cpp
Go to the documentation of this file.
1 
13 
14 #include <functional>
15 #include <limits>
16 #include <vector>
17 
18 #include "BaseLib/Algorithm.h"
19 #include "BaseLib/Logging.h"
20 
21 namespace NumLib
22 {
23 std::tuple<bool, double> EvolutionaryPIDcontroller::next(
24  double const solution_error, int const /*number_iterations*/)
25 {
26  const bool is_previous_step_accepted = _is_accepted;
27 
28  const double e_n = solution_error;
29  const double zero_threshlod = std::numeric_limits<double>::epsilon();
30  // step rejected.
31  if (e_n > _tol)
32  {
33  _is_accepted = false;
34 
35  double h_new = (e_n > zero_threshlod) ? _ts_current.dt() * _tol / e_n
36  : 0.5 * _ts_current.dt();
37 
38  h_new = limitStepSize(h_new, is_previous_step_accepted);
39 
40  WARN(
41  "This step is rejected due to the relative change from the"
42  " solution of the previous\n"
43  "\t time step to the current solution exceeds the given tolerance"
44  " of {:g}.\n"
45  "\t This time step will be repeated with a new time step size of"
46  " {:g}\n"
47  "\t or the simulation will be halted.",
48  _tol, h_new);
49 
50  return std::make_tuple(false, h_new);
51  }
52 
53  // step accepted.
54  _is_accepted = true;
55 
56  if (_ts_current.steps() == 0)
57  {
58  _e_n_minus1 = e_n;
59 
60  return std::make_tuple(true, _h0);
61  }
62  else
63  {
64  const double h_n = _ts_current.dt();
65  double h_new = h_n;
66 
67  if (e_n > zero_threshlod)
68  {
69  if (_e_n_minus1 > zero_threshlod)
70  {
71  if (_e_n_minus2 > zero_threshlod)
72  {
73  h_new = std::pow(_e_n_minus1 / e_n, _kP) *
74  std::pow(_tol / e_n, _kI) *
75  std::pow(
77  _kD) *
78  h_n;
79  }
80  else
81  {
82  h_new = std::pow(_e_n_minus1 / e_n, _kP) *
83  std::pow(_tol / e_n, _kI) * h_n;
84  }
85  }
86  else
87  {
88  h_new = std::pow(_tol / e_n, _kI) * h_n;
89  }
90  }
91 
92  h_new = limitStepSize(h_new, is_previous_step_accepted);
93 
95  _e_n_minus1 = e_n;
96 
97  return std::make_tuple(true, h_new);
98  }
99 
100  return {};
101 }
102 
104  const double h_new, const bool previous_step_accepted) const
105 {
106  const double h_n = _ts_current.dt();
107  // Forced the computed time step size in the given range
108  // (see the formulas in the documentation of the class)
109  const double h_in_range = std::max(_h_min, std::min(h_new, _h_max));
110  // Limit the step size change ratio.
111  double limited_h =
112  std::max(_rel_h_min * h_n, std::min(h_in_range, _rel_h_max * h_n));
113 
114  if (!previous_step_accepted)
115  {
116  // If the last time step was rejected and the new predicted time step
117  // size is identical to that of the previous rejected step, the new
118  // step size is then reduced by half.
119  if (std::abs(limited_h - _ts_current.dt()) <
120  std::numeric_limits<double>::min())
121  {
122  limited_h = std::max(_h_min, 0.5 * limited_h);
123  }
124 
125  // If the last time step was rejected and the new predicted time step
126  // size is larger than the step size of the rejected step, the new step
127  // size takes the half of the size of the rejected step. This could
128  // happen when a time step is rejected due to a diverged non-linear
129  // solver. In such case, this algorithm may give a large time step size
130  // by using the diverged solution.
131  if (limited_h > _ts_current.dt())
132  {
133  limited_h = std::max(_h_min, 0.5 * _ts_current.dt());
134  }
135  }
136  return limited_h;
137 }
138 
140 {
141  // If current and previous dt are both at minimum dt, then cannot reduce
142  // further.
143  return !(_ts_current.dt() == _h_min && _ts_prev.dt() == _h_min);
144 }
145 } // namespace NumLib
void WARN(char const *fmt, Args const &... args)
Definition: Logging.h:37
const double _h_max
maximum step size.
virtual bool canReduceTimestepSize() const override
Query the timestepper if further time step size reduction is possible.
const double _h_min
minimum step size.
std::tuple< bool, double > next(double solution_error, int number_iterations) override
double limitStepSize(const double h_new, const bool previous_step_accepted) const
const double _h0
initial time step size.
TimeStep _ts_current
current time step information
TimeStep _ts_prev
previous time step information
std::size_t steps() const
the number of time _steps
Definition: TimeStep.h:113
double dt() const
time step size from _previous
Definition: TimeStep.h:111