OGS
IterationNumberBasedTimeStepping.cpp
Go to the documentation of this file.
1
14
15#include <algorithm>
16#include <cassert>
17#include <cmath>
18#include <limits>
19#include <utility>
20
21#include "BaseLib/Algorithm.h"
25
26namespace NumLib
27{
29 double const t_initial, double const t_end, double const min_dt,
30 double const max_dt, double const initial_dt,
31 MultiplyerInterpolationType const multiplier_interpolation_type,
32 std::vector<int>&& iter_times_vector,
33 std::vector<double>&& multiplier_vector,
34 std::vector<double> const& fixed_times_for_output)
35 : TimeStepAlgorithm(t_initial, t_end),
36 _iter_times_vector(std::move(iter_times_vector)),
37 _multiplier_vector(std::move(multiplier_vector)),
38 _min_dt(min_dt),
39 _max_dt(max_dt),
40 _initial_dt(initial_dt),
41 _multiplier_interpolation_type(multiplier_interpolation_type),
42 _max_iter(_iter_times_vector.empty() ? 0 : _iter_times_vector.back()),
43 _fixed_times_for_output(fixed_times_for_output)
44{
45 if (_iter_times_vector.empty())
46 {
47 OGS_FATAL("Vector of iteration numbers must not be empty.");
48 }
49 if (_iter_times_vector.size() != _multiplier_vector.size())
50 {
52 "Vector of iteration numbers must be of the same size as the "
53 "vector of multipliers.");
54 }
55 if (!std::is_sorted(std::begin(_iter_times_vector),
56 std::end(_iter_times_vector)))
57 {
58 OGS_FATAL("Vector of iteration numbers must be sorted.");
59 }
60}
61
62double IterationNumberBasedTimeStepping::next(double const /*solution_error*/,
63 int const number_iterations,
64 NumLib::TimeStep& ts_previous,
65 NumLib::TimeStep& ts_current)
66{
67 _iter_times = number_iterations;
68
69 if (ts_previous.isAccepted())
70 {
71 ts_previous = ts_current;
72 }
73
74 // confirm current time and move to the next if accepted
75 if (ts_current.isAccepted())
76 {
77 ts_previous.setAccepted(true);
78 return getNextTimeStepSize(ts_previous, ts_current);
79 }
80
81 double dt = getNextTimeStepSize(ts_previous, ts_current);
82 // In case it is the first time be rejected, re-computed dt again with
83 // current dt
84 if (std::abs(dt - ts_current.dt()) < std::numeric_limits<double>::epsilon())
85 {
86 // time step was rejected, keep dt for the next dt computation.
87 ts_previous = // essentially equal to _ts_prev.dt = _ts_current.dt.
88 TimeStep{ts_previous.previous(), ts_previous.previous() + dt,
89 ts_previous.timeStepNumber()};
90 dt = getNextTimeStepSize(ts_previous, ts_current);
91 }
92
93 // time step was rejected, keep dt for the next dt computation.
94 ts_previous = // essentially equal to ts_previous.dt = _ts_current.dt.
95 TimeStep{ts_previous.previous(), ts_previous.previous() + dt,
96 ts_previous.timeStepNumber()};
97 ts_current = TimeStep{ts_current.previous(), ts_current.previous() + dt,
98 ts_current.timeStepNumber()};
99
100 return dt;
101}
102
104 int const number_iterations, bool const current_time_step_is_accepted,
105 std::vector<int> const& nonlinear_iteration_numbers,
106 std::vector<double> const& multipliers,
107 MultiplyerInterpolationType const multiplier_interpolation_type)
108{
109 double multiplier = multipliers.front();
110 switch (multiplier_interpolation_type)
111 {
113 {
114 auto const& pwli = MathLib::PiecewiseLinearInterpolation(
115 nonlinear_iteration_numbers, multipliers, false);
116 multiplier = pwli.getValue(number_iterations);
117 DBUG("Using piecewise linear iteration-based time stepping.");
118 break;
119 }
121 DBUG("Using piecewise constant iteration-based time stepping.");
122 auto const& piecewise_constant_interpolation =
124 nonlinear_iteration_numbers, multipliers);
125 multiplier =
126 piecewise_constant_interpolation.value(number_iterations);
127 break;
128 }
129
130 if (!current_time_step_is_accepted && (multiplier >= 1.0))
131 {
132 return *std::min_element(multipliers.begin(), multipliers.end());
133 }
134
135 return multiplier;
136}
137
139 NumLib::TimeStep const& ts_previous,
140 NumLib::TimeStep const& ts_current) const
141{
142 double dt = 0.0;
143
144 // In first time step and first non-linear iteration take the initial dt.
145 if (ts_previous.timeStepNumber() == 0 && _iter_times == 0)
146 {
147 dt = _initial_dt;
148 }
149 else
150 {
151 // Attention: for the first time step and second iteration the
152 // ts_prev.dt is 0 and 0*multiplier is the next dt, which will be
153 // clamped to the minimum dt.
154 dt = ts_previous.dt() *
158 }
159
160 if (_fixed_times_for_output.empty())
161 {
162 return std::clamp(dt, _min_dt, _max_dt);
163 }
164
165 // restrict dt to _max_dt before taking fixed times for output into account
166 dt = std::min(dt, _max_dt);
167
168 // find first fixed timestep for output larger than the current time, i.e.,
169 // current time < fixed output time
170 auto fixed_output_time_it = std::find_if(
172 [&ts_current](auto const fixed_output_time)
173 { return ts_current.current()() < fixed_output_time; });
174
175 if (fixed_output_time_it != _fixed_times_for_output.end())
176 {
177 // check if the fixed output time is in the interval
178 // (current time, current time + dt)
179 if (*fixed_output_time_it < ts_current.current()() + dt)
180 {
181 // check if the potential adjusted time step is larger than zero
182 if (std::abs(*fixed_output_time_it - ts_current.current()()) >
183 std::numeric_limits<double>::epsilon() * ts_current.current()())
184 {
185 return *fixed_output_time_it - ts_current.current()();
186 }
187 }
188 }
189 return std::clamp(dt, _min_dt, _max_dt);
190}
191
193 NumLib::TimeStep const& timestep_previous,
194 NumLib::TimeStep const& timestep_current) const
195{
196 return NumLib::canReduceTimestepSize(timestep_previous, timestep_current,
197 _min_dt);
198}
199
200} // namespace NumLib
#define OGS_FATAL(...)
Definition Error.h:26
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:30
Definition of the PiecewiseConstantInterpolation class.
Definition of the PiecewiseLinearInterpolation class.
IterationNumberBasedTimeStepping(double const t_initial, double const t_end, double const min_dt, double const max_dt, double const initial_dt, MultiplyerInterpolationType const multiplier_interpolation_type, std::vector< int > &&iter_times_vector, std::vector< double > &&multiplier_vector, std::vector< double > const &fixed_times_for_output)
double next(double solution_error, int number_iterations, NumLib::TimeStep &ts_previous, NumLib::TimeStep &ts_current) override
const MultiplyerInterpolationType _multiplier_interpolation_type
Interpolation type for the multiplier.
double getNextTimeStepSize(NumLib::TimeStep const &ts_previous, NumLib::TimeStep const &ts_current) const
Calculate the next time step size.
const double _max_dt
The maximum allowed time step size.
const std::vector< double > _multiplier_vector
This vector stores the multiplier coefficients.
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.
int _iter_times
The number of nonlinear iterations.
const double _min_dt
The minimum allowed time step size.
Interface of time stepping algorithms.
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
Time previous() const
return previous time step
Definition TimeStep.h:78
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
double findMultiplier(int const number_iterations, bool const current_time_step_is_accepted, std::vector< int > const &nonlinear_iteration_numbers, std::vector< double > const &multipliers, MultiplyerInterpolationType const multiplier_interpolation_type)
Find a multiplier for the given number of iterations.
bool canReduceTimestepSize(TimeStep const &timestep_previous, TimeStep const &timestep_current, double const min_dt)