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