impetproblem.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_IMPETPROBLEM_HH
#define DUMUX_IMPETPROBLEM_HH
 
#include <dune/common/float_cmp.hh>
#include "impetproperties.hh"
#include <dumux/io/vtkmultiwriter.hh>
#include <dumux/io/restart.hh>
#include <dumux/io/adaptivegridrestart.hh>
 
#include <dumux/porousmediumflow/sequential/gridadapt.hh>
 
namespace Dumux
{
template<class TypeTag>
class IMPETProblem
{
private:
    using Implementation = GetPropType<TypeTag, Properties::Problem>;
    using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
    using Grid = GetPropType<TypeTag, Properties::Grid>;
 
    using TimeManager = GetPropType<TypeTag, Properties::TimeManager>;
 
    using VtkMultiWriter = Dumux::VtkMultiWriter<GridView>;
 
    using Variables = GetPropType<TypeTag, Properties::Variables>;
 
    using SolutionTypes = GetProp<TypeTag, Properties::SolutionTypes>;
    using VertexMapper = typename SolutionTypes::VertexMapper;
    using ElementMapper = typename SolutionTypes::ElementMapper;
 
    using IMPETModel = GetPropType<TypeTag, Properties::Model>;
    using TransportSolutionType = GetPropType<TypeTag, Properties::TransportSolutionType>;
    using PressureModel = GetPropType<TypeTag, Properties::PressureModel>;
    using TransportModel = GetPropType<TypeTag, Properties::TransportModel>;
 
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
 
    enum
        {
            dim = GridView::dimension,
            dimWorld = GridView::dimensionworld
        };
    enum
        {
            wetting = 0, nonwetting = 1,
            adaptiveGrid = getPropValue<TypeTag, Properties::AdaptiveGrid>()
        };
 
    using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
    using Element = typename GridView::Traits::template Codim<0>::Entity;
    using Intersection = typename GridView::Intersection;
    // The module to adapt grid. If adaptiveGrid is false, this model does nothing.
    using GridAdaptModel = GetPropType<TypeTag, Properties::GridAdaptModel>;
 
    using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
    using BoundaryTypes = GetPropType<TypeTag, Properties::SequentialBoundaryTypes>;
 
    //private!! copy constructor
    IMPETProblem(const IMPETProblem&)
    {}
 
public:
    IMPETProblem(TimeManager &timeManager, Grid& grid)
        : IMPETProblem(timeManager, grid, grid.leafGridView())
    {}
 
    IMPETProblem(TimeManager &timeManager, Grid& grid, const GridView& gridView)
        : grid_(&grid),
          gridView_(gridView),
          bBoxMin_(std::numeric_limits<double>::max()),
          bBoxMax_(-std::numeric_limits<double>::max()),
          timeManager_(&timeManager),
          outputInterval_(1),
          outputTimeInterval_(0.0),
          vtkOutputLevel_(-1)
    {
        // calculate the bounding box of the grid view
        using std::min;
        using std::max;
        for (const auto& vertex : vertices(this->gridView())) {
            for (int i=0; i<dim; i++) {
                bBoxMin_[i] = min(bBoxMin_[i], vertex.geometry().center()[i]);
                bBoxMax_[i] = max(bBoxMax_[i], vertex.geometry().center()[i]);
            }
        }
 
        // communicate to get the bounding box of the whole domain
        if (this->gridView().comm().size() > 1)
            for (int i = 0; i < dim; ++i) {
                bBoxMin_[i] = this->gridView().comm().min(bBoxMin_[i]);
                bBoxMax_[i] = this->gridView().comm().max(bBoxMax_[i]);
            }
 
        variables_ = std::make_shared<Variables>(this->gridView());
        pressModel_ = std::make_shared<PressureModel>(asImp_());
 
        transportModel_ = std::make_shared<TransportModel>(asImp_());
        model_ = std::make_shared<IMPETModel>(asImp_()) ;
 
        // create an Object to handle adaptive grids
        if (adaptiveGrid)
            gridAdapt_ = std::make_shared<GridAdaptModel>(asImp_());
 
        vtkOutputLevel_ = getParam<int>("Vtk.OutputLevel");
        dtVariationRestrictionFactor_ = getParam<Scalar>("Impet.DtVariationRestrictionFactor", std::numeric_limits<Scalar>::max());
        maxTimeStepSize_ = getParam<Scalar>("TimeManager.MaxTimeStepSize", std::numeric_limits<Scalar>::max());
    }
 
    void boundaryTypes(BoundaryTypes &bcTypes,
                       const Intersection &intersection) const
    {
        // forward it to the method which only takes the global coordinate
        asImp_().boundaryTypesAtPos(bcTypes, intersection.geometry().center());
    }
 
    void boundaryTypesAtPos(BoundaryTypes &bcTypes,
                            const GlobalPosition &globalPos) const
    {
        // Throw an exception (there is no reasonable default value
        // for Dirichlet conditions)
        DUNE_THROW(Dune::InvalidStateException,
                   "The problem does not provide "
                   "a boundaryTypesAtPos() method.");
    }
 
    void dirichlet(PrimaryVariables &values,
                   const Intersection &intersection) const
    {
        // forward it to the method which only takes the global coordinate
        asImp_().dirichletAtPos(values, intersection.geometry().center());
    }
 
    void dirichletAtPos(PrimaryVariables &values,
                        const GlobalPosition &globalPos) const
    {
        // Throw an exception (there is no reasonable default value
        // for Dirichlet conditions)
        DUNE_THROW(Dune::InvalidStateException,
                   "The problem specifies that some boundary "
                   "segments are dirichlet, but does not provide "
                   "a dirichletAtPos() method.");
    }
 
    void neumann(PrimaryVariables &values,
                 const Intersection &intersection) const
    {
        // forward it to the interface with only the global position
        asImp_().neumannAtPos(values, intersection.geometry().center());
    }
 
    void neumannAtPos(PrimaryVariables &values,
                      const GlobalPosition &globalPos) const
    {
        // Throw an exception (there is no reasonable default value
        // for Neumann conditions)
        DUNE_THROW(Dune::InvalidStateException,
                   "The problem specifies that some boundary "
                   "segments are neumann, but does not provide "
                   "a neumannAtPos() method.");
    }
 
    void source(PrimaryVariables &values,
                const Element &element) const
    {
        // forward to generic interface
        asImp_().sourceAtPos(values, element.geometry().center());
    }
 
    void sourceAtPos(PrimaryVariables &values,
                     const GlobalPosition &globalPos) const
    {         // Throw an exception (there is no initial condition)
        DUNE_THROW(Dune::InvalidStateException,
                   "The problem does not provide "
                   "a sourceAtPos() method.");
    }
 
    void initial(PrimaryVariables &values,
                 const Element &element) const
    {
        // forward to generic interface
        asImp_().initialAtPos(values, element.geometry().center());
    }
 
    void initialAtPos(PrimaryVariables &values,
                      const GlobalPosition &globalPos) const
    {
        // initialize with 0 by default
        values = 0;
    }
 
    void init()
    {
        // set the initial condition of the model
        variables_->initialize();
        model().initialize();
        if (adaptiveGrid)
            gridAdapt().init();
    }
 
    void preTimeStep()
    {
        // if adaptivity is used, this method adapts the grid.
        // if it is not used, this method does nothing.
        if (adaptiveGrid && timeManager().timeStepIndex() > 0)
            this->gridAdapt().adaptGrid();
        asImp_().pressureModel().updateMaterialLaws();
    }
 
    void timeIntegration()
    {
        // allocate temporary vectors for the updates
        TransportSolutionType updateVector;//(this->transportModel().solutionSize());
 
        Scalar t = timeManager().time();
        Scalar dt = std::numeric_limits<Scalar>::max();
 
        // obtain the first update and the time step size
        model().update(t, dt, updateVector);
 
        using std::max;
        using std::min;
        if (t > 0 && timeManager().timeStepSize()> 0.)
        {
            Scalar oldDt = timeManager().timeStepSize();
            Scalar minDt = max(oldDt - dtVariationRestrictionFactor_ * oldDt, 0.0);
            Scalar maxDt = oldDt + dtVariationRestrictionFactor_ * oldDt;
            dt = max(min(dt, maxDt), minDt);
        }
 
        //make sure t_old + dt is not larger than tend
        dt = min(dt, timeManager().episodeMaxTimeStepSize());
 
        // check if we are in first TS and an initialDt was assigned
        if (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(t, 0.0, 1.0e-30)
            && Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(timeManager().timeStepSize(), 0.0, 1.0e-30))
        {
            if (gridView().comm().size() > 1)
                dt = gridView().comm().min(dt);
 
            // check if assigned initialDt is in accordance with dt from first transport step
            if (timeManager().timeStepSize() > dt
                && gridView().comm().rank() == 0)
                Dune::dwarn << "Initial timestep of size " << timeManager().timeStepSize()
                            << " is larger then dt=" << dt <<" from transport" << std::endl;
            // internally assign next timestep size
            dt = min(dt, timeManager().timeStepSize());
        }
 
        //make sure the right time-step is used by all processes in the parallel case
        if (gridView().comm().size() > 1)
            dt = gridView().comm().min(dt);
 
        // check maximum allowed time step size
        timeManager().setTimeStepSize(dt);
 
        // explicit Euler: Sat <- Sat + dt*N(Sat)
        transportModel().updateTransportedQuantity(updateVector);
    }
 
    void postTimeStep()
    {}
 
    void advanceTimeLevel()
    {}
 
    Scalar timeStepSize() const
    { return timeManager().timeStepSize(); }
 
    void setTimeStepSize(const Scalar dt)
    { timeManager().setTimeStepSize(dt); }
 
    Scalar nextTimeStepSize(const Scalar dt)
    { return timeManager().timeStepSize();}
 
    Scalar maxTimeStepSize() const
    { return maxTimeStepSize_; }
 
    bool shouldWriteRestartFile() const
    {
        if (outputInterval_ > 0)
        {
            return
                timeManager().timeStepIndex() > 0 &&
                (timeManager().timeStepIndex() % int(100*outputInterval_) == 0);
        }
        else
            return
                shouldWriteOutput();
    }
 
    void setOutputTimeInterval(const Scalar timeInterval)
    {
        outputTimeInterval_ = (timeInterval > 0.0) ? timeInterval : 1e100;
        timeManager().startNextEpisode(outputTimeInterval_);
    }
 
    void setOutputInterval(const int interval)
    {
        using std::max;
        outputInterval_ = max(interval, 0);
    }
 
    bool shouldWriteOutput() const
    {
        if (outputInterval_ > 0)
        {
            if (timeManager().timeStepIndex() % outputInterval_ == 0
                || timeManager().willBeFinished()
                || timeManager().episodeWillBeFinished())
            {
                return true;
            }
        }
        else if (timeManager().willBeFinished()
                 || timeManager().episodeWillBeFinished() || timeManager().timeStepIndex() == 0)
        {
            return true;
        }
        return false;
    }
 
    void episodeEnd()
    {
        if (outputTimeInterval_ > 0.0 && !timeManager().finished())
        {
            timeManager().startNextEpisode(outputTimeInterval_);
        }
        else if (!timeManager().finished())
        {
            std::cerr << "The end of an episode is reached, but the problem "
                      << "does not override the episodeEnd() method. "
                      << "Doing nothing!\n";
        }
    }
 
    // \}
 
    const std::string& name() const
    {
        return simname_;
    }
 
    void setName(std::string newName)
    {
        simname_ = newName;
    }
 
    const GridView& gridView() const
    {
        return gridView_;
    }
 
    Grid &grid()
    {
        if (!grid_)
        {
            DUNE_THROW(Dune::InvalidStateException, "Grid was called in problemclass, "
                       << "although it is not specified. Do so by using setGrid() method!");
        }
        return *grid_;
    }
    void setGrid(Grid &grid)
    {
        grid_ = &grid;
    }
 
    GridAdaptModel& gridAdapt()
    {
        return *gridAdapt_;
    }
 
    void preAdapt()
    {
    }
 
    void postAdapt()
    {
        // write out new grid
        //        Dune::VTKWriter<LeafGridView> vtkwriter(problem_.gridView());
        //        vtkwriter.write("latestgrid",Dune::VTK::appendedraw);
    }
 
    const VertexMapper &vertexMapper() const
    { return variables_->vertexMapper(); }
 
    const ElementMapper &elementMapper() const
    { return variables_->elementMapper(); }
 
    const GlobalPosition &bBoxMin() const
    { return bBoxMin_; }
 
    const GlobalPosition &bBoxMax() const
    { return bBoxMax_; }
 
 
 
    TimeManager &timeManager()
    { return *timeManager_; }
 
    const TimeManager &timeManager() const
    { return *timeManager_; }
 
    Variables& variables ()
    { return *variables_; }
 
    const Variables& variables () const
    { return *variables_; }
 
    IMPETModel &model()
    { return *model_; }
 
    const IMPETModel &model() const
    { return *model_; }
 
    PressureModel &pressureModel()
    { return *pressModel_; }
 
    const PressureModel &pressureModel() const
    { return *pressModel_; }
 
    TransportModel &transportModel()
    { return *transportModel_; }
 
    const TransportModel &transportModel() const
    { return *transportModel_; }
    // \}
 
    // \{
 
    void serialize()
    {
        using Restarter = Restart;
 
        Restarter res;
        res.serializeBegin(asImp_());
        std::cout << "Serialize to file " << res.fileName() << "\n";
 
        timeManager().serialize(res);
        resultWriter().serialize(res);
 
        // do the actual serialization process: write primary variables
        res.template serializeEntities<0> (*pressModel_, gridView());
        res.template serializeEntities<0> (*transportModel_, gridView());
 
        res.serializeEnd();
 
        if (adaptiveGrid)
        {
            AdaptiveGridRestart<Grid>::serializeGrid(asImp_());
        }
    }
 
    void restart(const double tRestart)
    {
        if (adaptiveGrid)
        {
            AdaptiveGridRestart<Grid>::restartGrid(asImp_());
            variables().initialize();
            model().initialize();
        }
 
        using Restarter = Restart;
 
        Restarter res;
        res.deserializeBegin(asImp_(), tRestart);
        std::cout << "Deserialize from file " << res.fileName() << "\n";
 
        timeManager().deserialize(res);
 
        resultWriter().deserialize(res);
 
        // do the actual serialization process: get primary variables
        res.template deserializeEntities<0> (*pressModel_, gridView());
        res.template deserializeEntities<0> (*transportModel_, gridView());
        pressureModel().updateMaterialLaws();
 
        res.deserializeEnd();
    }
 
    void addOutputVtkFields()
    {
    }
    int vtkOutputLevel() const
    {
        return vtkOutputLevel_;
    }
 
    void writeOutput(bool verbose = true)
    {
        if (verbose && gridView().comm().rank() == 0)
            std::cout << "Writing result file for current time step\n";
 
        if (!resultWriter_)
            resultWriter_ = std::make_shared<VtkMultiWriter>(gridView(), asImp_().name());
        if (adaptiveGrid)
            resultWriter_->gridChanged();
        resultWriter_->beginWrite(timeManager().time() + timeManager().timeStepSize());
        model().addOutputVtkFields(*resultWriter_);
        asImp_().addOutputVtkFields();
        resultWriter_->endWrite();
    }
 
    // \}
 
protected:
    VtkMultiWriter& resultWriter()
    {
        if (!resultWriter_)
            resultWriter_ = std::make_shared<VtkMultiWriter>(gridView(), asImp_().name());
        return *resultWriter_;
    }
    VtkMultiWriter& resultWriter() const
    {
        if (!resultWriter_)
            resultWriter_ = std::make_shared<VtkMultiWriter>(gridView(), asImp_().name());
        return *resultWriter_;
    }
 
private:
    Implementation &asImp_()
    { return *static_cast<Implementation *>(this); }
 
    const Implementation &asImp_() const
    { return *static_cast<const Implementation *>(this); }
 
    std::string simname_; // a string for the name of the current simulation,
    // which could be set by means of a program argument, for example.
 
    // pointer to a possibly adaptive grid.
    Grid *grid_;
    GridView gridView_;
 
    GlobalPosition bBoxMin_;
    GlobalPosition bBoxMax_;
 
    TimeManager *timeManager_;
    Scalar maxTimeStepSize_;
 
    std::shared_ptr<Variables> variables_;
 
    std::shared_ptr<PressureModel> pressModel_;
    std::shared_ptr<TransportModel> transportModel_;
    std::shared_ptr<IMPETModel> model_;
 
    std::shared_ptr<VtkMultiWriter> resultWriter_;
    int outputInterval_;
    Scalar outputTimeInterval_;
    int vtkOutputLevel_;
    std::shared_ptr<GridAdaptModel> gridAdapt_;
 
    Scalar dtVariationRestrictionFactor_;
};
}
#endif