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 <dune/common/fvector.hh>

#include <opm/models/discretization/common/fvbaseparameters.hh>
#include <opm/models/discretization/common/fvbaseproperties.hh>
#include <opm/models/discretization/common/restrictprolong.hh>

#include <opm/simulators/flow/NewtonIterationContext.hpp>

#include <opm/models/io/vtkmultiwriter.hh>
#include <opm/models/io/restart.hpp>

#include <opm/models/utils/simulatorutils.hpp>

#include <functional>
#include <iomanip>
#include <iostream>
#include <limits>
#include <memory>
#include <string>

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 constexpr auto 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:
    explicit 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)
    {
        // calculate the bounding box of the local partition of the grid view
        for (const auto& vertex : vertices(gridView_)) {
            for (unsigned i = 0; i < dim; ++i) {
                boundingBoxMin_[i] = std::min(boundingBoxMin_[i], vertex.geometry().corner(0)[i]);
                boundingBoxMax_[i] = std::max(boundingBoxMax_[i], vertex.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_()) {
            const 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.
            const 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");
            }

            defaultVtkWriter_ = std::make_unique<VtkMultiWriter>(asyncVtkOutput,
                                                                 gridView_,
                                                                 asImp_().outputDir(),
                                                                 asImp_().name());
        }
    }

    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; }

     unsigned intensiveQuantityHistorySize() const
     {
         if (Parameters::Get<Parameters::EnableStorageCache>() && asImp_().recycleFirstIterationStorage()) {
             return 1;
         }
         return 2;
     }

    std::string outputDir() const
    {
        return simulatorOutputDir();
    }

    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::function<void(const std::string&,
                                                            const std::string&)>,
                                         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()
    {
        elementMapper_.update(gridView_);
        vertexMapper_.update(gridView_);

        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();

        const Scalar executionTime = executionTimer.realTimeElapsed();
        const Scalar setupTime = simulator().setupTimer().realTimeElapsed();
        const Scalar prePostProcessTime = simulator().prePostProcessTimer().realTimeElapsed();
        const Scalar localCpuTime = executionTimer.cpuTimeElapsed();
        const Scalar globalCpuTime = executionTimer.globalCpuTimeElapsed();
        const Scalar writeTime = simulator().writeTimer().realTimeElapsed();
        const Scalar linearizeTime = simulator().linearizeTimer().realTimeElapsed();
        const Scalar solveTime = simulator().solveTimer().realTimeElapsed();
        const Scalar updateTime = simulator().updateTimer().realTimeElapsed();
        const unsigned numProcesses = static_cast<unsigned>(this->gridView().comm().size());
        const 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"
                      << humanReadableTime(setupTime)
                      << ", " << setupTime / (executionTime + setupTime) * 100 << "%\n"
                      << "Simulation time: " << executionTime << " seconds"
                      << humanReadableTime(executionTime)
                      << ", " << executionTime / (executionTime + setupTime) * 100 << "%\n"
                      << "    Linearization time: " << linearizeTime << " seconds"
                      << humanReadableTime(linearizeTime)
                      << ", " << linearizeTime / executionTime * 100 << "%\n"
                      << "    Linear solve time: "  << solveTime << " seconds"
                      << humanReadableTime(solveTime)
                      << ", " << solveTime / executionTime * 100 << "%\n"
                      << "    Newton update time: "  << updateTime << " seconds"
                      << humanReadableTime(updateTime)
                      << ", " << updateTime / executionTime * 100 << "%\n"
                      << "    Pre/postprocess time: "  << prePostProcessTime << " seconds"
                      << humanReadableTime(prePostProcessTime)
                      << ", " << prePostProcessTime / executionTime * 100 << "%\n"
                      << "    Output write time: "  << writeTime << " seconds"
                      << humanReadableTime(writeTime)
                      << ", " << writeTime / executionTime * 100 << "%\n"
                      << "First process' simulation CPU time: "  << localCpuTime << " seconds"
                      <<  humanReadableTime(localCpuTime) << "\n"
                      << "Number of processes: " << numProcesses << "\n"
                      << "Threads per processes: " << threadsPerProcess << "\n"
                      << "Total CPU time: " << globalCpuTime << " seconds"
                      << humanReadableTime(globalCpuTime) << "\n"
                      << "\n"
                      << "----------------------------------------------------------------\n"
                      << std::endl;
        }
    }

    void timeIntegration()
    {
        const unsigned maxFails = asImp_().maxTimeIntegrationFailures();
        Scalar minTimeStep = asImp_().minTimeStepSize();

        std::string errorMessage;
        for (unsigned i = 0; i < maxFails; ++i) {
            bool converged = model().update();
            if (converged) {
                return;
            }

            const Scalar dt = simulator().timeStepSize();
            Scalar nextDt = dt / 2.0;
            if (dt < minTimeStep * (1.0 + 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(minTimeStep)) + " seconds. Giving up!";
                    break; // give up: we can't make the time step smaller anymore!
                }
            }
            else if (nextDt < minTimeStep) {
                nextDt = minTimeStep;
            }
            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
        const Scalar t = simulator().time() + simulator().timeStepSize();

        defaultVtkWriter_->beginWrite(t);
        model().prepareOutputFields();
        model().appendOutputFields(*defaultVtkWriter_);
        defaultVtkWriter_->endWrite(false);
    }

    VtkMultiWriter& defaultVtkWriter() const
    { return *defaultVtkWriter_; }

    const NewtonIterationContext& iterationContext() const
    { return iterationContext_; }

    void resetIterationForNewTimestep()
    { iterationContext_.resetForNewTimestep(); }

    void advanceIteration()
    { iterationContext_.advanceIteration(); }

    void markTimestepInitialized()
    { iterationContext_.markTimestepInitialized(); }

private:
    // LocalContextGuard needs mutable access to save/restore the context
    // during NLDD domain-local solves.
    template<class P> friend class LocalContextGuard;

    NewtonIterationContext& mutableIterationContext()
    { return iterationContext_; }

public:

protected:
    Scalar nextTimeStepSize_;
    NewtonIterationContext iterationContext_;

    bool enableVtkOutput_() const
    { return Parameters::Get<Parameters::EnableVtkOutput>(); }

private:
    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_;
    std::unique_ptr<VtkMultiWriter> defaultVtkWriter_{};
};

} // namespace Opm

#endif