GCC Code Coverage Report


Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/common/timeloop.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 172 187 92.0%
Functions: 29 38 76.3%
Branches: 243 456 53.3%

Line Branch Exec Source
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 //
4 // SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5 // SPDX-License-Identifier: GPL-3.0-or-later
6 //
7 /*!
8 * \file
9 * \ingroup Core
10 * \brief Manages the handling of time dependent problems
11 */
12 #ifndef DUMUX_TIME_LOOP_HH
13 #define DUMUX_TIME_LOOP_HH
14
15 #include <algorithm>
16 #include <queue>
17 #include <iomanip>
18 #include <chrono>
19 #include <type_traits>
20 #include <initializer_list>
21 #include <optional>
22
23 #include <dune/common/float_cmp.hh>
24 #include <dune/common/timer.hh>
25
26 #include <dune/common/parallel/communication.hh>
27
28 #include <dune/common/parallel/mpihelper.hh>
29 #include <dune/common/exceptions.hh>
30
31 #include <dumux/common/typetraits/typetraits.hh>
32 #include <dumux/common/parameters.hh>
33 #include <dumux/io/format.hh>
34
35 namespace Dumux {
36
37
38 #ifndef DOXYGEN
39 namespace Detail::TimeLoop {
40
41 template<typename Scalar, typename R, typename P>
42 Scalar toSeconds(std::chrono::duration<R, P> duration)
43 {
44 using Second = std::chrono::duration<Scalar, std::ratio<1>>;
45 106 return std::chrono::duration_cast<Second>(duration).count();
46 }
47
48 template<typename Scalar, typename T, std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
49 Scalar toSeconds(T duration)
50 { return static_cast<Scalar>(duration); }
51
52 // overload for nice static_assert message in case of unuspported type
53 template<typename Scalar, typename T, std::enable_if_t<!std::is_floating_point_v<T>, bool> = true>
54 Scalar toSeconds(T duration)
55 { static_assert(AlwaysFalse<T>::value, "Given type not supported for representation of time values"); }
56
57 } // namespace Detail::TimeLoop
58 #endif // DOXYGEN
59
60 /*!
61 * \ingroup Core
62 * \brief Manages the handling of time dependent problems.
63 *
64 * This class facilitates the time management of the simulation.
65 * It doesn't manage any user data, but keeps track of what the
66 * current time, time step size and "episode" of the
67 * simulation is. It triggers the initialization of the problem and
68 * is responsible for the time control of a simulation run.
69 *
70 * The time manager allows to specify a sequence of "episodes" which
71 * determine the boundary conditions of a problem. This approach
72 * is handy if the problem is not static, i.e. the boundary
73 * conditions change over time.
74 *
75 * An episode is a span of simulated time in which
76 * the problem behaves in a specific way. It is characterized by
77 * the (simulation) time it starts, its length and a consecutive
78 * index starting at 0.
79 *
80 * \note Time and time step sizes are in units of seconds
81 */
82 template<class S>
83 class TimeLoopBase
84 {
85 public:
86 using Scalar = S;
87
88 //! Abstract base class needs virtual constructor
89 virtual ~TimeLoopBase() {};
90
91 /*!
92 * \brief Return the time \f$\mathrm{[s]}\f$ before the time integration.
93 * To get the time after the time integration you have to add timeStepSize() to
94 * time().
95 */
96 virtual Scalar time() const = 0;
97
98 /*!
99 * \brief Returns the suggested time step length \f$\mathrm{[s]}\f$
100 */
101 virtual Scalar timeStepSize() const = 0;
102
103 /*!
104 * \brief Get the maximum possible time step size \f$\mathrm{[s]}\f$
105 */
106 virtual Scalar maxTimeStepSize() const = 0;
107
108 /*!
109 * \brief Advance to the next time step.
110 */
111 virtual void advanceTimeStep() = 0;
112
113 /*!
114 * \brief Set the current time step size to a given value.
115 * \param dt The new value for the time step size \f$\mathrm{[s]}\f$
116 */
117 virtual void setTimeStepSize(Scalar dt) = 0;
118
119 /*!
120 * \brief Set the current time step size to a given value.
121 * \param dt The new value for the time step size \f$\mathrm{[s]}\f$
122 */
123 template<class Rep, class Period>
124 void setTimeStepSize(std::chrono::duration<Rep, Period> dt)
125
0/8
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
8 { setTimeStepSize(Detail::TimeLoop::toSeconds<Scalar>(dt)); }
126
127 /*!
128 * \brief Returns true if the simulation is finished.
129 */
130 virtual bool finished() const = 0;
131 };
132
133 /*!
134 * \ingroup Core
135 * \brief The default time loop for instationary simulations
136 */
137 template <class Scalar>
138
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
10 class TimeLoop : public TimeLoopBase<Scalar>
139 {
140 public:
141 356 TimeLoop(Scalar startTime, Scalar dt, Scalar tEnd, bool verbose = true)
142 356 : timer_(false)
143 {
144 356 reset(startTime, dt, tEnd, verbose);
145 356 }
146
147 template<class Rep1, class Period1,
148 class Rep2, class Period2,
149 class Rep3, class Period3>
150 8 TimeLoop(std::chrono::duration<Rep1, Period1> startTime,
151 std::chrono::duration<Rep2, Period2> dt,
152 std::chrono::duration<Rep3, Period3> tEnd,
153 bool verbose = true)
154 : TimeLoop(
155 Detail::TimeLoop::toSeconds<Scalar>(startTime),
156 Detail::TimeLoop::toSeconds<Scalar>(dt),
157 Detail::TimeLoop::toSeconds<Scalar>(tEnd),
158 verbose
159
8/16
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 1 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 1 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 1 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 1 times.
✗ Branch 26 not taken.
26 ){}
160
161 /*!
162 * \name Simulated time and time step management
163 * @{
164 */
165
166 /*!
167 * \brief Tells the time loop to start tracking the time.
168 */
169 void start()
170 {
171
6/12
✓ Branch 0 taken 321 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 4 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 4 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 4 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 4 times.
✗ Branch 15 not taken.
349 timer_.start();
172 }
173
174 /*!
175 * \brief Tells the time loop to stop tracking the time.
176 * \return the wall clock time (CPU time) spent until now
177 */
178 double stop()
179 {
180 4 return timer_.stop();
181 }
182
183 /*!
184 * \brief Reset the timer
185 */
186 void resetTimer()
187 {
188
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
4 timer_.reset();
189 }
190
191 /*!
192 * \brief Reset the time loop
193 */
194 template<class Rep1, class Period1,
195 class Rep2, class Period2,
196 class Rep3, class Period3>
197 void reset(std::chrono::duration<Rep1, Period1> startTime,
198 std::chrono::duration<Rep2, Period2> dt,
199 std::chrono::duration<Rep3, Period3> tEnd,
200 bool verbose = true)
201 {
202
0/16
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
16 reset(
203 Detail::TimeLoop::toSeconds<Scalar>(startTime),
204 Detail::TimeLoop::toSeconds<Scalar>(dt),
205 Detail::TimeLoop::toSeconds<Scalar>(tEnd)
206 );
207 }
208
209 /*!
210 * \brief Reset the time loop
211 */
212 360 void reset(Scalar startTime, Scalar dt, Scalar tEnd, bool verbose = true)
213 {
214 360 verbose_ =
215
3/4
✓ Branch 0 taken 360 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 19 times.
✓ Branch 3 taken 341 times.
720 verbose &&
216
2/2
✓ Branch 1 taken 19 times.
✓ Branch 2 taken 341 times.
360 Dune::MPIHelper::getCommunication().rank() == 0;
217
218 360 startTime_ = startTime;
219 360 time_ = startTime;
220 360 endTime_ = tEnd;
221
222 360 previousTimeStepSize_ = 0.0;
223 360 userSetMaxTimeStepSize_ = std::numeric_limits<Scalar>::max();
224 360 timeStepIdx_ = 0;
225 360 finished_ = false;
226 360 timeAfterLastTimeStep_ = 0.0;
227 360 timeStepWallClockTime_ = 0.0;
228
229 // ensure that dt is not greater than tEnd-startTime
230 360 setTimeStepSize(dt);
231
232 360 timer_.stop();
233 360 timer_.reset();
234 360 }
235
236 /*!
237 * \brief Advance time step.
238 */
239 20377 void advanceTimeStep() override
240 {
241 20377 timeStepIdx_++;
242 20377 time_ += timeStepSize_;
243 20377 previousTimeStepSize_ = timeStepSize_;
244
245 // compute how long the last time step took
246 40754 const auto cpuTime = wallClockTime();
247 20377 timeStepWallClockTime_ = cpuTime - timeAfterLastTimeStep_;
248 20377 timeAfterLastTimeStep_ = cpuTime;
249
250 // ensure that using current dt we don't exceed tEnd in next time step
251 20377 setTimeStepSize(timeStepSize_);
252 20377 }
253
254 /*!
255 * \brief Set the current simulated time, don't change the current
256 * time step index.
257 *
258 * \param t The time \f$\mathrm{[s]}\f$ which should be jumped to
259 */
260 template<class ScalarOrDuration>
261 void setTime(ScalarOrDuration t)
262
0/8
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
8 { time_ = Detail::TimeLoop::toSeconds<Scalar>(t); }
263
264 /*!
265 * \brief Set the current simulated time and the time step index.
266 *
267 * \param t The time \f$\mathrm{[s]}\f$ which should be jumped to
268 * \param stepIdx The new time step index
269 */
270 template<class ScalarOrDuration>
271 void setTime(ScalarOrDuration t, int stepIdx)
272 {
273 time_ = Detail::TimeLoop::toSeconds<Scalar>(t);
274 timeStepIdx_ = stepIdx;
275 }
276
277 /*!
278 * \brief Return the time \f$\mathrm{[s]}\f$ before the time integration.
279 * To get the time after the time integration you have to add timeStepSize() to
280 * time().
281 */
282 Scalar time() const final
283
17/27
✓ Branch 0 taken 98625452 times.
✓ Branch 1 taken 4797413 times.
✓ Branch 2 taken 421619 times.
✓ Branch 3 taken 728330 times.
✓ Branch 4 taken 99708 times.
✓ Branch 5 taken 87191 times.
✓ Branch 6 taken 40 times.
✓ Branch 7 taken 266 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 935 times.
✓ Branch 10 taken 137 times.
✓ Branch 11 taken 4 times.
✓ Branch 12 taken 9 times.
✓ Branch 13 taken 61 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 10 times.
✗ Branch 17 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✓ Branch 24 taken 2 times.
✓ Branch 25 taken 64 times.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✓ Branch 30 taken 4 times.
104761742 { return time_; }
284
285 /*!
286 * \brief Returns the number of (simulated) seconds which the simulation runs.
287 */
288 Scalar endTime() const
289 { return endTime_; }
290
291 /*!
292 * \brief Set the time of simulated seconds at which the simulation runs.
293 *
294 * \param t The time \f$\mathrm{[s]}\f$ at which the simulation is finished
295 */
296 4 void setEndTime(Scalar t)
297 {
298 4 endTime_ = t;
299
1/2
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
4 if (verbose_)
300
2/6
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
8 std::cout << Fmt::format("Set new end time to t = {:.5g} seconds.\n", t);
301 4 }
302
303 /*!
304 * \brief Returns the current wall clock time (cpu time) spend in this time loop
305 */
306 double wallClockTime() const
307 20377 { return timer_.elapsed(); }
308
309 using TimeLoopBase<Scalar>::setTimeStepSize;
310 /*!
311 * \brief Set the current time step size to a given value.
312 *
313 * If the step size would exceed the length of the current
314 * episode, the timeStep() method will take care that the step
315 * size won't exceed the episode or the end of the simulation,
316 * though.
317 *
318 * \param dt The new value for the time step size \f$\mathrm{[s]}\f$
319 */
320 57348 void setTimeStepSize(Scalar dt) final
321 {
322 using std::min;
323
2/2
✓ Branch 1 taken 10749 times.
✓ Branch 2 taken 46599 times.
57348 timeStepSize_ = min(dt, maxTimeStepSize());
324 // Warn if dt is so small w.r.t. current time that it renders float addition meaningless
325 // For instance, consider (may depend on architecture):
326 // double cien = 100;
327 // double mil = 1000;
328 // if (cien + 1e-14 == cien) std::cout << "Will not be printed" << std::endl;
329 // if (mil + 1e-14 == mil) std::cout << "Will be printed" << std::endl;
330
4/4
✓ Branch 1 taken 56690 times.
✓ Branch 2 taken 658 times.
✓ Branch 3 taken 4 times.
✓ Branch 4 taken 56686 times.
57348 if (!finished() && (time_ + timeStepSize_ == time_))
331 4 std::cerr << Fmt::format("You have set a very small timestep size (dt = {:.5g}).", timeStepSize_)
332
3/8
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 4 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 4 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
8 << " This might lead to numerical problems!\n";
333 57348 }
334
335 /*!
336 * \brief Set the maximum time step size to a given value.
337 *
338 * \param maxDt The new value for the maximum time step size \f$\mathrm{[s]}\f$
339 * \note This may also reduce the currently set timestep size if needed to comply with the set maximum
340 */
341 template<class ScalarOrDuration>
342 void setMaxTimeStepSize(ScalarOrDuration maxDt)
343 {
344
2/8
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
299 userSetMaxTimeStepSize_ = Detail::TimeLoop::toSeconds<Scalar>(maxDt);
345
2/8
✓ Branch 1 taken 294 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
991 setTimeStepSize(timeStepSize_);
346 }
347
348 /*!
349 * \brief Returns the suggested time step length \f$\mathrm{[s]}\f$ so that we
350 * don't miss the beginning of the next episode or cross
351 * the end of the simulation.
352 */
353 692442197 Scalar timeStepSize() const final
354
26/58
✓ Branch 0 taken 63 times.
✓ Branch 1 taken 16074 times.
✓ Branch 2 taken 54 times.
✓ Branch 3 taken 22 times.
✓ Branch 4 taken 3616 times.
✓ Branch 5 taken 53 times.
✓ Branch 6 taken 4 times.
✓ Branch 7 taken 479 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 5 times.
✓ Branch 10 taken 44 times.
✓ Branch 11 taken 2 times.
✓ Branch 12 taken 3 times.
✓ Branch 13 taken 10 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✓ Branch 25 taken 4 times.
✗ Branch 26 not taken.
✓ Branch 27 taken 4 times.
✓ Branch 28 taken 22 times.
✓ Branch 29 taken 666 times.
✓ Branch 30 taken 22 times.
✓ Branch 31 taken 666 times.
✗ Branch 32 not taken.
✓ Branch 33 taken 4 times.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✓ Branch 38 taken 4 times.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
✓ Branch 43 taken 4 times.
✗ Branch 45 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
✓ Branch 48 taken 4 times.
✗ Branch 50 not taken.
✗ Branch 51 not taken.
✓ Branch 52 taken 66 times.
✗ Branch 53 not taken.
✗ Branch 55 not taken.
✗ Branch 56 not taken.
✗ Branch 57 not taken.
✓ Branch 58 taken 4 times.
✗ Branch 60 not taken.
✗ Branch 61 not taken.
✗ Branch 62 not taken.
✓ Branch 63 taken 4 times.
29573902 { return timeStepSize_; }
355
356 /*!
357 * \brief Returns number of time steps which have been
358 * executed since the beginning of the simulation.
359 */
360 int timeStepIndex() const
361 { return timeStepIdx_; }
362
363 /*!
364 * \brief The previous time step size
365 */
366 Scalar previousTimeStepSize() const
367 { return previousTimeStepSize_; }
368
369 /*!
370 * \brief Specify whether the simulation is finished
371 *
372 * \param finished If true the simulation is considered finished
373 * before the end time is reached, else it is only
374 * considered finished if the end time is reached.
375 */
376 void setFinished(bool finished = true)
377 { finished_ = finished; }
378
379 /*!
380 * \brief Returns true if the simulation is finished.
381 *
382 * This is the case if either setFinished(true) has been called or
383 * if the end time is reached.
384 */
385 157443 bool finished() const override
386 {
387
9/12
✓ Branch 0 taken 157443 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 155627 times.
✓ Branch 3 taken 1816 times.
✓ Branch 4 taken 692 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 688 times.
✓ Branch 7 taken 4 times.
✓ Branch 8 taken 70 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 66 times.
✓ Branch 11 taken 4 times.
158205 return finished_ || (endTime_ - time_) < baseEps_*(time_ - startTime_);
388 }
389
390 /*!
391 * \brief Returns true if the simulation is finished after the
392 * time level is incremented by the current time step size.
393 */
394 bool willBeFinished() const
395 {
396
4/4
✓ Branch 0 taken 14393 times.
✓ Branch 1 taken 123 times.
✓ Branch 2 taken 14286 times.
✓ Branch 3 taken 107 times.
14516 return finished() || (endTime_ - time_ - timeStepSize_) < baseEps_*timeStepSize_;
397 }
398
399 /*!
400 * \brief The current maximum time step size
401 * \note This gets aligned on every setTimeStepSize call to end time
402 * and other possible check points
403 */
404 57353 Scalar maxTimeStepSize() const override
405 {
406
2/2
✓ Branch 1 taken 56695 times.
✓ Branch 2 taken 658 times.
57353 if (finished())
407 return 0.0;
408
409 using std::min; using std::max;
410
3/4
✓ Branch 0 taken 56695 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 20057 times.
✓ Branch 3 taken 36638 times.
133447 return min(userSetMaxTimeStepSize_, max<Scalar>(0.0, endTime_ - time_));
411 }
412
413 /*!
414 * \brief State info on cpu time.
415 * \note Always call this after TimeLoop::advanceTimeStep()
416 */
417 20377 void reportTimeStep() const
418 {
419
2/2
✓ Branch 0 taken 19654 times.
✓ Branch 1 taken 723 times.
20377 if (verbose_)
420 {
421 19654 const auto cpuTime = wallClockTime();
422 using std::round;
423 19654 const auto percent = round( (time_ - startTime_) / (endTime_ - startTime_) * 100 );
424 std::cout << Fmt::format("[{:3.0f}%] ", percent)
425 19654 << Fmt::format("Time step {} done in {:.2g} seconds. ", timeStepIdx_, timeStepWallClockTime_)
426
8/22
✓ Branch 3 taken 19654 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 19654 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 19654 times.
✗ Branch 11 not taken.
✓ Branch 14 taken 19654 times.
✗ Branch 15 not taken.
✓ Branch 17 taken 19654 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 19654 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 19654 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✓ Branch 24 taken 19654 times.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
58962 << Fmt::format("Wall clock time: {:.5g}, time: {:.5g}, time step size: {:.5g}\n", cpuTime, time_, previousTimeStepSize_);
427 }
428 20377 }
429
430 /*!
431 * \brief Print final status and stops tracking the time.
432 */
433 template< class Communicator = Dune::Communication<typename Dune::MPIHelper::MPICommunicator> >
434
3/4
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
328 void finalize(const Communicator& comm = Dune::MPIHelper::getCommunication())
435 {
436 328 auto cpuTime = timer_.stop();
437
438
2/2
✓ Branch 0 taken 310 times.
✓ Branch 1 taken 18 times.
328 if (verbose_)
439
4/7
✓ Branch 3 taken 31 times.
✓ Branch 4 taken 279 times.
✓ Branch 5 taken 31 times.
✓ Branch 6 taken 279 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
899 std::cout << Fmt::format("Simulation took {:.5g} seconds on {} processes.\n", cpuTime, comm.size());
440
441
4/4
✓ Branch 0 taken 36 times.
✓ Branch 1 taken 261 times.
✓ Branch 2 taken 36 times.
✓ Branch 3 taken 261 times.
625 if (comm.size() > 1)
442 36 cpuTime = comm.sum(cpuTime);
443
444
2/2
✓ Branch 0 taken 310 times.
✓ Branch 1 taken 18 times.
328 if (verbose_)
445
2/6
✓ Branch 3 taken 310 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 310 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
620 std::cout << Fmt::format("The cumulative CPU time was {:.5g} seconds.\n", cpuTime);
446 328 }
447
448 //! If the time loop has verbose output
449 bool verbose() const
450 { return verbose_; }
451
452 //! Sets time loop verbosity
453 void setVerbose(bool verbose = true)
454 { verbose_ = verbose; }
455
456 /*
457 * @}
458 */
459
460 protected:
461 static constexpr Scalar baseEps_ = 1e-10;
462
463 Dune::Timer timer_;
464 Scalar time_;
465 Scalar endTime_;
466 Scalar startTime_;
467
468 Scalar timeStepSize_;
469 Scalar previousTimeStepSize_;
470 Scalar userSetMaxTimeStepSize_;
471 Scalar timeAfterLastTimeStep_, timeStepWallClockTime_;
472 int timeStepIdx_;
473 bool finished_;
474 bool verbose_;
475 };
476
477 // always fall back to floating-point representation
478 template<class Rep1, class Period1,
479 class Rep2, class Period2,
480 class Rep3, class Period3>
481 TimeLoop(std::chrono::duration<Rep1, Period1>,
482 std::chrono::duration<Rep2, Period2>,
483 std::chrono::duration<Rep3, Period3>,
484 bool verbose = true) -> TimeLoop<std::conditional_t<
485 std::is_floating_point_v<std::common_type_t<Rep1, Rep2, Rep3>>,
486 std::common_type_t<Rep1, Rep2, Rep3>,
487 double
488 >>;
489
490 /*!
491 * \ingroup Core
492 * \brief A time loop with a check point mechanism
493 */
494 template <class Scalar>
495
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
15 class CheckPointTimeLoop : public TimeLoop<Scalar>
496 {
497 29140 class CheckPointType {
498 static constexpr std::size_t manualIdx = 0;
499 static constexpr std::size_t periodicIdx = 1;
500
501 public:
502 29140 bool isPeriodic() const { return set_[periodicIdx]; }
503 29032 bool isManual() const { return set_[manualIdx]; }
504 29032 bool isAny() const { return set_.any(); }
505
506 14570 CheckPointType& withPeriodic(bool value) { set_[periodicIdx] = value; return *this; }
507 14570 CheckPointType& withManual(bool value) { set_[manualIdx] = value; return *this; }
508
509 private:
510 std::bitset<2> set_;
511 };
512
513 public:
514 template<class... Args>
515 150 CheckPointTimeLoop(Args&&... args)
516 150 : TimeLoop<Scalar>(std::forward<Args>(args)...)
517 {
518 150 periodicCheckPoints_ = false;
519 150 deltaPeriodicCheckPoint_ = 0.0;
520 150 lastPeriodicCheckPoint_ = this->startTime_;
521 150 isCheckPoint_ = false;
522 150 }
523
524 /*!
525 * \brief Advance time step.
526 */
527 14516 void advanceTimeStep() override
528 {
529 29032 const auto dt = this->timeStepSize();
530 29032 const auto newTime = this->time()+dt;
531
532 //! Check point management, TimeLoop::isCheckPoint() has to be called after this!
533 14516 const auto cpType = nextCheckPointType_(newTime);
534
4/4
✓ Branch 0 taken 79 times.
✓ Branch 1 taken 14437 times.
✓ Branch 2 taken 79 times.
✓ Branch 3 taken 14437 times.
29032 if (cpType.isManual()) checkPoints_.pop();
535
4/4
✓ Branch 0 taken 701 times.
✓ Branch 1 taken 13815 times.
✓ Branch 2 taken 701 times.
✓ Branch 3 taken 13815 times.
29032 if (cpType.isPeriodic()) lastPeriodicCheckPoint_ += deltaPeriodicCheckPoint_;
536 14516 isCheckPoint_ = cpType.isAny();
537
538 14516 const auto previousTimeStepSize = this->previousTimeStepSize();
539
540 // advance the time step like in the parent class
541 14516 TimeLoop<Scalar>::advanceTimeStep();
542
543 // if this is a check point we might have reduced the time step to reach this check point
544 // reset the time step size to the time step size before this time step
545 14516 if (!this->willBeFinished())
546 {
547 using std::max;
548
2/2
✓ Branch 0 taken 712 times.
✓ Branch 1 taken 13574 times.
14286 if (isCheckPoint_)
549
2/2
✓ Branch 0 taken 99 times.
✓ Branch 1 taken 613 times.
811 this->setTimeStepSize(max(dt, previousTimeStepSize));
550
551 // if there is a check point soon check if the time step after the next time step would be smaller
552 // than 20% of the next time step, if yes increase the suggested next time step to exactly reach the check point
553 // (in the limits of the maximum time step size)
554 28572 auto nextDt = this->timeStepSize();
555 14286 const auto threshold = 0.2*nextDt;
556 28572 const auto nextTime = this->time() + nextDt;
557
558 14286 const auto nextDtToCheckPoint = maxDtToCheckPoint_(nextTime);
559
4/4
✓ Branch 0 taken 13577 times.
✓ Branch 1 taken 709 times.
✓ Branch 2 taken 13560 times.
✓ Branch 3 taken 17 times.
14286 if (nextDtToCheckPoint > Scalar{0} && Dune::FloatCmp::le(nextDtToCheckPoint, threshold))
560 17 nextDt += nextDtToCheckPoint;
561
562
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 14286 times.
14286 assert(nextDt > 0.0);
563 14286 this->setTimeStepSize(nextDt);
564 }
565 14516 }
566
567 /*!
568 * \brief The current maximum time step size
569 * \note This gets aligned on every setTimeStepSize call to end time
570 * and other possible check points
571 */
572 44275 Scalar maxTimeStepSize() const override
573 {
574 using std::min;
575 88550 const auto maxCheckPointDt = timeStepSizeToNextCheckPoint();
576 44275 const auto maxDtParent = TimeLoop<Scalar>::maxTimeStepSize();
577
2/2
✓ Branch 0 taken 10552 times.
✓ Branch 1 taken 33723 times.
44275 return min(maxDtParent, maxCheckPointDt);
578 }
579
580 /*!
581 * \brief Set a periodic check point
582 * \note You can query if we are at a time check point with isCheckPoint()
583 * \param interval Set a periodic checkout every [interval] seconds
584 * \param offset time from which the periodic check points are supposed to start (simulation time)
585 * the first checkpoint will be at time = offset.
586 * \note If offset is in the past the first check point will be at the next
587 * periodic check point greater or equal than time
588 * \note This also updates the time step size and potentially reduces the time step size to meet the next check point
589 */
590 template<class ScalarOrDuration1, class ScalarOrDuration2 = Scalar>
591 void setPeriodicCheckPoint(ScalarOrDuration1 interval, ScalarOrDuration2 offset = 0.0)
592 {
593
6/12
✓ Branch 1 taken 38 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 4 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 4 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 4 times.
✗ Branch 17 not taken.
55 setPeriodicCheckPoint_(
594 Detail::TimeLoop::toSeconds<Scalar>(interval),
595 Detail::TimeLoop::toSeconds<Scalar>(offset)
596 );
597 }
598
599 //! disable periodic check points
600 void disablePeriodicCheckPoints()
601 { periodicCheckPoints_ = false; }
602
603 //! remove all check points
604 void removeAllCheckPoints()
605 {
606 4 periodicCheckPoints_ = false;
607
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
8 while (!checkPoints_.empty())
608 checkPoints_.pop();
609 }
610
611 /*!
612 * \brief Whether now is a time checkpoint
613 * \note has to be called after TimeLoop::advanceTimeStep()
614 */
615 bool isCheckPoint() const
616 { return isCheckPoint_; }
617
618 /*!
619 * \brief add a checkpoint to the queue
620 * \note checkpoints have to be provided in ascending order
621 * \param t the check point (in seconds)
622 * \note This also updates the time step size and potentially reduces the time step size to meet the next check point
623 */
624 template<class ScalarOrDuration>
625 void setCheckPoint(ScalarOrDuration t)
626 {
627 // set the check point
628
4/8
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 4 times.
✗ Branch 11 not taken.
30 setCheckPoint_(Detail::TimeLoop::toSeconds<Scalar>(t));
629
630 // make sure we respect this check point on the next time step
631
4/8
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 4 times.
✗ Branch 11 not taken.
30 this->setTimeStepSize(this->timeStepSize());
632 }
633
634 /*!
635 * \brief add checkpoints to the queue from a list of time points
636 * \note checkpoints have to be provided in ascending order
637 * \param checkPoints the list of check points
638 * \note This also updates the time step size and potentially reduces the time step size to meet the next check point
639 */
640 template<class ScalarOrDuration>
641 void setCheckPoint(const std::initializer_list<ScalarOrDuration>& checkPoints)
642
2/4
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
10 { setCheckPoint(checkPoints.begin(), checkPoints.end()); }
643
644 template<class ScalarOrDuration>
645 [[deprecated("Use setCheckpoint(begin, end) instead. Will be removed after release 3.9.")]]
646 void setCheckPoint(const std::vector<ScalarOrDuration>& checkPoints)
647 { setCheckPoint(checkPoints.begin(), checkPoints.end()); }
648
649 /*!
650 * \brief add checkpoints to the queue from a container from the first iterator to the last iterator
651 * \note checkpoints have to be provided in ascending order
652 * \param first iterator to the first element to be inserted
653 * \param last iterator to the one-after-last element to be inserted
654 * \note This also updates the time step size and potentially reduces the time step size to meet the next check point
655 */
656 template<class ForwardIterator>
657 19 void setCheckPoint(ForwardIterator first, ForwardIterator last)
658 {
659 // set the check points
660
4/4
✓ Branch 0 taken 277 times.
✓ Branch 1 taken 19 times.
✓ Branch 2 taken 267 times.
✓ Branch 3 taken 14 times.
310 for (; first != last; ++first)
661 277 setCheckPoint_(Detail::TimeLoop::toSeconds<Scalar>(*first));
662
663 // make sure we respect this check point on the next time step
664 19 this->setTimeStepSize(this->timeStepSize());
665 19 }
666
667 /*!
668 * \brief Return the time step size to exactly reach the next check point
669 * In case there is no check point in the future, the largest representable scalar is returned
670 */
671 Scalar timeStepSizeToNextCheckPoint() const
672 44275 { return maxDtToCheckPoint_(this->time()); }
673
674 private:
675 bool fuzzyEqual_(const Scalar t0, const Scalar t1) const
676 36336 { return Dune::FloatCmp::eq(t0, t1, this->baseEps_*this->timeStepSize()); }
677
678 54 void setPeriodicCheckPoint_(Scalar interval, Scalar offset = 0.0)
679 {
680 using std::signbit;
681
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 54 times.
108 if (signbit(interval))
682 DUNE_THROW(Dune::InvalidStateException, "Interval has to be positive!");
683
684 54 periodicCheckPoints_ = true;
685 54 deltaPeriodicCheckPoint_ = interval;
686 54 lastPeriodicCheckPoint_ = offset;
687
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 54 times.
62 while (lastPeriodicCheckPoint_ + interval - this->time() < 1e-14*interval)
688 8 lastPeriodicCheckPoint_ += interval;
689
690
1/2
✓ Branch 0 taken 54 times.
✗ Branch 1 not taken.
54 if (this->verbose())
691 std::cout << Fmt::format("Enabled periodic check points every {:.5g} seconds ", interval)
692
5/14
✓ Branch 3 taken 54 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 54 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 54 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 54 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 54 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
162 << Fmt::format("with the next check point at {:.5g} seconds.\n", lastPeriodicCheckPoint_ + interval);
693
694 // check if the current time point is a periodic check point
695
4/4
✓ Branch 1 taken 48 times.
✓ Branch 2 taken 6 times.
✓ Branch 3 taken 48 times.
✓ Branch 4 taken 6 times.
54 if (nextCheckPointType_(this->time() + deltaPeriodicCheckPoint_).isPeriodic())
696 48 isCheckPoint_ = true;
697
698 // make sure we respect this check point on the next time step
699 54 this->setTimeStepSize(this->timeStepSize());
700 54 }
701
702 //! Adds a check point to the queue
703 307 void setCheckPoint_(Scalar t)
704 {
705
2/2
✓ Branch 0 taken 303 times.
✓ Branch 1 taken 4 times.
307 if (Dune::FloatCmp::le(t - this->time(), 0.0, this->timeStepSize()*this->baseEps_))
706 {
707
1/2
✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
8 if (this->verbose())
708 std::cerr << Fmt::format("Couldn't insert checkpoint at t = {:.5g} ", t)
709
5/14
✓ Branch 3 taken 8 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 8 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 8 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 8 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 8 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
24 << Fmt::format("because that's not in the future! (current simulation time is {:.5g})\n", this->time());
710 8 return;
711 }
712
713
4/4
✓ Branch 0 taken 264 times.
✓ Branch 1 taken 35 times.
✓ Branch 2 taken 264 times.
✓ Branch 3 taken 35 times.
598 if (!checkPoints_.empty())
714 {
715
4/4
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 261 times.
✓ Branch 2 taken 6 times.
✓ Branch 3 taken 258 times.
267 if (t <= checkPoints_.back())
716 {
717
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2 times.
6 if (this->verbose())
718 std::cerr << Fmt::format("Couldn't insert checkpoint at t = {:.5g} ", t)
719 8 << Fmt::format("because it's earlier than or equal to the last check point (t = {:.5g}) in the queue.\n", checkPoints_.back())
720
7/18
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 4 times.
✓ Branch 9 taken 4 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 4 times.
✗ Branch 13 not taken.
✓ Branch 16 taken 4 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 4 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 4 times.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
16 << "Checkpoints can only be inserted in ascending order." << std::endl;
721 6 return;
722 }
723 }
724
725
2/2
✓ Branch 0 taken 290 times.
✓ Branch 1 taken 3 times.
293 checkPoints_.push(t);
726
2/2
✓ Branch 0 taken 278 times.
✓ Branch 1 taken 15 times.
293 if (this->verbose())
727
2/6
✓ Branch 3 taken 278 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 278 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
556 std::cout << Fmt::format("Set check point at t = {:.5g} seconds.\n", t);
728 }
729
730 /*!
731 * \brief Return the type of (next) check point at the given time
732 */
733 14570 CheckPointType nextCheckPointType_(Scalar t)
734 {
735 return CheckPointType{}
736
4/4
✓ Branch 0 taken 6575 times.
✓ Branch 1 taken 731 times.
✓ Branch 2 taken 749 times.
✓ Branch 3 taken 6557 times.
13881 .withPeriodic(periodicCheckPoints_ && fuzzyEqual_(t - lastPeriodicCheckPoint_, deltaPeriodicCheckPoint_))
737
14/14
✓ Branch 0 taken 7306 times.
✓ Branch 1 taken 7264 times.
✓ Branch 2 taken 7306 times.
✓ Branch 3 taken 7264 times.
✓ Branch 4 taken 4806 times.
✓ Branch 5 taken 9764 times.
✓ Branch 6 taken 4806 times.
✓ Branch 7 taken 9764 times.
✓ Branch 8 taken 4727 times.
✓ Branch 9 taken 79 times.
✓ Branch 10 taken 4727 times.
✓ Branch 11 taken 79 times.
✓ Branch 12 taken 4727 times.
✓ Branch 13 taken 79 times.
43752 .withManual(!checkPoints_.empty() && fuzzyEqual_(t - checkPoints_.front(), 0.0));
738 }
739
740 /*!
741 * \brief Compute a time step size respecting upcoming checkpoints, starting from the given time t.
742 */
743 58561 Scalar maxDtToCheckPoint_(Scalar t) const
744 {
745 static constexpr auto unset = std::numeric_limits<Scalar>::max();
746 117122 const auto dtToPeriodic = dtToNextPeriodicCheckPoint_(t);
747 117122 const auto dtToManual = dtToNextManualCheckPoint_(t);
748
749 using std::min;
750
6/6
✓ Branch 0 taken 19292 times.
✓ Branch 1 taken 39269 times.
✓ Branch 2 taken 29576 times.
✓ Branch 3 taken 28985 times.
✓ Branch 4 taken 19288 times.
✓ Branch 5 taken 39273 times.
107429 return min(
751 dtToPeriodic.value_or(unset),
752 dtToManual.value_or(unset)
753 117122 );
754 }
755
756 /*!
757 * \brief Compute the time step size to the next periodic check point.
758 */
759 std::optional<Scalar> dtToNextPeriodicCheckPoint_(Scalar t) const
760 {
761
2/2
✓ Branch 0 taken 29576 times.
✓ Branch 1 taken 28985 times.
58561 if (periodicCheckPoints_)
762 59152 return lastPeriodicCheckPoint_ + deltaPeriodicCheckPoint_ - t;
763 57970 return {};
764 }
765
766 /*!
767 * \brief Compute a time step size respecting the next manual check point.
768 */
769 std::optional<Scalar> dtToNextManualCheckPoint_(Scalar t) const
770 {
771
4/4
✓ Branch 0 taken 19292 times.
✓ Branch 1 taken 39269 times.
✓ Branch 2 taken 19292 times.
✓ Branch 3 taken 39269 times.
117122 if (!checkPoints_.empty())
772 57876 return checkPoints_.front() - t;
773 78538 return {};
774 }
775
776 bool periodicCheckPoints_;
777 Scalar deltaPeriodicCheckPoint_;
778 Scalar lastPeriodicCheckPoint_;
779 std::queue<Scalar> checkPoints_;
780 bool isCheckPoint_;
781 };
782
783 template<class... Args>
784 CheckPointTimeLoop(Args&&... args) -> CheckPointTimeLoop<
785 typename decltype(TimeLoop{std::forward<Args>(args)...})::Scalar
786 >;
787
788 } // end namespace Dumux
789
790 #endif
791