OGS
FixedTimeStepping.cpp
Go to the documentation of this file.
1
13#include "FixedTimeStepping.h"
14
15#include <spdlog/fmt/bundled/ranges.h>
16
17#include <algorithm>
18#include <cassert>
19#include <limits>
20#include <numeric>
21
23
24namespace
25{
27NumLib::Time addTimeIncrement(std::vector<double>& delta_ts,
28 std::size_t const repeat, double const delta_t,
29 NumLib::Time const t_curr)
30{
31 auto const new_size = delta_ts.size() + repeat;
32 try
33 {
34 delta_ts.resize(new_size, delta_t);
35 }
36 catch (std::exception const& e)
37 {
39 "Resize of the time steps vector failed for the requested new size "
40 "{:d}. Probably there is not enough memory ({:g} GiB "
41 "requested).\nThrown exception: {:s}",
42 new_size,
43 new_size * sizeof(double) / 1024. / 1024. / 1024.,
44 e.what());
45 }
46
47 // Multiplying dt * repeat is not the same as in the current
48 // implementation of the time loop, where the dt's are added.
49 // Therefore the sum of all dt is taken here.
50 return std::accumulate(end(delta_ts) - repeat, end(delta_ts), t_curr);
51}
52} // namespace
53
54namespace NumLib
55{
56std::size_t findDeltatInterval(Time const& t_initial,
57 std::vector<double> const& delta_ts,
58 Time const& fixed_output_time)
59{
60 if (fixed_output_time < t_initial)
61 {
62 return std::numeric_limits<std::size_t>::max();
63 }
64
65 auto timestepper_time = t_initial;
66 for (std::size_t k = 0; k < delta_ts.size(); ++k)
67 {
68 if (timestepper_time <= fixed_output_time &&
69 fixed_output_time < timestepper_time + delta_ts[k])
70 {
71 return k;
72 }
73 timestepper_time += delta_ts[k];
74 }
75
76 return std::numeric_limits<std::size_t>::max();
77}
78
80 NumLib::Time const t_initial, NumLib::Time const t_end,
81 std::vector<double>& delta_ts,
82 std::vector<double> const& fixed_times_for_output)
83{
84 if (fixed_times_for_output.empty())
85 {
86 return;
87 }
88
89 if (auto lower_bound = std::lower_bound(
90 begin(fixed_times_for_output), end(fixed_times_for_output),
91 t_initial,
92 [](auto const time, NumLib::Time const& initial_time)
93 { return NumLib::Time(time) < initial_time; });
94 lower_bound != begin(fixed_times_for_output))
95 {
96 WARN(
97 "Request for output at times {}, but the simulation's start time "
98 "is {}. Output will be skipped.",
99 fmt::join(begin(fixed_times_for_output), lower_bound, ", "),
100 t_initial());
101 }
102
103 if (auto upper_bound =
104 std::upper_bound(begin(fixed_times_for_output),
105 end(fixed_times_for_output), t_end());
106 upper_bound != end(fixed_times_for_output))
107 {
108 WARN(
109 "Request for output at times {}, but simulation's end time is {}. "
110 "Output will be skipped.",
111 fmt::join(upper_bound, end(fixed_times_for_output), ", "),
112 t_end());
113 }
114
115 if (delta_ts.empty())
116 {
117 WARN("No timesteps specified.");
118 return;
119 }
120
121 // incorporate fixed output times into dts vector
122 for (auto const fixed_time_for_output : fixed_times_for_output)
123 {
124 auto const interval_number = findDeltatInterval(
125 t_initial, delta_ts, Time(fixed_time_for_output));
126 if (interval_number == std::numeric_limits<std::size_t>::max())
127 {
128 WARN("Did not find interval for fixed output time {}",
129 fixed_time_for_output);
130 continue;
131 }
132
133 auto const lower_bound = std::accumulate(
134 begin(delta_ts), begin(delta_ts) + interval_number, t_initial);
135 auto const upper_bound = lower_bound + delta_ts[interval_number];
136 // in order to use the comparison from struct Time
137 if (NumLib::Time(fixed_time_for_output) == lower_bound)
138 {
139 continue;
140 }
141 if (upper_bound == NumLib::Time(fixed_time_for_output))
142 {
143 continue;
144 }
145 delta_ts[interval_number] = fixed_time_for_output - lower_bound();
146
147 delta_ts.insert(delta_ts.begin() + interval_number + 1,
148 upper_bound() - fixed_time_for_output);
149 }
150}
151
153 double t0, double tn, std::vector<RepeatDtPair> const& repeat_dt_pairs,
154 std::vector<double> const& fixed_times_for_output)
155 : TimeStepAlgorithm(t0, tn)
156{
157 Time t_curr = _t_initial;
158
159 if (!areRepeatDtPairsValid(repeat_dt_pairs))
160 {
161 OGS_FATAL("FixedTimeStepping: Couldn't construct object from data");
162 }
163 for (auto const& [repeat, delta_t] : repeat_dt_pairs)
164 {
165 if (t_curr <= _t_end)
166 {
167 t_curr = addTimeIncrement(_dt_vector, repeat, delta_t, t_curr);
168 }
169 }
170
171 // append last delta_t until t_end is reached
172 if (t_curr <= _t_end)
173 {
174 auto const delta_t = std::get<1>(repeat_dt_pairs.back());
175 auto const repeat = static_cast<std::size_t>(
176 std::ceil((_t_end() - t_curr()) / delta_t));
177 addTimeIncrement(_dt_vector, repeat, delta_t, t_curr);
178 }
179
181 fixed_times_for_output);
182}
183
184FixedTimeStepping::FixedTimeStepping(double t0, double t_end, double dt)
185 : TimeStepAlgorithm(t0, t_end)
186{
187 auto const new_size =
188 static_cast<std::size_t>(std::ceil((t_end - t0) / dt));
189 try
190 {
191 _dt_vector = std::vector<double>(new_size, dt);
192 }
193 catch (std::length_error const& e)
194 {
195 OGS_FATAL(
196 "Resize of the time steps vector failed for the requested new "
197 "size {}. Probably there is not enough memory ({:g} GiB "
198 "requested).\n"
199 "Thrown exception: {}",
200 new_size, new_size * sizeof(double) / 1024. / 1024. / 1024.,
201 e.what());
202 }
203 catch (std::bad_alloc const& e)
204 {
205 OGS_FATAL(
206 "Allocation of the time steps vector failed for the requested "
207 "size {}. Probably there is not enough memory ({:g} GiB "
208 "requested).\n"
209 "Thrown exception: {}",
210 new_size,
211 new_size * sizeof(double) / 1024. / 1024. / 1024.,
212 e.what());
213 }
214}
215
216std::tuple<bool, double> FixedTimeStepping::next(
217 double const /*solution_error*/, int const /*number_iterations*/,
218 NumLib::TimeStep& /*ts_previous*/, NumLib::TimeStep& ts_current)
219{
220 // check if last time step
221 if (ts_current.timeStepNumber() == _dt_vector.size() ||
222 ts_current.current() >= end())
223 {
224 return std::make_tuple(true, 0.0);
225 }
226
227 double dt = _dt_vector[ts_current.timeStepNumber()];
228 if (ts_current.current() + dt > end())
229 { // upper bound by t_end
230 dt = end()() - ts_current.current()();
231 }
232
233 return std::make_tuple(true, dt);
234}
235
237 std::vector<RepeatDtPair> const& repeat_dt_pairs)
238{
239 if (repeat_dt_pairs.empty())
240 {
241 return false;
242 }
243
244 for (auto const& [repeat, delta_t] : repeat_dt_pairs)
245 {
246 if (repeat == 0)
247 {
248 ERR("FixedTimeStepping: <repeat> is zero.");
249 return false;
250 }
251 if (delta_t <= 0.0)
252 {
253 ERR("FixedTimeStepping: timestep <delta_t> is <= 0.0.");
254 return false;
255 }
256 }
257 return true;
258}
259
260} // namespace NumLib
#define OGS_FATAL(...)
Definition Error.h:26
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:45
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:40
std::tuple< bool, double > next(double solution_error, int number_iterations, NumLib::TimeStep &ts_previous, NumLib::TimeStep &ts_current) override
static bool areRepeatDtPairsValid(std::vector< RepeatDtPair > const &repeat_dt_pairs)
std::vector< double > _dt_vector
a vector of time step sizes
FixedTimeStepping(double t0, double t_end, double dt)
Interface of time stepping algorithms.
const Time _t_initial
initial time
Time end() const
return the end of time steps
Time step object.
Definition TimeStep.h:31
Time current() const
return current time step
Definition TimeStep.h:78
std::size_t timeStepNumber() const
the time step number
Definition TimeStep.h:82
void incorporateFixedTimesForOutput(NumLib::Time const t_initial, NumLib::Time const t_end, std::vector< double > &delta_ts, std::vector< double > const &fixed_times_for_output)
std::size_t findDeltatInterval(Time const &t_initial, std::vector< double > const &delta_ts, Time const &fixed_output_time)
NumLib::Time addTimeIncrement(std::vector< double > &delta_ts, std::size_t const repeat, double const delta_t, NumLib::Time const t_curr)
Returns sum of the newly added time increments.