simulator.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:
/*
  Copyright (C) 2009-2013 by Andreas Lauser
  Copyright (C) 2011-2012 by Markus Wolff
  Copyright (C) 2011 by Benjamin Faigle

  This file is part of the Open Porous Media project (OPM).

  OPM 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 2 of the License, or
  (at your option) any later version.

  OPM 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 OPM.  If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef EWOMS_SIMULATOR_HH
#define EWOMS_SIMULATOR_HH

#include <ewoms/io/restart.hh>
#include <ewoms/common/parametersystem.hh>

#include <ewoms/common/propertysystem.hh>
#include <ewoms/common/timer.hh>

#include <dune/common/version.hh>
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2,3)
#include <dune/common/parallel/mpihelper.hh>
#else
#include <dune/common/mpihelper.hh>
#endif

#include <iostream>
#include <iomanip>
#include <memory>

namespace Ewoms {
namespace Properties {
NEW_PROP_TAG(Scalar);
NEW_PROP_TAG(GridManager);
NEW_PROP_TAG(GridView);
NEW_PROP_TAG(Model);
NEW_PROP_TAG(Problem);
NEW_PROP_TAG(EndTime);
NEW_PROP_TAG(RestartTime);
NEW_PROP_TAG(InitialTimeStepSize);
}

template <class TypeTag>
class Simulator
{
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, GridManager) GridManager;
    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
    typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;


public:
    // do not allow to copy simulators around
    Simulator(const Simulator &) = delete;

    Simulator(bool verbose = true)
    {
        setupTimer_.start();

        verbose_ = verbose && Dune::MPIHelper::getCollectiveCommunication().rank() == 0;

        timeStepIdx_ = 0;
        startTime_ = 0.0;
        time_ = 0.0;
        endTime_ = EWOMS_GET_PARAM(TypeTag, Scalar, EndTime);
        timeStepSize_ = EWOMS_GET_PARAM(TypeTag, Scalar, InitialTimeStepSize);

        episodeIdx_ = 0;
        episodeStartTime_ = 0;
        episodeLength_ = 1e100;

        finished_ = false;

        if (verbose_)
            std::cout << "Allocating the grid\n" << std::flush;
        gridManager_.reset(new GridManager(*this));

        if (verbose_)
            std::cout << "Distributing the grid\n" << std::flush;
        gridManager_->loadBalance();

        if (verbose_)
            std::cout << "Allocating the model\n" << std::flush;
        model_.reset(new Model(*this));

        if (verbose_)
            std::cout << "Allocating the problem\n" << std::flush;
        problem_.reset(new Problem(*this));

        if (verbose_)
            std::cout << "Finish init of the model\n" << std::flush;
        model_->finishInit();

        if (verbose_)
            std::cout << "Finish init of the problem\n" << std::flush;
        problem_->finishInit();

        if (verbose_)
            std::cout << "Construction of simulation done\n" << std::flush;
    }

    static void registerParameters()
    {
        EWOMS_REGISTER_PARAM(TypeTag, Scalar, EndTime,
                             "The simulation time at which the simulation is finished [s]");
        EWOMS_REGISTER_PARAM(TypeTag, Scalar, InitialTimeStepSize,
                             "The size of the initial time step [s]");
        EWOMS_REGISTER_PARAM(TypeTag, Scalar, RestartTime,
                             "The simulation time at which a restart should be attempted [s]");

        GridManager::registerParameters();
        Model::registerParameters();
        Problem::registerParameters();
    }

    GridManager &gridManager()
    { return *gridManager_; }

    const GridManager &gridManager() const
    { return *gridManager_; }

    GridView &gridView()
    { return gridManager_->gridView(); }

    const GridView &gridView() const
    { return gridManager_->gridView(); }

    Model &model()
    { return *model_; }

    const Model &model() const
    { return *model_; }

    Problem &problem()
    { return *problem_; }

    const Problem &problem() const
    { return *problem_; }

    void setStartTime(Scalar t)
    { startTime_ = t; }

    Scalar startTime() const
    { return startTime_; }

    void setTime(Scalar t)
    { time_ = t; }

    void setTime(Scalar t, int stepIdx)
    {
        time_ = t;
        timeStepIdx_ = stepIdx;
    }

    Scalar time() const
    { return time_; }

    void setEndTime(Scalar t)
    { endTime_ = t; }

    Scalar endTime() const
    { return endTime_; }

    const Ewoms::Timer& timer() const
    { return executionTimer_; }

    Scalar totalWriteTime() const
    { return totalWriteTime_; }

    double setupTime() const
    { return setupTimer_.realTimeElapsed(); }

    void setTimeStepSize(Scalar timeStepSize)
    { timeStepSize_ = timeStepSize; }

    Scalar timeStepSize() const
    {
        Scalar maximumTimeStepSize =
            std::min(episodeMaxTimeStepSize(),
                     std::max<Scalar>(0.0, endTime() - this->time()));

        return std::min(timeStepSize_, maximumTimeStepSize);
    }

    int timeStepIndex() const
    { return timeStepIdx_; }

    void setFinished(bool yesno = true)
    { finished_ = yesno; }

    bool finished() const
    {
        assert(timeStepSize_ >= 0.0);
        Scalar eps = std::max(std::abs(this->time()), timeStepSize())*1e-8;
        return finished_ || (this->time() + eps >= endTime());
    }

    bool willBeFinished() const
    {
        return
            finished_ ||
            (this->time() + std::max(std::abs(this->time()), timeStepSize())*1e-8 + timeStepSize_
             >= endTime());
    }

    Scalar maxTimeStepSize() const
    {
        if (finished())
            return 0.0;

        return std::min(episodeMaxTimeStepSize(),
                        std::max<Scalar>(0.0, endTime() - this->time()));
    }

    void startNextEpisode(Scalar episodeStartTime, Scalar episodeLength)
    {
        ++episodeIdx_;
        episodeStartTime_ = episodeStartTime;
        episodeLength_ = episodeLength;
    }

    void startNextEpisode(Scalar len = 1e100)
    {
        ++episodeIdx_;
        episodeStartTime_ = startTime_ + time_;
        episodeLength_ = len;
    }

    void setEpisodeIndex(int episodeIdx)
    { episodeIdx_ = episodeIdx; }

    int episodeIndex() const
    { return episodeIdx_; }

    Scalar episodeStartTime() const
    { return episodeStartTime_; }

    Scalar episodeLength() const
    { return episodeLength_; }

    bool episodeIsOver() const
    { return startTime() + this->time() > episodeStartTime_ + episodeLength()*(1 - 1e-8); }

    bool episodeWillBeOver() const
    {
        return startTime() + this->time() + timeStepSize()
            >=  episodeStartTime_ + episodeLength()*(1 - 1e-8);
    }

    Scalar episodeMaxTimeStepSize() const
    {
        // if the current episode is over and the simulation
        // wants to give it some extra time, we will return
        // the time step size it suggested instead of trying
        // to align it to the end of the episode.
        if (episodeIsOver())
            return 0.0;

        // make sure that we don't exceed the end of the
        // current episode.
        return std::max<Scalar>(0.0,
                                (episodeStartTime() + episodeLength())
                                - (this->time() + this->startTime()));
    }

    /*
     * \}
     */

    void run()
    {
        Scalar restartTime = EWOMS_GET_PARAM(TypeTag, Scalar, RestartTime);
        if (restartTime > -1e100) {
            // try to restart a previous simulation
            time_ = restartTime;

            Ewoms::Restart res;
            res.deserializeBegin(*this, time_);
            if (verbose_)
                std::cout << "Deserialize from file '" << res.fileName() << "'\n" << std::flush;
            this->deserialize(res);
            problem_->deserialize(res);
            model_->deserialize(res);
            res.deserializeEnd();
            if (verbose_)
                std::cout << "Deserialization done."
                          << " Simulator time: " << time() << humanReadableTime(time())
                          << " Time step index: " << timeStepIndex()
                          << " Episode index: " << episodeIndex()
                          << "\n" << std::flush;
        }
        else {
            // if no restart is done, apply the initial solution
            if (verbose_)
                std::cout << "Applying the initial solution of the \"" << problem_->name()
                          << "\" problem\n" << std::flush;

            Scalar oldTimeStepSize = timeStepSize_;
            int oldTimeStepIdx = timeStepIdx_;
            timeStepSize_ = 0.0;
            timeStepIdx_ = -1;

            model_->applyInitialSolution();

            // write initial condition
            if (problem_->shouldWriteOutput())
                problem_->writeOutput();

            timeStepSize_ = oldTimeStepSize;
            timeStepIdx_ = oldTimeStepIdx;
        }

        setupTimer_.stop();

        totalWriteTime_ = 0.0;
        executionTimer_.start();

        prePostProcessTime_ = 0;
        Timer prePostProcessTimer;
        bool episodeBegins = episodeIsOver() || (timeStepIdx_ == 0);
        // do the time steps
        while (!finished()) {
            prePostProcessTimer.start();
            if (episodeBegins) {
                // notify the problem that a new episode has just been
                // started.
                problem_->beginEpisode();

                if (finished()) {
                    // the problem can chose to terminate the simulation in
                    // beginEpisode(), so we have handle this case.
                    problem_->endEpisode();
                    prePostProcessTimer.stop();
                    prePostProcessTime_ += prePostProcessTimer.realTimeElapsed();
                    break;
                }
            }
            episodeBegins = false;

            if (verbose_) {
                std::cout << "Begin time step " << timeStepIndex() << ". "
                          << "Start time: " << this->time() << " seconds" << humanReadableTime(this->time())
                          << ", step size: " << timeStepSize() << " seconds" << humanReadableTime(timeStepSize())
                          << "\n";
            }

            // pre-process the current solution
            problem_->beginTimeStep();
            if (finished()) {
                // the problem can chose to terminate the simulation in
                // beginTimeStep(), so we have handle this case.
                problem_->endTimeStep();
                problem_->endEpisode();
                prePostProcessTimer.stop();
                prePostProcessTime_ += prePostProcessTimer.realTimeElapsed();
                break;
            }
            prePostProcessTimer.stop();
            prePostProcessTime_ += prePostProcessTimer.realTimeElapsed();

            // execute the time integration scheme
            problem_->timeIntegration();

            // post-process the current solution
            prePostProcessTimer.start();
            problem_->endTimeStep();
            prePostProcessTimer.stop();
            prePostProcessTime_ += prePostProcessTimer.realTimeElapsed();

            // write the result to disk
            writeTimer_.start();
            if (problem_->shouldWriteOutput())
                problem_->writeOutput();
            writeTimer_.stop();
            totalWriteTime_ += writeTimer_.realTimeElapsed();

            // do the next time integration
            Scalar oldDt = timeStepSize();
            problem_->advanceTimeLevel();

            // advance the simulated time by the current time step size
            time_ += oldDt;
            ++timeStepIdx_;

            if (verbose_) {
                std::cout << "Time step " << timeStepIndex() << " done. "
                          << "CPU time: " << executionTimer_.realTimeElapsed() << " seconds" << humanReadableTime(executionTimer_.realTimeElapsed())
                          << ", end time: " << this->time() << " seconds" << humanReadableTime(this->time())
                          << ", step size: " << oldDt << " seconds" << humanReadableTime(oldDt)
                          << "\n" << std::flush;
            }

            prePostProcessTimer.start();
            // notify the problem if an episode is finished
            if (episodeIsOver()) {
                // make the linearization not usable for the first iteration of the next
                // time step at the end of each episode. Strictly speaking, this is a
                // layering violation as the simulator is not supposed to know any
                // specifics of the non-linear solver, but this call is too easy to
                // forget within the problem's endEpisode(), so we do it here anyway!
                model_->linearizer().setLinearizationReusable(false);

                // Notify the problem about the end of the current episode...
                problem_->endEpisode();
                episodeBegins = true;
            }
            else
                // ask the problem to provide the next time step size
                setTimeStepSize(problem_->nextTimeStepSize());
            prePostProcessTimer.stop();
            prePostProcessTime_ += prePostProcessTimer.realTimeElapsed();

            // write restart file if mandated by the problem
            writeTimer_.start();
            if (problem_->shouldWriteRestartFile())
                serialize();
            writeTimer_.stop();
            totalWriteTime_ += writeTimer_.realTimeElapsed();
        }

        executionTimer_.stop();

        problem_->finalize();
    }

    Scalar prePostProcessTime() const
    { return prePostProcessTime_; }

    void addPrePostProcessTime(Scalar value)
    { prePostProcessTime_ += value; }

    static std::string humanReadableTime(Scalar timeInSeconds, bool isAmendment=true)
    {
        std::ostringstream oss;
        oss << std::setprecision(4);
        if (isAmendment)
            oss << " (";
        if (timeInSeconds >= 365.25*24*60*60) {
            int years = timeInSeconds/(365.25*24*60*60);
            int days = (timeInSeconds - years*(365.25*24*60*60))/(24*60*60);
            double hours = (timeInSeconds - years*(365.25*24*60*60) - days*(24*60*60))/(60*60);
            oss << years << " years, " << days << " days, "  << hours << " hours";
        }
        else if (timeInSeconds >= 24.0*60*60) {
            int days = timeInSeconds/(24*60*60);
            int hours = (timeInSeconds - days*(24*60*60))/(60*60);
            double minutes = (timeInSeconds - days*(24*60*60) - hours*(60*60))/60;
            oss << days << " days, " << hours << " hours, " << minutes << " minutes";
        }
        else if (timeInSeconds >= 60.0*60) {
            int hours = timeInSeconds/(60*60);
            int minutes = (timeInSeconds - hours*(60*60))/60;
            double seconds = (timeInSeconds - hours*(60*60) - minutes*60);
            oss << hours << " hours, " << minutes << " minutes, "  << seconds << " seconds";
        }
        else if (timeInSeconds >= 60.0) {
            int minutes = timeInSeconds/60;
            double seconds = (timeInSeconds - minutes*60);
            oss << minutes << " minutes, "  << seconds << " seconds";
        }
        else if (!isAmendment)
            oss << timeInSeconds << " seconds";
        else
            return "";
        if (isAmendment)
            oss << ")";

        return oss.str();
    }

    void serialize()
    {
        typedef Ewoms::Restart Restarter;
        Restarter res;
        res.serializeBegin(*this);
        if (gridView().comm().rank() == 0)
            std::cout << "Serialize to file '" << res.fileName() << "'"
                      << ", next time step size: " << timeStepSize()
                      << "\n" << std::flush;

        this->serialize(res);
        problem_->serialize(res);
        model_->serialize(res);
        res.serializeEnd();
    }

    template <class Restarter>
    void serialize(Restarter &restarter)
    {
        restarter.serializeSectionBegin("Simulator");
        restarter.serializeStream()
            << episodeIdx_ << " "
            << episodeStartTime_ << " "
            << episodeLength_ << " "
            << startTime_ << " "
            << time_ << " "
            << timeStepIdx_ << " ";
        restarter.serializeSectionEnd();
    }

    template <class Restarter>
    void deserialize(Restarter &restarter)
    {
        restarter.deserializeSectionBegin("Simulator");
        restarter.deserializeStream()
            >> episodeIdx_
            >> episodeStartTime_
            >> episodeLength_
            >> startTime_
            >> time_
            >> timeStepIdx_;
        restarter.deserializeSectionEnd();
    }

private:
    std::unique_ptr<GridManager> gridManager_;
    std::unique_ptr<Model> model_;
    std::unique_ptr<Problem> problem_;

    int episodeIdx_;
    Scalar episodeStartTime_;
    Scalar episodeLength_;

    Ewoms::Timer setupTimer_;
    Ewoms::Timer executionTimer_;
    Ewoms::Timer writeTimer_;
    Scalar totalWriteTime_;
    Scalar startTime_;
    Scalar prePostProcessTime_;
    Scalar time_;
    Scalar endTime_;

    Scalar timeStepSize_;
    int timeStepIdx_;
    bool finished_;
    bool verbose_;
};
} // namespace Ewoms

#endif