onemodelproblem.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_ONE_MODEL_PROBLEM_HH
#define DUMUX_ONE_MODEL_PROBLEM_HH
 
#include <dune/common/shared_ptr.hh>
#include <dumux/common/properties.hh>
#include <dumux/porousmediumflow/sequential/properties.hh>
#include <dumux/io/vtkmultiwriter.hh>
#include <dumux/io/restart.hh>
 
 
namespace Dumux
{
 
template<class TypeTag>
class OneModelProblem
{
private:
    using Implementation = GetPropType<TypeTag, Properties::Problem>;
    using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
    using Grid = typename GridView::Grid;
 
    using TimeManager = GetPropType<TypeTag, Properties::TimeManager>;
 
    using VtkMultiWriter = Dumux::VtkMultiWriter<GridView>;
 
    using Variables = GetPropType<TypeTag, Properties::Variables>;
 
    using Model = GetPropType<TypeTag, Properties::Model>;
 
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
 
    using SolutionTypes = GetProp<TypeTag, Properties::SolutionTypes>;
    using VertexMapper = typename SolutionTypes::VertexMapper;
    using ElementMapper = typename SolutionTypes::ElementMapper;
 
    enum
    {
        dim = GridView::dimension,
        dimWorld = GridView::dimensionworld
    };
    enum
    {
        wetting = 0, nonwetting = 1
    };
 
    using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
    using Element = typename GridView::Traits::template Codim<0>::Entity;
    using Intersection = typename GridView::Intersection;
 
    using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
    using BoundaryTypes = GetPropType<TypeTag, Properties::SequentialBoundaryTypes>;
 
    // private!! copy constructor
    OneModelProblem(const OneModelProblem&)
    {}
 
public:
 
 
    OneModelProblem(Grid& grid, bool verbose = true)
        : gridView_(grid.leafGridView()),
          bBoxMin_(std::numeric_limits<double>::max()),
          bBoxMax_(-std::numeric_limits<double>::max()),
          variables_(grid.leafGridView()),
          outputInterval_(1),
          outputTimeInterval_(0)
    {
        // calculate the bounding box of the grid view
        using std::max;
        using std::min;
        for (const auto& vertex : vertices(grid.leafGridView())) {
            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]);
            }
        }
 
        timeManager_ = std::make_shared<TimeManager>(verbose);
 
        model_ = std::make_shared<Model>(asImp_()) ;
        maxTimeStepSize_ = getParam<Scalar>("TimeManager.MaxTimeStepSize", std::numeric_limits<Scalar>::max());
    }
 
 
    OneModelProblem(TimeManager& timeManager, Grid& grid)
        : gridView_(grid.leafGridView()),
          bBoxMin_(std::numeric_limits<double>::max()),
          bBoxMax_(-std::numeric_limits<double>::max()),
          variables_(grid.leafGridView()),
          outputInterval_(1),
          outputTimeInterval_(0)
    {
        // calculate the bounding box of the grid view
        using std::max;
        using std::min;
        for (const auto& vertex : vertices(grid.leafGridView())) {
            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]);
            }
        }
 
        timeManager_ = Dune::stackobject_to_shared_ptr<TimeManager>(timeManager);
 
        model_ = std::make_shared<Model>(asImp_()) ;
        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();
    }
 
    void preTimeStep()
    {}
 
    void timeIntegration()
    {}
 
    void postTimeStep()
    {}
 
    void advanceTimeLevel()
    {}
 
    Scalar maxTimeStepSize() const
    { return maxTimeStepSize_; }
 
    Scalar timeStepSize() const
    { return timeManager().timeStepSize(); }
 
    void setTimeStepSize(Scalar dt)
    { timeManager().setTimeStepSize(dt); }
 
    Scalar nextTimeStepSize(Scalar dt)
    { return timeManager().timeStepSize();}
 
    bool shouldWriteRestartFile() const
    {
        return
            timeManager().timeStepIndex() > 0 &&
            (timeManager().timeStepIndex() % 5 == 0);
    }
 
    void setOutputTimeInterval(const Scalar timeInterval)
    {
        outputTimeInterval_ = (timeInterval > 0.0) ? timeInterval : 1e100;
        timeManager().startNextEpisode(outputTimeInterval_);
    }
 
    void setOutputInterval(int interval)
    { outputInterval_ = interval; }
 
    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 addOutputVtkFields()
    {}
 
    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());
        resultWriter_->beginWrite(timeManager().time() + timeManager().timeStepSize());
        model().addOutputVtkFields(*resultWriter_);
        asImp_().addOutputVtkFields();
        resultWriter_->endWrite();
    }
 
    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(const std::string& newName)
    {
        simname_ = newName;
    }
 
    const GridView &gridView() const
    { return gridView_; }
 
    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_; }
 
    Model &model()
    { return *model_; }
 
    const Model &model() const
    { return *model_; }
    // \}
 
 
    // \{
 
    void serialize()
    {
        using Restarter = Restart;
 
        Restarter res;
        res.serializeBegin(asImp_());
        std::cout << "Serialize to file " << res.fileName() << "\n";
 
        timeManager().serialize(res);
        resultWriter().serialize(res);
        res.template deserializeEntities<0> (model(), gridView_);
 
        res.serializeEnd();
    }
 
    void restart(double tRestart)
    {
        using Restarter = Restart;
 
        Restarter res;
        res.deserializeBegin(asImp_(), tRestart);
        std::cout << "Deserialize from file " << res.fileName() << "\n";
 
        timeManager().deserialize(res);
        resultWriter().deserialize(res);
        res.template deserializeEntities<0> (model(), gridView_);
 
        res.deserializeEnd();
    }
 
    // \}
 
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 an program argument,
                                 // for example.
    const GridView gridView_;
 
    GlobalPosition bBoxMin_;
    GlobalPosition bBoxMax_;
 
    std::shared_ptr<TimeManager> timeManager_;
    Scalar maxTimeStepSize_;
 
    Variables variables_;
 
    std::shared_ptr<Model> model_;
 
    std::shared_ptr<VtkMultiWriter> resultWriter_;
    int outputInterval_;
    Scalar outputTimeInterval_;
};
 
}
#endif