timeloop.hh 20.8 KB
Newer Older
1
2
3
4
5
6
7
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   This program is free software: you can redistribute it and/or modify    *
 *   it under the terms of the GNU General Public License as published by    *
8
 *   the Free Software Foundation, either version 3 of the License, or       *
9
10
11
12
13
14
15
16
17
18
19
20
 *   (at your option) any later version.                                     *
 *                                                                           *
 *   This program is distributed in the hope that it will be useful,         *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
 *   GNU General Public License for more details.                            *
 *                                                                           *
 *   You should have received a copy of the GNU General Public License       *
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
/*!
 * \file
21
 * \ingroup Common
22
23
24
25
26
27
 * \brief Manages the handling of time dependent problems
 */
#ifndef DUMUX_TIME_LOOP_HH
#define DUMUX_TIME_LOOP_HH

#include <algorithm>
28
#include <queue>
29
#include <iomanip>
30
31
32

#include <dune/common/float_cmp.hh>
#include <dune/common/timer.hh>
33
34
35
36

#include <dune/common/version.hh>
#include <dune/common/parallel/communication.hh>

37
#include <dune/common/parallel/mpihelper.hh>
38
#include <dune/common/exceptions.hh>
39

Timo Koch's avatar
Timo Koch committed
40
#include <dumux/common/parameters.hh>
41
#include <dumux/io/format.hh>
42

43
44
namespace Dumux {

45
/*!
46
 * \ingroup Common
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
 * \brief Manages the handling of time dependent problems.
 *
 * This class facilitates the time management of the simulation.
 * It doesn't manage any user data, but keeps track of what the
 * current time, time step size and "episode" of the
 * simulation is. It triggers the initialization of the problem and
 * is responsible for the time control of a simulation run.
 *
 * The time manager allows to specify a sequence of "episodes" which
 * determine the boundary conditions of a problem. This approach
 * is handy if the problem is not static, i.e. the boundary
 * conditions change over time.
 *
 * An episode is a span of simulated time in which
 * the problem behaves in a specific way. It is characterized by
 * the (simulation) time it starts, its length and a consecutive
 * index starting at 0.
64
65
 *
 * \note Time and time step sizes are in units of seconds
66
 */
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
template<class Scalar>
class TimeLoopBase
{
public:
    //! Abstract base class needs virtual constructor
    virtual ~TimeLoopBase() {};

    /*!
     * \brief Return the time \f$\mathrm{[s]}\f$ before the time integration.
     * To get the time after the time integration you have to add timeStepSize() to
     * time().
     */
    virtual Scalar time() const = 0;

    /*!
82
     * \brief Returns the suggested time step length \f$\mathrm{[s]}\f$
83
84
     */
    virtual Scalar timeStepSize() const = 0;
85
86

    /*!
87
88
89
90
91
92
     * \brief Get the maximum possible time step size \f$\mathrm{[s]}\f$
     */
    virtual Scalar maxTimeStepSize() const = 0;

    /*!
     * \brief Advance to the next time step.
93
94
95
96
97
98
99
100
     */
    virtual void advanceTimeStep() = 0;

    /*!
     * \brief Set the current time step size to a given value.
     * \param dt The new value for the time step size \f$\mathrm{[s]}\f$
     */
    virtual void setTimeStepSize(Scalar dt) = 0;
101
102
103
104
105

    /*!
     * \brief Returns true if the simulation is finished.
     */
    virtual bool finished() const = 0;
106
107
};

108
109
110
111
/*!
 * \ingroup Common
 * \brief The default time loop for instationary simulations
 */
112
template <class Scalar>
113
class TimeLoop : public TimeLoopBase<Scalar>
114
115
{
public:
116
117
    TimeLoop(Scalar startTime, Scalar dt, Scalar tEnd, bool verbose = true)
    : timer_(false)
118
    {
119
        reset(startTime, dt, tEnd, verbose);
120
121
122
123
124
125
126
127
128
129
130
131
132
    }

    /*!
     *  \name Simulated time and time step management
     * @{
     */

    /*!
     * \brief Tells the time loop to start tracking the time.
     */
    void start()
    {
        timer_.start();
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
    }

    /*!
     * \brief Tells the time loop to stop tracking the time.
     * \return the wall clock time (CPU time) spent until now
     */
    double stop()
    {
        return timer_.stop();
    }

    /*!
     * \brief Reset the timer
     */
    void resetTimer()
    {
        timer_.reset();
    }

    /*!
     * \brief Reset the time loop
     */
    void reset(Scalar startTime, Scalar dt, Scalar tEnd, bool verbose = true)
    {
        verbose_ =
            verbose &&
            Dune::MPIHelper::getCollectiveCommunication().rank() == 0;

        time_ = startTime;
        endTime_ = tEnd;

164
        previousTimeStepSize_ = 0.0;
165
        userSetMaxTimeStepSize_ = std::numeric_limits<Scalar>::max();
166
167
        timeStepIdx_ = 0;
        finished_ = false;
168
169
        timeAfterLastTimeStep_ = 0.0;
        timeStepWallClockTime_ = 0.0;
170

171
172
173
        // ensure that dt is not greater than tEnd-startTime
        setTimeStepSize(dt);

174
175
        timer_.stop();
        timer_.reset();
176
177
178
179
180
    }

    /*!
     * \brief Advance time step.
     */
181
    void advanceTimeStep() override
182
183
184
    {
        timeStepIdx_++;
        time_ += timeStepSize_;
185
        previousTimeStepSize_ = timeStepSize_;
186
187
188
189
190

        // compute how long the last time step took
        const auto cpuTime = wallClockTime();
        timeStepWallClockTime_ = cpuTime - timeAfterLastTimeStep_;
        timeAfterLastTimeStep_ = cpuTime;
191
192

        // ensure that using current dt we don't exceed tEnd in next time step
193
        setTimeStepSize(timeStepSize_);
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
    }

    /*!
     * \brief Set the current simulated time, don't change the current
     *        time step index.
     *
     * \param t The time \f$\mathrm{[s]}\f$ which should be jumped to
     */
    void setTime(Scalar t)
    { time_ = t; }

    /*!
     * \brief Set the current simulated time and the time step index.
     *
     * \param t The time \f$\mathrm{[s]}\f$ which should be jumped to
     * \param stepIdx The new time step index
     */
    void setTime(Scalar t, int stepIdx)
    { time_ = t; timeStepIdx_ = stepIdx; }

    /*!
     * \brief Return the time \f$\mathrm{[s]}\f$ before the time integration.
     * To get the time after the time integration you have to add timeStepSize() to
     * time().
     */
219
    Scalar time() const final
220
221
222
223
224
225
226
227
228
229
230
231
232
233
    { return time_; }

    /*!
     * \brief Returns the number of (simulated) seconds which the simulation runs.
     */
    Scalar endTime() const
    { return endTime_; }

    /*!
     * \brief Set the time of simulated seconds at which the simulation runs.
     *
     * \param t The time \f$\mathrm{[s]}\f$ at which the simulation is finished
     */
    void setEndTime(Scalar t)
234
235
236
    {
        endTime_ = t;
        if (verbose_)
237
            std::cout << Fmt::format("Set new end time to t = {:.5g} seconds.\n", t);
238
    }
239
240

    /*!
241
     * \brief Returns the current wall clock time (cpu time) spend in this time loop
242
     */
243
    double wallClockTime() const
244
    { return timer_.elapsed(); }
245
246
247
248
249
250
251
252
253
254
255

    /*!
     * \brief Set the current time step size to a given value.
     *
     * If the step size would exceed the length of the current
     * episode, the timeStep() method will take care that the step
     * size won't exceed the episode or the end of the simulation,
     * though.
     *
     * \param dt The new value for the time step size \f$\mathrm{[s]}\f$
     */
256
    void setTimeStepSize(Scalar dt) final
257
258
    {
        using std::min;
259
        timeStepSize_ = min(dt, maxTimeStepSize());
260
        if (!finished() && Dune::FloatCmp::le(timeStepSize_, 0.0, 1e-14*endTime_))
261
            std::cerr << Fmt::format("You have set a very small timestep size (dt = {:.5g}).", timeStepSize_)
262
                      << " This might lead to numerical problems!\n";
263
264
265
266
267
    }

    /*!
     * \brief Set the maximum time step size to a given value.
     *
268
     * \param maxDt The new value for the maximum time step size \f$\mathrm{[s]}\f$
269
     * \note This also updates the time step size
270
271
     */
    void setMaxTimeStepSize(Scalar maxDt)
272
    {
273
        userSetMaxTimeStepSize_ = maxDt;
274
        setTimeStepSize(timeStepSize_);
275
    }
276
277
278
279
280
281

    /*!
     * \brief Returns the suggested time step length \f$\mathrm{[s]}\f$ so that we
     *        don't miss the beginning of the next episode or cross
     *        the end of the simulation.
     */
282
    Scalar timeStepSize() const final
283
284
285
286
287
288
289
290
291
    { return timeStepSize_; }

    /*!
     * \brief Returns number of time steps which have been
     *        executed since the beginning of the simulation.
     */
    int timeStepIndex() const
    { return timeStepIdx_; }

292
293
294
295
296
297
    /*!
     * \brief The previous time step size
     */
    Scalar previousTimeStepSize() const
    { return previousTimeStepSize_; }

298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
    /*!
     * \brief Specify whether the simulation is finished
     *
     * \param finished If true the simulation is considered finished
     *                 before the end time is reached, else it is only
     *                 considered finished if the end time is reached.
     */
    void setFinished(bool finished = true)
    { finished_ = finished; }

    /*!
     * \brief Returns true if the simulation is finished.
     *
     * This is the case if either setFinished(true) has been called or
     * if the end time is reached.
     */
314
    bool finished() const override
315
    {
316
        return finished_ || endTime_-time_ < 1e-10*time_;
317
    }
318
319
320
321
322
323

    /*!
     * \brief Returns true if the simulation is finished after the
     *        time level is incremented by the current time step size.
     */
    bool willBeFinished() const
324
    {
325
        return finished() || endTime_-time_-timeStepSize_ < 1e-10*timeStepSize_;
326
    }
327
328

    /*!
329
330
331
     * \brief The current maximum time step size
     * \note This gets aligned on every setTimeStepSize call to end time
     *       and other possible check points
332
     */
333
    Scalar maxTimeStepSize() const override
334
    {
335
336
337
338
339
        if (finished())
            return 0.0;

        using std::min; using std::max;
        return min(userSetMaxTimeStepSize_, max<Scalar>(0.0, endTime_ - time_));
340
    }
341

342
343
    /*!
     * \brief State info on cpu time.
344
     * \note Always call this after TimeLoop::advanceTimeStep()
345
     */
346
    void reportTimeStep() const
347
348
349
    {
        if (verbose_)
        {
350
            const auto cpuTime = wallClockTime();
351
            using std::round;
352
            const auto percent = round( time_ / endTime_ * 100 );
353
            std::cout << Fmt::format("[{:3.0f}%] ", percent)
354
                      << Fmt::format("Time step {} done in {:.2g} seconds. ", timeStepIdx_, timeStepWallClockTime_)
355
                      << Fmt::format("Wall clock time: {:.5g}, time: {:.5g}, time step size: {:.5g}\n", cpuTime, time_, previousTimeStepSize_);
356
357
358
359
360
361
        }
    }

    /*!
     * \brief Print final status and stops tracking the time.
     */
362
    template< class Communicator = Dune::CollectiveCommunication<typename Dune::MPIHelper::MPICommunicator> >
363
364
    void finalize(const Communicator& comm = Dune::MPIHelper::getCollectiveCommunication())
    {
365
        auto cpuTime = timer_.stop();
366
367

        if (verbose_)
368
            std::cout << Fmt::format("Simulation took {:.5g} seconds on {} processes.\n", cpuTime, comm.size());
369
370

        if (comm.size() > 1)
371
            cpuTime = comm.sum(cpuTime);
372
373

        if (verbose_)
374
            std::cout << Fmt::format("The cumulative CPU time was {:.5g} seconds.\n", cpuTime);
375
376
    }

377
378
379
380
    //! If the time loop has verbose output
    bool verbose() const
    { return verbose_; }

381
382
383
384
    //! Sets time loop verbosity
    void setVerbose(bool verbose = true)
    { verbose_ = verbose; }

385
386
387
388
    /*
     * @}
     */

389
protected:
390
391
392
393
394
    Dune::Timer timer_;
    Scalar time_;
    Scalar endTime_;

    Scalar timeStepSize_;
395
    Scalar previousTimeStepSize_;
396
    Scalar userSetMaxTimeStepSize_;
397
    Scalar timeAfterLastTimeStep_, timeStepWallClockTime_;
398
399
400
401
    int timeStepIdx_;
    bool finished_;
    bool verbose_;
};
402

403
404
405
406
/*!
 * \ingroup Common
 * \brief A time loop with a check point mechanism
 */
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
template <class Scalar>
class CheckPointTimeLoop : public TimeLoop<Scalar>
{
public:
    CheckPointTimeLoop(Scalar startTime, Scalar dt, Scalar tEnd, bool verbose = true)
    : TimeLoop<Scalar>(startTime, dt, tEnd, verbose)
    {
        periodicCheckPoints_ = false;
        deltaPeriodicCheckPoint_ = 0.0;
        lastPeriodicCheckPoint_ = startTime;
        isCheckPoint_ = false;
    }

    /*!
     * \brief Advance time step.
     */
423
    void advanceTimeStep() override
424
    {
425
        const auto dt = this->timeStepSize();
426
        const auto newTime = this->time()+dt;
427
428
429

        //! Check point management, TimeLoop::isCheckPoint() has to be called after this!
        // if we reached a periodic check point
430
        if (periodicCheckPoints_ && Dune::FloatCmp::eq(newTime - lastPeriodicCheckPoint_, deltaPeriodicCheckPoint_, 1e-7))
431
432
433
434
435
436
        {
            lastPeriodicCheckPoint_ += deltaPeriodicCheckPoint_;
            isCheckPoint_ = true;
        }

        // or a manually set check point
437
        else if (!checkPoints_.empty() && Dune::FloatCmp::eq(newTime - checkPoints_.front(), 0.0, 1e-7))
438
439
440
441
442
443
444
445
446
447
        {
            checkPoints_.pop();
            isCheckPoint_ = true;
        }

        // if not reset the check point flag
        else
        {
            isCheckPoint_ = false;
        }
448

449
450
        const auto previousTimeStepSize = this->previousTimeStepSize();

451
452
        // advance the time step like in the parent class
        TimeLoop<Scalar>::advanceTimeStep();
453
454
455

        // if this is a check point we might have reduced the time step to reach this check point
        // reset the time step size to the time step size before this time step
456
457
458
459
460
        if (!this->willBeFinished())
        {
            using std::max;
            if (isCheckPoint_)
                this->setTimeStepSize(max(dt, previousTimeStepSize));
461

462
463
464
465
466
467
            // if there is a check point soon check if the time step after the next time step would be smaller
            // than 20% of the next time step, if yes increase the suggested next time step to exactly reach the check point
            // (in the limits of the maximum time step size)
            auto nextDt = this->timeStepSize();
            const auto threshold = 0.2*nextDt;
            const auto nextTime = this->time() + nextDt;
468

469
470
            if (periodicCheckPoints_ && Dune::FloatCmp::le(lastPeriodicCheckPoint_ + deltaPeriodicCheckPoint_ - nextTime, threshold))
                nextDt = lastPeriodicCheckPoint_ + deltaPeriodicCheckPoint_ - this->time();
471

472
473
            if (!checkPoints_.empty() && Dune::FloatCmp::le(checkPoints_.front() - nextTime, threshold))
                nextDt = checkPoints_.front() - this->time();
474

475
476
477
            assert(nextDt > 0.0);
            this->setTimeStepSize(nextDt);
        }
478
479
480
    }

    /*!
481
482
483
     * \brief The current maximum time step size
     * \note This gets aligned on every setTimeStepSize call to end time
     *       and other possible check points
484
     */
485
    Scalar maxTimeStepSize() const override
486
487
    {
        using std::min;
488
489
490
        const auto maxCheckPointDt = computeStepSizeRespectingCheckPoints_();
        const auto maxDtParent = TimeLoop<Scalar>::maxTimeStepSize();
        return min(maxDtParent, maxCheckPointDt);
491
492
493
    }

    /*!
494
495
     * \brief Set a periodic check point
     * \note You can query if we are at a time check point with isCheckPoint()
496
     * \param interval Set a periodic checkout every [interval] seconds
497
498
     * \param offset time from which the periodic check points are supposed to start (simulation time)
     *        the first checkpoint will be at time = offset.
499
500
     * \note If offset is in the past the first check point will be at the next
     *       periodic check point greater or equal than time
501
     * \note This also updates the time step size and potentially reduces the time step size to meet the next check point
502
     */
503
    void setPeriodicCheckPoint(Scalar interval, Scalar offset = 0.0)
504
    {
505
506
507
508
        using std::signbit;
        if (signbit(interval))
            DUNE_THROW(Dune::InvalidStateException, "Interval has to be positive!");

509
510
        periodicCheckPoints_ = true;
        deltaPeriodicCheckPoint_ = interval;
511
512
513
514
        lastPeriodicCheckPoint_ = offset;
        while (lastPeriodicCheckPoint_ + interval - this->time() < 1e-14*interval)
            lastPeriodicCheckPoint_ += interval;

515
        if (this->verbose())
516
517
            std::cout << Fmt::format("Enabled periodic check points every {:.5g} seconds ", interval)
                      << Fmt::format("with the next check point at {:.5g} seconds.\n", lastPeriodicCheckPoint_ + interval);
518

519
        // check if the current time point is a check point
520
        if (Dune::FloatCmp::eq(this->time()-lastPeriodicCheckPoint_, 0.0, 1e-7))
521
            isCheckPoint_ = true;
522
523
524

        // make sure we respect this check point on the next time step
        this->setTimeStepSize(this->timeStepSize());
525
526
527
528
529
530
531
532
533
534
535
536
    }

    //! disable periodic check points
    void disablePeriodicCheckPoints()
    { periodicCheckPoints_ = false; }

    //! remove all check points
    void removeAllCheckPoints()
    {
        periodicCheckPoints_ = false;
        while (!checkPoints_.empty())
            checkPoints_.pop();
537
538
    }

539
540
541
542
    /*!
     * \brief Whether now is a time checkpoint
     * \note has to be called after TimeLoop::advanceTimeStep()
     */
543
544
545
    bool isCheckPoint() const
    { return isCheckPoint_; }

546
547
548
549
    /*!
     * \brief add a checkpoint to the queue
     * \note checkpoints have to be provided in ascending order
     * \param t the check point (in seconds)
550
     * \note This also updates the time step size and potentially reduces the time step size to meet the next check point
551
     */
552
    void setCheckPoint(Scalar t)
553
554
555
556
557
    {
        // set the check point
        setCheckPoint_(t);

        // make sure we respect this check point on the next time step
558
        this->setTimeStepSize(this->timeStepSize());
559
560
561
562
563
564
    }

    /*!
     * \brief add checkpoints to the queue from a vector of time points
     * \note checkpoints have to be provided in ascending order
     * \param checkPoints the vector of check points
565
     * \note This also updates the time step size and potentially reduces the time step size to meet the next check point
566
567
568
569
570
571
572
573
574
     */
    void setCheckPoint(const std::vector<Scalar>& checkPoints)
    { setCheckPoint(checkPoints.begin(), checkPoints.end()); }

    /*!
     * \brief add checkpoints to the queue from a container from the first iterator to the last iterator
     * \note checkpoints have to be provided in ascending order
     * \param first iterator to the first element to be inserted
     * \param last iterator to the one-after-last element to be inserted
575
     * \note This also updates the time step size and potentially reduces the time step size to meet the next check point
576
577
578
579
580
581
582
583
584
     */
    template<class ForwardIterator>
    void setCheckPoint(ForwardIterator first, ForwardIterator last)
    {
        // set the check points
        for (; first != last; ++first)
            setCheckPoint_(*first);

        // make sure we respect this check point on the next time step
585
        this->setTimeStepSize(this->timeStepSize());
586
587
588
589
590
    }

private:
    //! Adds a check point to the queue
    void setCheckPoint_(Scalar t)
591
    {
592
        if (Dune::FloatCmp::le(t - this->time(), 0.0, this->timeStepSize()*1e-7))
593
594
        {
            if (this->verbose())
595
596
                std::cerr << Fmt::format("Couldn't insert checkpoint at t = {:.5g} ", t)
                          << Fmt::format("because that's in the past! (current simulation time is {:.5g})\n", this->time());
597
598
599
            return;
        }

600
601
602
603
604
        if (!checkPoints_.empty())
        {
            if (t < checkPoints_.back())
            {
                if (this->verbose())
605
606
607
                    std::cerr << "Couldn't insert checkpoint as it is earlier than the last check point in the queue.\n"
                              << "Checkpoints can only be inserted in ascending order." << std::endl;
                return;
608
609
            }
        }
610
611
612

        checkPoints_.push(t);
        if (this->verbose())
613
            std::cout << Fmt::format("Set check point at t = {:.5g} seconds.\n", t);
614
615
    }

616
    /*!
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
     * \brief Aligns dt to the next check point
     */
    Scalar computeStepSizeRespectingCheckPoints_() const
    {
        using std::min;
        auto maxDt = std::numeric_limits<Scalar>::max();

        if (periodicCheckPoints_)
            maxDt = min(maxDt, lastPeriodicCheckPoint_ + deltaPeriodicCheckPoint_ - this->time());

        if (!checkPoints_.empty())
            maxDt = min(maxDt, checkPoints_.front() - this->time());

        return maxDt;
    }

    bool periodicCheckPoints_;
    Scalar deltaPeriodicCheckPoint_;
    Scalar lastPeriodicCheckPoint_;
    std::queue<Scalar> checkPoints_;
    bool isCheckPoint_;
};

} // end namespace Dumux
641
642

#endif