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:
/*
  Copyright (C) 2008-2013 by Andreas Lauser
  Copyright (C) 2010-2011 by Markus Wolff

  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_FV_BASE_PROBLEM_HH
#define EWOMS_FV_BASE_PROBLEM_HH

#include "fvbaseproperties.hh"

#include <ewoms/io/vtkmultiwriter.hh>
#include <ewoms/io/restart.hh>
#include <dune/common/fvector.hh>

#include <iostream>
#include <limits>
#include <string>

namespace Ewoms {

template<class TypeTag>
class FvBaseProblem
{
private:
    typedef typename GET_PROP_TYPE(TypeTag, Problem) Implementation;
    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;

    static const int vtkOutputFormat = GET_PROP_VALUE(TypeTag, VtkOutputFormat);
    typedef Ewoms::VtkMultiWriter<GridView, vtkOutputFormat> VtkMultiWriter;

    typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
    typedef typename GET_PROP_TYPE(TypeTag, ThreadManager) ThreadManager;
    typedef typename GET_PROP_TYPE(TypeTag, NewtonMethod) NewtonMethod;

    typedef typename GET_PROP_TYPE(TypeTag, VertexMapper) VertexMapper;
    typedef typename GET_PROP_TYPE(TypeTag, ElementMapper) ElementMapper;

    typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
    typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector;
    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
    typedef typename GET_PROP_TYPE(TypeTag, Constraints) Constraints;

    enum {
        dim = GridView::dimension,
        dimWorld = GridView::dimensionworld
    };

    typedef typename GridView::template Codim<0>::Entity Element;
    typedef typename GridView::template Codim<dim>::Entity Vertex;
    typedef typename GridView::template Codim<dim>::Iterator VertexIterator;

    typedef typename GridView::Grid::ctype CoordScalar;
    typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;

    // copying a problem is not a good idea
    FvBaseProblem(const FvBaseProblem &) = delete;

public:
    FvBaseProblem(Simulator &simulator)
        : gridView_(simulator.gridView())
        , elementMapper_(gridView_)
        , vertexMapper_(gridView_)
        , 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 (int 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 (int i = 0; i < dim; ++i) {
            boundingBoxMin_[i] = gridView_.comm().min(boundingBoxMin_[i]);
            boundingBoxMax_[i] = gridView_.comm().max(boundingBoxMax_[i]);
        }

        if (enableVtkOutput_())
            defaultVtkWriter_ = new VtkMultiWriter(gridView_, asImp_().name());
    }

    static void registerParameters()
    {
        Model::registerParameters();
        EWOMS_REGISTER_PARAM(TypeTag, Scalar, MaxTimeStepSize,
                             "The maximum size to which all time steps are limited to [s]");
        EWOMS_REGISTER_PARAM(TypeTag, Scalar, MinTimeStepSize,
                             "The minimum size to which all time steps are limited to [s]");
        EWOMS_REGISTER_PARAM(TypeTag, unsigned, MaxTimeStepDivisions,
                             "The maximum number of divisions by two of the timestep size "
                             "before the simulation bails out");
    }

    void finishInit()
    {
        linearizeTime_ = 0.0;
        solveTime_ = 0.0;
        updateTime_ = 0.0;
    }

    Scalar solveTime() const
    { return solveTime_; }

    Scalar updateTime() const
    { return updateTime_; }

    template <class Context>
    void boundary(BoundaryRateVector &values,
                  const Context &context,
                  int spaceIdx, int timeIdx) const
    { OPM_THROW(std::logic_error, "Problem does not provide a boundary() method"); }

    template <class Context>
    void constraints(Constraints &constraints,
                     const Context &context,
                     int spaceIdx, int timeIdx) const
    { OPM_THROW(std::logic_error, "Problem does not provide a constraints() method"); }

    template <class Context>
    void source(RateVector &rate,
                const Context &context,
                int spaceIdx, int timeIdx) const
    { OPM_THROW(std::logic_error, "Problem does not provide a source() method"); }

    template <class Context>
    void initial(PrimaryVariables &values,
                 const Context &context,
                 int spaceIdx, int timeIdx) const
    { OPM_THROW(std::logic_error, "Problem does not provide a initial() method"); }

    template <class Context>
    Scalar extrusionFactor(const Context &context,
                           int spaceIdx, int timeIdx) 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& timer = simulator().timer();
        Scalar simulationTime = timer.realTimeElapsed();
        Scalar setupTime = simulator().setupTime();
        Scalar prePostProcessTime = simulator().prePostProcessTime();
        Scalar localCpuTime = timer.cpuTimeElapsed();
        Scalar globalCpuTime = timer.globalCpuTimeElapsed();
        Scalar totalWriteTime = simulator().totalWriteTime();
        int numProcesses = this->gridView().comm().size();
        int threadsPerProcess = ThreadManager::maxThreads();
        if (gridView().comm().rank() == 0) {
            std::cout << std::setprecision(3)
                      << "Simulation of problem '" << asImp_().name() << "' finished.\n"
                      << "\n"
                      << "------------------------ Timing receipt ------------------------\n"
                      << "Setup time: "<< Simulator::humanReadableTime(setupTime, /*isAmendment=*/false)
                      << ", " << setupTime/(simulationTime + setupTime)*100 << "%\n"
                      << "Simulation time: "<< Simulator::humanReadableTime(simulationTime, /*isAmendment=*/false)
                      << ", " << simulationTime/(simulationTime + setupTime)*100 << "%\n"
                      << "    Linearization time: "<< Simulator::humanReadableTime(linearizeTime_, /*isAmendment=*/false)
                      << ", " << linearizeTime_/simulationTime*100 << "%\n"
                      << "    Linear solve time: " << Simulator::humanReadableTime(solveTime_, /*isAmendment=*/false)
                      << ", " << solveTime_/simulationTime*100 << "%\n"
                      << "    Newton update time: " << Simulator::humanReadableTime(updateTime_, /*isAmendment=*/false)
                      << ", " << updateTime_/simulationTime*100 << "%\n"
                      << "    Pre/postprocess time: " << Simulator::humanReadableTime(prePostProcessTime, /*isAmendment=*/false)
                      << ", " << prePostProcessTime/simulationTime*100 << "%\n"
                      << "    Output write time: " << Simulator::humanReadableTime(totalWriteTime, /*isAmendment=*/false)
                      << ", " << totalWriteTime/simulationTime*100 << "%\n"
                      << "First process' simulation CPU time: " <<  Simulator::humanReadableTime(localCpuTime, /*isAmendment=*/false) << "\n"
                      << "Number of processes: " << numProcesses << "\n"
                      << "Threads per processes: " << threadsPerProcess << "\n"
                      << "Total CPU time: " << Simulator::humanReadableTime(globalCpuTime, /*isAmendment=*/false) << "\n"
                      << "\n"
                      << "Note 1: If not stated otherwise, all times are wall clock times\n"
                      << "Note 2: Taxes and administrative overhead are "
                      << (simulationTime - (linearizeTime_+solveTime_+updateTime_+prePostProcessTime+totalWriteTime))/simulationTime*100
                      << "%\n"
                      << "\n"
                      << "Our simulation hours are 24/7. Thank you for choosing us.\n"
                      << "----------------------------------------------------------------\n"
                      << "\n"
                      << std::flush;
        }
    }

    void timeIntegration()
    {
        int maxFails = EWOMS_GET_PARAM(TypeTag, unsigned, MaxTimeStepDivisions);
        Scalar minTimeStepSize = EWOMS_GET_PARAM(TypeTag, Scalar, MinTimeStepSize);

        // if the time step size of the simulator is smaller than
        // the specified minimum size and we're not going to finish
        // the simulation or an episode, try with the minimum size.
        if (simulator().timeStepSize() < minTimeStepSize &&
            !simulator().episodeWillBeOver() &&
            !simulator().willBeFinished())
        {
            simulator().setTimeStepSize(minTimeStepSize);
        }

        for (int i = 0; i < maxFails; ++i) {
            bool converged = model().update(newtonMethod());

            linearizeTime_ += newtonMethod().linearizeTime();
            solveTime_ += newtonMethod().solveTime();
            updateTime_ += newtonMethod().updateTime();

            if (converged)
                return;

            Scalar dt = simulator().timeStepSize();
            Scalar nextDt = dt / 2;
            if (nextDt < minTimeStepSize)
                break; // give up: we can't make the time step smaller anymore!
            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;
        }

        OPM_THROW(std::runtime_error,
                   "Newton solver didn't converge after "
                   << maxFails << " time-step divisions. dt="
                   << simulator().timeStepSize());
    }

    Scalar nextTimeStepSize()
    {
        Scalar dtNext = std::min(EWOMS_GET_PARAM(TypeTag, Scalar, MaxTimeStepSize),
                                 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(); }
    // \}

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

        if (enableVtkOutput_())
            defaultVtkWriter_->beginWrite(t);

        model().prepareOutputFields();

        if (enableVtkOutput_()) {
            model().appendOutputFields(*defaultVtkWriter_);
            defaultVtkWriter_->endWrite();
        }
    }

    VtkMultiWriter &defaultVtkWriter() const
    { return defaultVtkWriter_; }

private:
    bool enableVtkOutput_() const
    { return EWOMS_GET_PARAM(TypeTag, bool, 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_;

    // CPU time keeping
    Scalar linearizeTime_;
    Scalar solveTime_;
    Scalar updateTime_;
};

} // namespace Ewoms

#endif