fvbaseproblem.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:
/*
  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/>.
 
  Consult the COPYING file in the top-level source directory of this
  module for the precise wording of the license and the list of
  copyright holders.
*/
#ifndef EWOMS_FV_BASE_PROBLEM_HH
#define EWOMS_FV_BASE_PROBLEM_HH
 
#include <opm/models/discretization/common/fvbaseparameters.hh>
#include <opm/models/discretization/common/fvbaseproperties.hh>
 
#include <opm/models/io/vtkmultiwriter.hh>
#include <opm/models/io/restart.hh>
#include <opm/models/discretization/common/restrictprolong.hh>
 
#include <dune/common/fvector.hh>
 
#include <iostream>
#include <limits>
#include <string>
 
#include <sys/stat.h>
 
namespace Opm::Properties {
 
template <class TypeTag, class MyTypeTag>
struct NewtonMethod;
 
} // namespace Opm::Properties
 
namespace Opm {
 
template<class TypeTag>
class FvBaseProblem
{
private:
    using Implementation = GetPropType<TypeTag, Properties::Problem>;
    using GridView = GetPropType<TypeTag, Properties::GridView>;
 
    static const int vtkOutputFormat = getPropValue<TypeTag, Properties::VtkOutputFormat>();
    using VtkMultiWriter = ::Opm::VtkMultiWriter<GridView, vtkOutputFormat>;
 
    using Model = GetPropType<TypeTag, Properties::Model>;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using Simulator = GetPropType<TypeTag, Properties::Simulator>;
    using ThreadManager = GetPropType<TypeTag, Properties::ThreadManager>;
    using NewtonMethod = GetPropType<TypeTag, Properties::NewtonMethod>;
 
    using VertexMapper = GetPropType<TypeTag, Properties::VertexMapper>;
    using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
 
    using RateVector = GetPropType<TypeTag, Properties::RateVector>;
    using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
    using Constraints = GetPropType<TypeTag, Properties::Constraints>;
 
    enum {
        dim = GridView::dimension,
        dimWorld = GridView::dimensionworld
    };
 
    using Element = typename GridView::template Codim<0>::Entity;
    using Vertex = typename GridView::template Codim<dim>::Entity;
    using VertexIterator = typename GridView::template Codim<dim>::Iterator;
 
    using CoordScalar = typename GridView::Grid::ctype;
    using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
 
public:
    // the default restriction and prolongation for adaptation is simply an empty one
    using RestrictProlongOperator = EmptyRestrictProlong ;
 
private:
    // copying a problem is not a good idea
    FvBaseProblem(const FvBaseProblem& ) = delete;
 
public:
    FvBaseProblem(Simulator& simulator)
        : nextTimeStepSize_(0.0)
        , gridView_(simulator.gridView())
        , elementMapper_(gridView_, Dune::mcmgElementLayout())
        , vertexMapper_(gridView_, Dune::mcmgVertexLayout())
        , boundingBoxMin_(std::numeric_limits<double>::max())
        , boundingBoxMax_(-std::numeric_limits<double>::max())
        , simulator_(simulator)
        , defaultVtkWriter_(0)
    {
        // calculate the bounding box of the local partition of the grid view
        VertexIterator vIt = gridView_.template begin<dim>();
        const VertexIterator vEndIt = gridView_.template end<dim>();
        for (; vIt!=vEndIt; ++vIt) {
            for (unsigned i=0; i<dim; i++) {
                boundingBoxMin_[i] = std::min(boundingBoxMin_[i], vIt->geometry().corner(0)[i]);
                boundingBoxMax_[i] = std::max(boundingBoxMax_[i], vIt->geometry().corner(0)[i]);
            }
        }
 
        // communicate to get the bounding box of the whole domain
        for (unsigned i = 0; i < dim; ++i) {
            boundingBoxMin_[i] = gridView_.comm().min(boundingBoxMin_[i]);
            boundingBoxMax_[i] = gridView_.comm().max(boundingBoxMax_[i]);
        }
 
        if (enableVtkOutput_()) {
            bool asyncVtkOutput =
                simulator_.gridView().comm().size() == 1 &&
                Parameters::Get<Parameters::EnableAsyncVtkOutput>();
 
            // asynchonous VTK output currently does not work in conjunction with grid
            // adaptivity because the async-IO code assumes that the grid stays
            // constant. complain about that case.
            bool enableGridAdaptation = Parameters::Get<Parameters::EnableGridAdaptation>();
            if (asyncVtkOutput && enableGridAdaptation)
                throw std::runtime_error("Asynchronous VTK output currently cannot be used "
                                         "at the same time as grid adaptivity");
 
            std::string outputDir = asImp_().outputDir();
 
            defaultVtkWriter_ =
                new VtkMultiWriter(asyncVtkOutput, gridView_, outputDir, asImp_().name());
        }
    }
 
    ~FvBaseProblem()
    { delete defaultVtkWriter_; }
 
    static void registerParameters()
    {
        Model::registerParameters();
        Parameters::Register<Parameters::MaxTimeStepSize<Scalar>>
            ("The maximum size to which all time steps are limited to [s]");
        Parameters::Register<Parameters::MinTimeStepSize<Scalar>>
            ("The minimum size to which all time steps are limited to [s]");
        Parameters::Register<Parameters::MaxTimeStepDivisions>
            ("The maximum number of divisions by two of the timestep size "
             "before the simulation bails out");
        Parameters::Register<Parameters::EnableAsyncVtkOutput>
            ("Dispatch a separate thread to write the VTK output");
        Parameters::Register<Parameters::ContinueOnConvergenceError>
            ("Continue with a non-converged solution instead of giving up "
             "if we encounter a time step size smaller than the minimum time "
             "step size.");
    }
 
    bool recycleFirstIterationStorage() const
    { return true; }
 
    std::string outputDir() const
    {
        std::string outputDir = Parameters::Get<Parameters::OutputDir>();
 
        if (outputDir.empty())
            outputDir = ".";
 
        // TODO: replace this by std::filesystem once we require c++-2017
        struct stat st;
        if (::stat(outputDir.c_str(), &st) != 0)
            throw std::runtime_error("Could not access output directory '"+outputDir+"':"
                                     +strerror(errno));
        if (!S_ISDIR(st.st_mode))
            throw std::runtime_error("Path to output directory '"+outputDir+"' exists but is not a directory");
 
        if (access(outputDir.c_str(), W_OK) != 0)
            throw std::runtime_error("Output directory '"+outputDir+"' exists but is not writeable");
 
        return outputDir;
    }
 
    static std::string helpPreamble(int,
                                    const char **argv)
    {
        std::string desc = Implementation::briefDescription();
        if (!desc.empty())
            desc = desc + "\n";
 
        return
            "Usage: "+std::string(argv[0]) + " [OPTIONS]\n"
            + desc;
    }
 
    static std::string briefDescription()
    { return ""; }
 
    // TODO (?): detailedDescription()
 
    static int handlePositionalParameter(std::set<std::string>&,
                                         std::string& errorMsg,
                                         int,
                                         const char** argv,
                                         int paramIdx,
                                         int)
    {
        errorMsg = std::string("Illegal parameter \"")+argv[paramIdx]+"\".";
        return 0;
    }
 
    void finishInit()
    { }
 
    void prefetch(const Element&) const
    {
        // do nothing by default
    }
 
    void gridChanged()
    {
#if DUNE_VERSION_NEWER(DUNE_GRID, 2, 8)
                elementMapper_.update(gridView_);
                vertexMapper_.update(gridView_);
#else
                elementMapper_.update();
                vertexMapper_.update();
#endif
 
        if (enableVtkOutput_())
            defaultVtkWriter_->gridChanged();
    }
 
    template <class Context>
    void boundary(BoundaryRateVector&,
                  const Context&,
                  unsigned,
                  unsigned) const
    { throw std::logic_error("Problem does not provide a boundary() method"); }
 
    template <class Context>
    void constraints(Constraints&,
                     const Context&,
                     unsigned,
                     unsigned) const
    { throw std::logic_error("Problem does not provide a constraints() method"); }
 
    template <class Context>
    void source(RateVector&,
                const Context&,
                unsigned,
                unsigned) const
    { throw std::logic_error("Problem does not provide a source() method"); }
 
    template <class Context>
    void initial(PrimaryVariables&,
                 const Context&,
                 unsigned,
                 unsigned) const
    { throw std::logic_error("Problem does not provide a initial() method"); }
 
    template <class Context>
    Scalar extrusionFactor(const Context&,
                           unsigned,
                           unsigned) const
    { return asImp_().extrusionFactor(); }
 
    Scalar extrusionFactor() const
    { return 1.0; }
 
    void initialSolutionApplied()
    {}
 
    void beginEpisode()
    { }
 
    void beginTimeStep()
    { }
 
    void beginIteration()
    { }
 
    void endIteration()
    { }
 
    void endTimeStep()
    { }
 
    void endEpisode()
    {
        std::cerr << "The end of episode " << simulator().episodeIndex() + 1 << " has been "
                  << "reached, but the problem does not override the endEpisode() method. "
                  << "Doing nothing!\n";
    }
 
    void finalize()
    {
        const auto& executionTimer = simulator().executionTimer();
 
        Scalar executionTime = executionTimer.realTimeElapsed();
        Scalar setupTime = simulator().setupTimer().realTimeElapsed();
        Scalar prePostProcessTime = simulator().prePostProcessTimer().realTimeElapsed();
        Scalar localCpuTime = executionTimer.cpuTimeElapsed();
        Scalar globalCpuTime = executionTimer.globalCpuTimeElapsed();
        Scalar writeTime = simulator().writeTimer().realTimeElapsed();
        Scalar linearizeTime = simulator().linearizeTimer().realTimeElapsed();
        Scalar solveTime = simulator().solveTimer().realTimeElapsed();
        Scalar updateTime = simulator().updateTimer().realTimeElapsed();
        unsigned numProcesses = static_cast<unsigned>(this->gridView().comm().size());
        unsigned threadsPerProcess = ThreadManager::maxThreads();
        if (gridView().comm().rank() == 0) {
            std::cout << std::setprecision(3)
                      << "Simulation of problem '" << asImp_().name() << "' finished.\n"
                      << "\n"
                      << "------------------------ Timing ------------------------\n"
                      << "Setup time: " << setupTime << " seconds" << Simulator::humanReadableTime(setupTime)
                      << ", " << setupTime/(executionTime + setupTime)*100 << "%\n"
                      << "Simulation time: " << executionTime << " seconds" << Simulator::humanReadableTime(executionTime)
                      << ", " << executionTime/(executionTime + setupTime)*100 << "%\n"
                      << "    Linearization time: " << linearizeTime << " seconds" << Simulator::humanReadableTime(linearizeTime)
                      << ", " << linearizeTime/executionTime*100 << "%\n"
                      << "    Linear solve time: "  << solveTime << " seconds" << Simulator::humanReadableTime(solveTime)
                      << ", " << solveTime/executionTime*100 << "%\n"
                      << "    Newton update time: "  << updateTime << " seconds" << Simulator::humanReadableTime(updateTime)
                      << ", " << updateTime/executionTime*100 << "%\n"
                      << "    Pre/postprocess time: "  << prePostProcessTime << " seconds" << Simulator::humanReadableTime(prePostProcessTime)
                      << ", " << prePostProcessTime/executionTime*100 << "%\n"
                      << "    Output write time: "  << writeTime << " seconds" << Simulator::humanReadableTime(writeTime)
                      << ", " << writeTime/executionTime*100 << "%\n"
                      << "First process' simulation CPU time: "  << localCpuTime << " seconds" <<  Simulator::humanReadableTime(localCpuTime) << "\n"
                      << "Number of processes: " << numProcesses << "\n"
                      << "Threads per processes: " << threadsPerProcess << "\n"
                      << "Total CPU time: " << globalCpuTime << " seconds" << Simulator::humanReadableTime(globalCpuTime) << "\n"
                      << "\n"
                      << "----------------------------------------------------------------\n"
                      << std::endl;
        }
    }
 
    void timeIntegration()
    {
        unsigned maxFails = asImp_().maxTimeIntegrationFailures();
        Scalar minTimeStepSize = asImp_().minTimeStepSize();
 
        std::string errorMessage;
        for (unsigned i = 0; i < maxFails; ++i) {
            bool converged = model().update();
            if (converged)
                return;
 
            Scalar dt = simulator().timeStepSize();
            Scalar nextDt = dt / 2.0;
            if (dt < minTimeStepSize*(1 + 1e-9)) {
                if (asImp_().continueOnConvergenceError()) {
                    if (gridView().comm().rank() == 0)
                        std::cout << "Newton solver did not converge with minimum time step of "
                                  << dt << " seconds. Continuing with unconverged solution!\n"
                                  << std::flush;
                    return;
                }
                else {
                    errorMessage =
                        "Time integration did not succeed with the minumum time step size of "
                        + std::to_string(double(minTimeStepSize)) + " seconds. Giving up!";
                    break; // give up: we can't make the time step smaller anymore!
                }
            }
            else if (nextDt < minTimeStepSize)
                nextDt = minTimeStepSize;
            simulator().setTimeStepSize(nextDt);
 
            // update failed
            if (gridView().comm().rank() == 0)
                std::cout << "Newton solver did not converge with "
                          << "dt=" << dt << " seconds. Retrying with time step of "
                          << nextDt << " seconds\n" << std::flush;
        }
 
        if (errorMessage.empty())
            errorMessage =
                "Newton solver didn't converge after "
                +std::to_string(maxFails)+" time-step divisions. dt="
                +std::to_string(double(simulator().timeStepSize()));
        throw std::runtime_error(errorMessage);
    }
 
    Scalar minTimeStepSize() const
    { return Parameters::Get<Parameters::MinTimeStepSize<Scalar>>(); }
 
    unsigned maxTimeIntegrationFailures() const
    { return Parameters::Get<Parameters::MaxTimeStepDivisions>(); }
 
    bool continueOnConvergenceError() const
    { return Parameters::Get<Parameters::ContinueOnConvergenceError>(); }
 
    void setNextTimeStepSize(Scalar dt)
    { nextTimeStepSize_ = dt; }
 
    Scalar nextTimeStepSize() const
    {
        if (nextTimeStepSize_ > 0.0)
            return nextTimeStepSize_;
 
        Scalar dtNext = std::min(Parameters::Get<Parameters::MaxTimeStepSize<Scalar>>(),
                                 newtonMethod().suggestTimeStepSize(simulator().timeStepSize()));
 
        if (dtNext < simulator().maxTimeStepSize()
            && simulator().maxTimeStepSize() < dtNext*2)
        {
            dtNext = simulator().maxTimeStepSize()/2 * 1.01;
        }
 
        return dtNext;
    }
 
    bool shouldWriteRestartFile() const
    {
        return simulator().timeStepIndex() > 0 &&
            (simulator().timeStepIndex() % 10 == 0);
    }
 
    bool shouldWriteOutput() const
    { return true; }
 
    void advanceTimeLevel()
    { model().advanceTimeLevel(); }
 
    std::string name() const
    { return "sim"; }
 
    const GridView& gridView() const
    { return gridView_; }
 
    const GlobalPosition& boundingBoxMin() const
    { return boundingBoxMin_; }
 
    const GlobalPosition& boundingBoxMax() const
    { return boundingBoxMax_; }
 
    const VertexMapper& vertexMapper() const
    { return vertexMapper_; }
 
    const ElementMapper& elementMapper() const
    { return elementMapper_; }
 
    Simulator& simulator()
    { return simulator_; }
 
    const Simulator& simulator() const
    { return simulator_; }
 
    Model& model()
    { return simulator_.model(); }
 
    const Model& model() const
    { return simulator_.model(); }
 
    NewtonMethod& newtonMethod()
    { return model().newtonMethod(); }
 
    const NewtonMethod& newtonMethod() const
    { return model().newtonMethod(); }
    // \}
 
    RestrictProlongOperator restrictProlongOperator()
    {
        return RestrictProlongOperator();
    }
 
    unsigned markForGridAdaptation()
    {
        return 0;
    }
 
    template <class Restarter>
    void serialize(Restarter& res)
    {
        if (enableVtkOutput_())
            defaultVtkWriter_->serialize(res);
    }
 
    template <class Restarter>
    void deserialize(Restarter& res)
    {
        if (enableVtkOutput_())
            defaultVtkWriter_->deserialize(res);
    }
 
    void writeOutput(bool verbose = true)
    {
        if (!enableVtkOutput_())
            return;
 
        if (verbose && gridView().comm().rank() == 0)
            std::cout << "Writing visualization results for the current time step.\n"
                      << std::flush;
 
        // calculate the time _after_ the time was updated
        Scalar t = simulator().time() + simulator().timeStepSize();
 
        defaultVtkWriter_->beginWrite(t);
        model().prepareOutputFields();
        model().appendOutputFields(*defaultVtkWriter_);
        defaultVtkWriter_->endWrite();
 
    }
 
    VtkMultiWriter& defaultVtkWriter() const
    { return defaultVtkWriter_; }
 
protected:
    Scalar nextTimeStepSize_;
 
private:
    bool enableVtkOutput_() const
    { return Parameters::Get<Parameters::EnableVtkOutput>(); }
 
    Implementation& asImp_()
    { return *static_cast<Implementation *>(this); }
 
    const Implementation& asImp_() const
    { return *static_cast<const Implementation *>(this); }
 
    // Grid management stuff
    const GridView gridView_;
    ElementMapper elementMapper_;
    VertexMapper vertexMapper_;
    GlobalPosition boundingBoxMin_;
    GlobalPosition boundingBoxMax_;
 
    // Attributes required for the actual simulation
    Simulator& simulator_;
    mutable VtkMultiWriter *defaultVtkWriter_;
};
 
} // namespace Opm
 
#endif