timeloop.hh 21 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

42
43
namespace Dumux {

44
/*!
45
 * \ingroup Common
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
 * \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.
63
64
 *
 * \note Time and time step sizes are in units of seconds
65
 */
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
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;

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

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

    /*!
     * \brief Advance to the next time step.
92
93
94
95
96
97
98
99
     */
    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;
100
101
102
103
104

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

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

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

    /*!
     * \brief Tells the time loop to start tracking the time.
     */
    void start()
    {
        timer_.start();
132
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
    }

    /*!
     * \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;

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

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

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

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

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

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

    /*!
     * \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().
     */
218
    Scalar time() const final
219
220
221
222
223
224
225
226
227
228
229
230
231
232
    { 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)
233
234
235
236
237
    {
        endTime_ = t;
        if (verbose_)
            std::cout << "Set new end time to t = " << t << " seconds." << std::endl;
    }
238
239

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

    /*!
     * \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$
     */
255
    void setTimeStepSize(Scalar dt) final
256
257
    {
        using std::min;
258
        timeStepSize_ = min(dt, maxTimeStepSize());
259
260
261
262
263
264
        if (!finished() && Dune::FloatCmp::le(timeStepSize_, 0.0, 1e-14*endTime_))
        {
            std::cerr << "You have set a very small timestep size (dt = "
                      << timeStepSize_ << "). This might lead to numerical problems!"
                      << std::endl;
        }
265
266
267
268
269
    }

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

    /*!
     * \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.
     */
284
    Scalar timeStepSize() const final
285
286
287
288
289
290
291
292
293
    { return timeStepSize_; }

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

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

300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
    /*!
     * \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.
     */
316
    bool finished() const override
317
    {
318
        return finished_ || endTime_-time_ < 1e-10*time_;
319
    }
320
321
322
323
324
325

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

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

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

344
345
    /*!
     * \brief State info on cpu time.
346
     * \note Always call this after TimeLoop::advanceTimeStep()
347
     */
348
    void reportTimeStep() const
349
350
351
    {
        if (verbose_)
        {
352
353
354
355
356
357
358
359
            const auto cpuTime = wallClockTime();
            const auto percent = std::round( time_ / endTime_ * 100 );
            std::cout << "[" << std::fixed << std::setw( 3 ) << std::setfill( ' ' )
                      << std::setprecision( 0 )  << percent << "%] "
                      << "Time step " << timeStepIdx_ << " done in "
                      << std::setprecision( 6 ) << timeStepWallClockTime_ << " seconds. "
                      << "Wall clock time: " << std::setprecision( 3 ) << cpuTime
                      << ", time: " << std::setprecision( 5 ) << time_
360
                      << ", time step size: " << std::setprecision( 8 ) << previousTimeStepSize_
361
362
363
364
365
366
367
                      << std::endl;
        }
    }

    /*!
     * \brief Print final status and stops tracking the time.
     */
368
    template< class Communicator = Dune::CollectiveCommunication<typename Dune::MPIHelper::MPICommunicator> >
369
370
    void finalize(const Communicator& comm = Dune::MPIHelper::getCollectiveCommunication())
    {
371
        auto cpuTime = timer_.stop();
372
373
374

        if (verbose_)
        {
375
            std::cout << "Simulation took " << cpuTime << " seconds on "
376
377
378
379
                      << comm.size() << " processes.\n";
        }

        if (comm.size() > 1)
380
            cpuTime = comm.sum(cpuTime);
381
382
383

        if (verbose_)
        {
384
            std::cout << "The cumulative CPU time was " << cpuTime << " seconds.\n";
385
386
387
        }
    }

388
389
390
391
    //! If the time loop has verbose output
    bool verbose() const
    { return verbose_; }

392
393
394
395
    //! Sets time loop verbosity
    void setVerbose(bool verbose = true)
    { verbose_ = verbose; }

396
397
398
399
    /*
     * @}
     */

400
protected:
401
402
403
404
405
    Dune::Timer timer_;
    Scalar time_;
    Scalar endTime_;

    Scalar timeStepSize_;
406
    Scalar previousTimeStepSize_;
407
    Scalar userSetMaxTimeStepSize_;
408
    Scalar timeAfterLastTimeStep_, timeStepWallClockTime_;
409
410
411
412
    int timeStepIdx_;
    bool finished_;
    bool verbose_;
};
413

414
415
416
417
/*!
 * \ingroup Common
 * \brief A time loop with a check point mechanism
 */
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
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.
     */
434
    void advanceTimeStep() override
435
    {
436
        const auto dt = this->timeStepSize();
437
        const auto newTime = this->time()+dt;
438
439
440

        //! Check point management, TimeLoop::isCheckPoint() has to be called after this!
        // if we reached a periodic check point
441
        if (periodicCheckPoints_ && Dune::FloatCmp::eq(newTime - lastPeriodicCheckPoint_, deltaPeriodicCheckPoint_, 1e-7))
442
443
444
445
446
447
        {
            lastPeriodicCheckPoint_ += deltaPeriodicCheckPoint_;
            isCheckPoint_ = true;
        }

        // or a manually set check point
448
        else if (!checkPoints_.empty() && Dune::FloatCmp::eq(newTime - checkPoints_.front(), 0.0, 1e-7))
449
450
451
452
453
454
455
456
457
458
        {
            checkPoints_.pop();
            isCheckPoint_ = true;
        }

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

460
461
        const auto previousTimeStepSize = this->previousTimeStepSize();

462
463
        // advance the time step like in the parent class
        TimeLoop<Scalar>::advanceTimeStep();
464
465
466

        // 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
467
468
469
470
471
        if (!this->willBeFinished())
        {
            using std::max;
            if (isCheckPoint_)
                this->setTimeStepSize(max(dt, previousTimeStepSize));
472

473
474
475
476
477
478
            // 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;
479

480
481
            if (periodicCheckPoints_ && Dune::FloatCmp::le(lastPeriodicCheckPoint_ + deltaPeriodicCheckPoint_ - nextTime, threshold))
                nextDt = lastPeriodicCheckPoint_ + deltaPeriodicCheckPoint_ - this->time();
482

483
484
            if (!checkPoints_.empty() && Dune::FloatCmp::le(checkPoints_.front() - nextTime, threshold))
                nextDt = checkPoints_.front() - this->time();
485

486
487
488
            assert(nextDt > 0.0);
            this->setTimeStepSize(nextDt);
        }
489
490
491
    }

    /*!
492
493
494
     * \brief The current maximum time step size
     * \note This gets aligned on every setTimeStepSize call to end time
     *       and other possible check points
495
     */
496
    Scalar maxTimeStepSize() const override
497
498
    {
        using std::min;
499
500
501
        const auto maxCheckPointDt = computeStepSizeRespectingCheckPoints_();
        const auto maxDtParent = TimeLoop<Scalar>::maxTimeStepSize();
        return min(maxDtParent, maxCheckPointDt);
502
503
504
    }

    /*!
505
506
     * \brief Set a periodic check point
     * \note You can query if we are at a time check point with isCheckPoint()
507
     * \param interval Set a periodic checkout every [interval] seconds
508
509
     * \param offset time from which the periodic check points are supposed to start (simulation time)
     *        the first checkpoint will be at time = offset.
510
511
     * \note If offset is in the past the first check point will be at the next
     *       periodic check point greater or equal than time
512
     * \note This also updates the time step size and potentially reduces the time step size to meet the next check point
513
     */
514
    void setPeriodicCheckPoint(Scalar interval, Scalar offset = 0.0)
515
    {
516
517
518
519
        using std::signbit;
        if (signbit(interval))
            DUNE_THROW(Dune::InvalidStateException, "Interval has to be positive!");

520
521
        periodicCheckPoints_ = true;
        deltaPeriodicCheckPoint_ = interval;
522
523
524
525
        lastPeriodicCheckPoint_ = offset;
        while (lastPeriodicCheckPoint_ + interval - this->time() < 1e-14*interval)
            lastPeriodicCheckPoint_ += interval;

526
        if (this->verbose())
527
528
            std::cout << "Enabled periodic check points every " << interval
                      << " seconds with the next check point at " << lastPeriodicCheckPoint_ + interval << " seconds." << std::endl;
529

530
        // check if the current time point is a check point
531
        if (Dune::FloatCmp::eq(this->time()-lastPeriodicCheckPoint_, 0.0, 1e-7))
532
            isCheckPoint_ = true;
533
534
535

        // make sure we respect this check point on the next time step
        this->setTimeStepSize(this->timeStepSize());
536
537
538
539
540
541
542
543
544
545
546
547
    }

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

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

550
551
552
553
    /*!
     * \brief Whether now is a time checkpoint
     * \note has to be called after TimeLoop::advanceTimeStep()
     */
554
555
556
    bool isCheckPoint() const
    { return isCheckPoint_; }

557
558
559
560
    /*!
     * \brief add a checkpoint to the queue
     * \note checkpoints have to be provided in ascending order
     * \param t the check point (in seconds)
561
     * \note This also updates the time step size and potentially reduces the time step size to meet the next check point
562
     */
563
    void setCheckPoint(Scalar t)
564
565
566
567
568
    {
        // set the check point
        setCheckPoint_(t);

        // make sure we respect this check point on the next time step
569
        this->setTimeStepSize(this->timeStepSize());
570
571
572
573
574
575
    }

    /*!
     * \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
576
     * \note This also updates the time step size and potentially reduces the time step size to meet the next check point
577
578
579
580
581
582
583
584
585
     */
    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
586
     * \note This also updates the time step size and potentially reduces the time step size to meet the next check point
587
588
589
590
591
592
593
594
595
     */
    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
596
        this->setTimeStepSize(this->timeStepSize());
597
598
599
600
601
    }

private:
    //! Adds a check point to the queue
    void setCheckPoint_(Scalar t)
602
    {
603
        if (Dune::FloatCmp::le(t - this->time(), 0.0, this->timeStepSize()*1e-7))
604
605
606
607
608
609
610
        {
            if (this->verbose())
                std::cerr << "Couldn't insert checkpoint at t = " << t
                          << " because that's in the past! (current simulation time is " << this->time() << ")" << std::endl;
            return;
        }

611
612
613
614
615
        if (!checkPoints_.empty())
        {
            if (t < checkPoints_.back())
            {
                if (this->verbose())
616
617
618
                    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;
619
620
            }
        }
621
622
623
624

        checkPoints_.push(t);
        if (this->verbose())
            std::cout << "Set check point at t = " << t << " seconds." << std::endl;
625
626
    }

627
    /*!
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
     * \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
652
653

#endif