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
21namespace NumLib
22{
23std::tuple<bool, double> EvolutionaryPIDcontroller::next(
24 double const solution_error, int const /*number_iterations*/,
25 NumLib::TimeStep& /*timestep_previous*/, NumLib::TimeStep& timestep_current)
26{
27 const bool is_previous_step_accepted = timestep_current.isAccepted();
28
29 const double e_n = solution_error;
30 const double zero_threshold = std::numeric_limits<double>::epsilon();
31 // step rejected.
32 if (e_n > _tol)
33 {
34 timestep_current.setAccepted(false);
35
36 double h_new = (e_n > zero_threshold)
37 ? timestep_current.dt() * _tol / e_n
38 : 0.5 * timestep_current.dt();
39
40 h_new =
41 limitStepSize(h_new, is_previous_step_accepted, timestep_current);
42
43 WARN(
44 "This step is rejected due to the relative change from the"
45 " solution of the previous\n"
46 "\t time step to the current solution exceeds the given tolerance"
47 " of {:g}.\n"
48 "\t This time step will be repeated with a new time step size of"
49 " {:g}\n"
50 "\t or the simulation will be halted.",
51 _tol, h_new);
52
53 return std::make_tuple(timestep_current.isAccepted(), h_new);
54 }
55
56 // step accepted.
57 timestep_current.setAccepted(true);
58
59 if (timestep_current.timeStepNumber() == 0)
60 {
61 _e_n_minus1 = e_n;
62
63 return std::make_tuple(timestep_current.isAccepted(), _h0);
64 }
65 else
66 {
67 const double h_n = timestep_current.dt();
68 double h_new = h_n;
69
70 if (e_n > zero_threshold)
71 {
72 if (_e_n_minus1 > zero_threshold)
73 {
74 if (_e_n_minus2 > zero_threshold)
75 {
76 h_new = std::pow(_e_n_minus1 / e_n, _kP) *
77 std::pow(_tol / e_n, _kI) *
78 std::pow(
80 _kD) *
81 h_n;
82 }
83 else
84 {
85 h_new = std::pow(_e_n_minus1 / e_n, _kP) *
86 std::pow(_tol / e_n, _kI) * h_n;
87 }
88 }
89 else
90 {
91 h_new = std::pow(_tol / e_n, _kI) * h_n;
92 }
93 }
94
95 h_new =
96 limitStepSize(h_new, is_previous_step_accepted, timestep_current);
97
99 _e_n_minus1 = e_n;
100
101 return std::make_tuple(timestep_current.isAccepted(), h_new);
102 }
103}
104
106 const double h_new, const bool previous_step_accepted,
107 NumLib::TimeStep const& timestep_current) const
108{
109 const double h_n = timestep_current.dt();
110 // Forced the computed time step size in the given range
111 // (see the formulas in the documentation of the class)
112 const double h_in_range = std::max(_h_min, std::min(h_new, _h_max));
113 // Limit the step size change ratio.
114 double limited_h =
115 std::max(_rel_h_min * h_n, std::min(h_in_range, _rel_h_max * h_n));
116
117 if (!previous_step_accepted)
118 {
119 // If the last time step was rejected and the new predicted time step
120 // size is identical to that of the previous rejected step, the new
121 // step size is then reduced by half.
122 if (std::abs(limited_h - timestep_current.dt()) <
123 std::numeric_limits<double>::min())
124 {
125 limited_h = std::max(_h_min, 0.5 * limited_h);
126 }
127
128 // If the last time step was rejected and the new predicted time step
129 // size is larger than the step size of the rejected step, the new step
130 // size takes the half of the size of the rejected step. This could
131 // happen when a time step is rejected due to a diverged non-linear
132 // solver. In such case, this algorithm may give a large time step size
133 // by using the diverged solution.
134 if (limited_h > timestep_current.dt())
135 {
136 limited_h = std::max(_h_min, 0.5 * timestep_current.dt());
137 }
138 }
139
140 if (_fixed_times_for_output.empty())
141 {
142 return limited_h;
143 }
144
145 // find first fixed output time larger than the current time, i.e.,
146 // current time < fixed output time
147 auto fixed_output_time_it = std::find_if(
149 [&timestep_current](auto const fixed_output_time)
150 { return timestep_current.current() < Time{fixed_output_time}; });
151
152 if (fixed_output_time_it != _fixed_times_for_output.end())
153 {
154 // check if the fixed output time is in the interval
155 // (current time, current time + limited_h)
156 if (Time{*fixed_output_time_it} <
157 timestep_current.current() + limited_h)
158 {
159 // check if the potential adjusted time step is larger than zero
160 if (std::abs(*fixed_output_time_it - timestep_current.current()()) >
161 std::numeric_limits<double>::epsilon() *
162 timestep_current.current()())
163 {
164 return *fixed_output_time_it - timestep_current.current()();
165 }
166 }
167 }
168 return limited_h;
169}
170
172 NumLib::TimeStep const& timestep_previous,
173 NumLib::TimeStep const& timestep_current) const
174{
175 return NumLib::canReduceTimestepSize(timestep_previous, timestep_current,
176 _h_min);
177}
178} // namespace NumLib
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:40
const double _h_max
maximum step size.
std::vector< double > const _fixed_times_for_output
std::tuple< bool, double > next(double solution_error, int number_iterations, NumLib::TimeStep &timestep_previous, NumLib::TimeStep &timestep_current) override
const double _h_min
minimum step size.
virtual bool canReduceTimestepSize(NumLib::TimeStep const &timestep_previous, NumLib::TimeStep const &timestep_current) const override
Query the timestepper if further time step size reduction is possible.
double limitStepSize(const double h_new, const bool previous_step_accepted, NumLib::TimeStep const &timestep_current) const
const double _h0
initial time step size.
Time step object.
Definition TimeStep.h:31
void setAccepted(bool const accepted)
Definition TimeStep.h:84
Time current() const
return current time step
Definition TimeStep.h:78
bool isAccepted() const
Definition TimeStep.h:85
double dt() const
time step size from _previous
Definition TimeStep.h:80
std::size_t timeStepNumber() const
the time step number
Definition TimeStep.h:82
bool canReduceTimestepSize(TimeStep const &timestep_previous, TimeStep const &timestep_current, double const min_dt)