common/timeloop.hh

Go to the documentation of this file.

// -*- 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    *
 *   the Free Software Foundation, either version 3 of the License, or       *
 *   (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/>.   *
 *****************************************************************************/
#ifndef DUMUX_TIME_LOOP_HH
#define DUMUX_TIME_LOOP_HH
 
#include <algorithm>
#include <queue>
#include <iomanip>
 
#include <dune/common/float_cmp.hh>
#include <dune/common/timer.hh>
 
#include <dune/common/parallel/communication.hh>
 
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/exceptions.hh>
 
#include <dumux/common/parameters.hh>
#include <dumux/io/format.hh>
 
namespace Dumux {
 
template<class Scalar>
class TimeLoopBase
{
public:
    virtual ~TimeLoopBase() {};
 
    virtual Scalar time() const = 0;
 
    virtual Scalar timeStepSize() const = 0;
 
    virtual Scalar maxTimeStepSize() const = 0;
 
    virtual void advanceTimeStep() = 0;
 
    virtual void setTimeStepSize(Scalar dt) = 0;
 
    virtual bool finished() const = 0;
};
 
template <class Scalar>
class TimeLoop : public TimeLoopBase<Scalar>
{
public:
    TimeLoop(Scalar startTime, Scalar dt, Scalar tEnd, bool verbose = true)
    : timer_(false)
    {
        reset(startTime, dt, tEnd, verbose);
    }
 
    void start()
    {
        timer_.start();
    }
 
    double stop()
    {
        return timer_.stop();
    }
 
    void resetTimer()
    {
        timer_.reset();
    }
 
    void reset(Scalar startTime, Scalar dt, Scalar tEnd, bool verbose = true)
    {
        verbose_ =
            verbose &&
            Dune::MPIHelper::getCommunication().rank() == 0;
 
        time_ = startTime;
        endTime_ = tEnd;
 
        previousTimeStepSize_ = 0.0;
        userSetMaxTimeStepSize_ = std::numeric_limits<Scalar>::max();
        timeStepIdx_ = 0;
        finished_ = false;
        timeAfterLastTimeStep_ = 0.0;
        timeStepWallClockTime_ = 0.0;
 
        // ensure that dt is not greater than tEnd-startTime
        setTimeStepSize(dt);
 
        timer_.stop();
        timer_.reset();
    }
 
    void advanceTimeStep() override
    {
        timeStepIdx_++;
        time_ += timeStepSize_;
        previousTimeStepSize_ = timeStepSize_;
 
        // compute how long the last time step took
        const auto cpuTime = wallClockTime();
        timeStepWallClockTime_ = cpuTime - timeAfterLastTimeStep_;
        timeAfterLastTimeStep_ = cpuTime;
 
        // ensure that using current dt we don't exceed tEnd in next time step
        setTimeStepSize(timeStepSize_);
    }
 
    void setTime(Scalar t)
    { time_ = t; }
 
    void setTime(Scalar t, int stepIdx)
    { time_ = t; timeStepIdx_ = stepIdx; }
 
    Scalar time() const final
    { return time_; }
 
    Scalar endTime() const
    { return endTime_; }
 
    void setEndTime(Scalar t)
    {
        endTime_ = t;
        if (verbose_)
            std::cout << Fmt::format("Set new end time to t = {:.5g} seconds.\n", t);
    }
 
    double wallClockTime() const
    { return timer_.elapsed(); }
 
    void setTimeStepSize(Scalar dt) final
    {
        using std::min;
        timeStepSize_ = min(dt, maxTimeStepSize());
        if (!finished() && Dune::FloatCmp::le(timeStepSize_, 0.0, 1e-14*endTime_))
            std::cerr << Fmt::format("You have set a very small timestep size (dt = {:.5g}).", timeStepSize_)
                      << " This might lead to numerical problems!\n";
    }
 
    void setMaxTimeStepSize(Scalar maxDt)
    {
        userSetMaxTimeStepSize_ = maxDt;
        setTimeStepSize(timeStepSize_);
    }
 
    Scalar timeStepSize() const final
    { return timeStepSize_; }
 
    int timeStepIndex() const
    { return timeStepIdx_; }
 
    Scalar previousTimeStepSize() const
    { return previousTimeStepSize_; }
 
    void setFinished(bool finished = true)
    { finished_ = finished; }
 
    bool finished() const override
    {
        return finished_ || endTime_-time_ < 1e-10*time_;
    }
 
    bool willBeFinished() const
    {
        return finished() || endTime_-time_-timeStepSize_ < 1e-10*timeStepSize_;
    }
 
    Scalar maxTimeStepSize() const override
    {
        if (finished())
            return 0.0;
 
        using std::min; using std::max;
        return min(userSetMaxTimeStepSize_, max<Scalar>(0.0, endTime_ - time_));
    }
 
    void reportTimeStep() const
    {
        if (verbose_)
        {
            const auto cpuTime = wallClockTime();
            using std::round;
            const auto percent = round( time_ / endTime_ * 100 );
            std::cout << Fmt::format("[{:3.0f}%] ", percent)
                      << Fmt::format("Time step {} done in {:.2g} seconds. ", timeStepIdx_, timeStepWallClockTime_)
                      << Fmt::format("Wall clock time: {:.5g}, time: {:.5g}, time step size: {:.5g}\n", cpuTime, time_, previousTimeStepSize_);
        }
    }
 
    template< class Communicator = Dune::Communication<typename Dune::MPIHelper::MPICommunicator> >
    void finalize(const Communicator& comm = Dune::MPIHelper::getCommunication())
    {
        auto cpuTime = timer_.stop();
 
        if (verbose_)
            std::cout << Fmt::format("Simulation took {:.5g} seconds on {} processes.\n", cpuTime, comm.size());
 
        if (comm.size() > 1)
            cpuTime = comm.sum(cpuTime);
 
        if (verbose_)
            std::cout << Fmt::format("The cumulative CPU time was {:.5g} seconds.\n", cpuTime);
    }
 
    bool verbose() const
    { return verbose_; }
 
    void setVerbose(bool verbose = true)
    { verbose_ = verbose; }
 
    /*
     * @}
     */
 
protected:
    Dune::Timer timer_;
    Scalar time_;
    Scalar endTime_;
 
    Scalar timeStepSize_;
    Scalar previousTimeStepSize_;
    Scalar userSetMaxTimeStepSize_;
    Scalar timeAfterLastTimeStep_, timeStepWallClockTime_;
    int timeStepIdx_;
    bool finished_;
    bool verbose_;
};
 
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;
    }
 
    void advanceTimeStep() override
    {
        const auto dt = this->timeStepSize();
        const auto newTime = this->time()+dt;
 
        // if we reached a periodic check point
        if (periodicCheckPoints_ && Dune::FloatCmp::eq(newTime - lastPeriodicCheckPoint_, deltaPeriodicCheckPoint_, 1e-7))
        {
            lastPeriodicCheckPoint_ += deltaPeriodicCheckPoint_;
            isCheckPoint_ = true;
        }
 
        // or a manually set check point
        else if (!checkPoints_.empty() && Dune::FloatCmp::eq(newTime - checkPoints_.front(), 0.0, 1e-7))
        {
            checkPoints_.pop();
            isCheckPoint_ = true;
        }
 
        // if not reset the check point flag
        else
        {
            isCheckPoint_ = false;
        }
 
        const auto previousTimeStepSize = this->previousTimeStepSize();
 
        // advance the time step like in the parent class
        TimeLoop<Scalar>::advanceTimeStep();
 
        // 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
        if (!this->willBeFinished())
        {
            using std::max;
            if (isCheckPoint_)
                this->setTimeStepSize(max(dt, previousTimeStepSize));
 
            // 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;
 
            if (periodicCheckPoints_ && Dune::FloatCmp::le(lastPeriodicCheckPoint_ + deltaPeriodicCheckPoint_ - nextTime, threshold))
                nextDt = lastPeriodicCheckPoint_ + deltaPeriodicCheckPoint_ - this->time();
 
            if (!checkPoints_.empty() && Dune::FloatCmp::le(checkPoints_.front() - nextTime, threshold))
                nextDt = checkPoints_.front() - this->time();
 
            assert(nextDt > 0.0);
            this->setTimeStepSize(nextDt);
        }
    }
 
    Scalar maxTimeStepSize() const override
    {
        using std::min;
        const auto maxCheckPointDt = computeStepSizeRespectingCheckPoints_();
        const auto maxDtParent = TimeLoop<Scalar>::maxTimeStepSize();
        return min(maxDtParent, maxCheckPointDt);
    }
 
    void setPeriodicCheckPoint(Scalar interval, Scalar offset = 0.0)
    {
        using std::signbit;
        if (signbit(interval))
            DUNE_THROW(Dune::InvalidStateException, "Interval has to be positive!");
 
        periodicCheckPoints_ = true;
        deltaPeriodicCheckPoint_ = interval;
        lastPeriodicCheckPoint_ = offset;
        while (lastPeriodicCheckPoint_ + interval - this->time() < 1e-14*interval)
            lastPeriodicCheckPoint_ += interval;
 
        if (this->verbose())
            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);
 
        // check if the current time point is a check point
        if (Dune::FloatCmp::eq(this->time()-lastPeriodicCheckPoint_, 0.0, 1e-7))
            isCheckPoint_ = true;
 
        // make sure we respect this check point on the next time step
        this->setTimeStepSize(this->timeStepSize());
    }
 
    void disablePeriodicCheckPoints()
    { periodicCheckPoints_ = false; }
 
    void removeAllCheckPoints()
    {
        periodicCheckPoints_ = false;
        while (!checkPoints_.empty())
            checkPoints_.pop();
    }
 
    bool isCheckPoint() const
    { return isCheckPoint_; }
 
    void setCheckPoint(Scalar t)
    {
        // set the check point
        setCheckPoint_(t);
 
        // make sure we respect this check point on the next time step
        this->setTimeStepSize(this->timeStepSize());
    }
 
    void setCheckPoint(const std::vector<Scalar>& checkPoints)
    { setCheckPoint(checkPoints.begin(), checkPoints.end()); }
 
    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
        this->setTimeStepSize(this->timeStepSize());
    }
 
private:
    void setCheckPoint_(Scalar t)
    {
        if (Dune::FloatCmp::le(t - this->time(), 0.0, this->timeStepSize()*1e-7))
        {
            if (this->verbose())
                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());
            return;
        }
 
        if (!checkPoints_.empty())
        {
            if (t < checkPoints_.back())
            {
                if (this->verbose())
                    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;
            }
        }
 
        checkPoints_.push(t);
        if (this->verbose())
            std::cout << Fmt::format("Set check point at t = {:.5g} seconds.\n", t);
    }
 
    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
 
#endif