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{
23double EvolutionaryPIDcontroller::next(double const solution_error,
24 int const /*number_iterations*/,
25 NumLib::TimeStep& /*timestep_previous*/,
26 NumLib::TimeStep& timestep_current)
27{
28 const bool is_previous_step_accepted = timestep_current.isAccepted();
29
30 const double e_n = solution_error;
31 const double zero_threshold = std::numeric_limits<double>::epsilon();
32 // step rejected.
33 if (e_n > _tol)
34 {
35 timestep_current.setAccepted(false);
36
37 double h_new = (e_n > zero_threshold)
38 ? timestep_current.dt() * _tol / e_n
39 : 0.5 * timestep_current.dt();
40
41 h_new =
42 limitStepSize(h_new, is_previous_step_accepted, timestep_current);
43
44 WARN(
45 "This step is rejected due to the relative change from the"
46 " solution of the previous\n"
47 "\t time step to the current solution exceeds the given tolerance"
48 " of {:g}.\n"
49 "\t This time step will be repeated with a new time step size of"
50 " {:g}\n"
51 "\t or the simulation will be halted.",
52 _tol, h_new);
53
54 return h_new;
55 }
56
57 // step accepted.
58 timestep_current.setAccepted(true);
59
60 if (timestep_current.timeStepNumber() == 0)
61 {
62 _e_n_minus1 = e_n;
63
64 return _h0;
65 }
66 const double h_n = timestep_current.dt();
67 double h_new = h_n;
68
69 if (e_n > zero_threshold)
70 {
71 if (_e_n_minus1 > zero_threshold)
72 {
73 if (_e_n_minus2 > zero_threshold)
74 {
75 h_new =
76 std::pow(_e_n_minus1 / e_n, _kP) *
77 std::pow(_tol / e_n, _kI) *
78 std::pow(_e_n_minus1 * _e_n_minus1 / (e_n * _e_n_minus2),
79 _kD) *
80 h_n;
81 }
82 else
83 {
84 h_new = std::pow(_e_n_minus1 / e_n, _kP) *
85 std::pow(_tol / e_n, _kI) * h_n;
86 }
87 }
88 else
89 {
90 h_new = std::pow(_tol / e_n, _kI) * h_n;
91 }
92 }
93
94 h_new = limitStepSize(h_new, is_previous_step_accepted, timestep_current);
95
97 _e_n_minus1 = e_n;
98
99 return h_new;
100}
101
103 const double h_new, const bool previous_step_accepted,
104 NumLib::TimeStep const& timestep_current) const
105{
106 const double h_n = timestep_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 - timestep_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 > timestep_current.dt())
132 {
133 limited_h = std::max(_h_min, 0.5 * timestep_current.dt());
134 }
135 }
136
137 if (_fixed_times_for_output.empty())
138 {
139 return limited_h;
140 }
141
142 // find first fixed output time larger than the current time, i.e.,
143 // current time < fixed output time
144 auto fixed_output_time_it = std::find_if(
146 [&timestep_current](auto const fixed_output_time)
147 { return timestep_current.current() < Time{fixed_output_time}; });
148
149 if (fixed_output_time_it != _fixed_times_for_output.end())
150 {
151 // check if the fixed output time is in the interval
152 // (current time, current time + limited_h)
153 if (Time{*fixed_output_time_it} <
154 timestep_current.current() + limited_h)
155 {
156 // check if the potential adjusted time step is larger than zero
157 if (std::abs(*fixed_output_time_it - timestep_current.current()()) >
158 std::numeric_limits<double>::epsilon() *
159 timestep_current.current()())
160 {
161 return *fixed_output_time_it - timestep_current.current()();
162 }
163 }
164 }
165 return limited_h;
166}
167
169 NumLib::TimeStep const& timestep_previous,
170 NumLib::TimeStep const& timestep_current) const
171{
172 return NumLib::canReduceTimestepSize(timestep_previous, timestep_current,
173 _h_min);
174}
175} // 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
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.
double next(double solution_error, int number_iterations, NumLib::TimeStep &timestep_previous, NumLib::TimeStep &timestep_current) override
Time step object.
Definition TimeStep.h:33
void setAccepted(bool const accepted)
Definition TimeStep.h:86
Time current() const
return current time step
Definition TimeStep.h:80
bool isAccepted() const
Definition TimeStep.h:87
double dt() const
time step size from _previous
Definition TimeStep.h:82
std::size_t timeStepNumber() const
the time step number
Definition TimeStep.h:84
bool canReduceTimestepSize(TimeStep const &timestep_previous, TimeStep const &timestep_current, double const min_dt)