common/fvproblem.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_COMMON_FV_PROBLEM_HH
#define DUMUX_COMMON_FV_PROBLEM_HH
 
#include <memory>
#include <map>
 
#include <dune/common/fvector.hh>
#include <dune/grid/common/gridenums.hh>
 
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/discretization/method.hh>
 
namespace Dumux {
 
template<class TypeTag>
class FVProblem
{
    using Implementation = GetPropType<TypeTag, Properties::Problem>;
 
    using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
    using FVElementGeometry = typename GridGeometry::LocalView;
    using GridView = typename GridGeometry::GridView;
    using SubControlVolume = typename FVElementGeometry::SubControlVolume;
    using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
    using Element = typename GridView::template Codim<0>::Entity;
    using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
 
    enum { dim = GridView::dimension };
 
    using PointSource = GetPropType<TypeTag, Properties::PointSource>;
    using PointSourceHelper = GetPropType<TypeTag, Properties::PointSourceHelper>;
    using PointSourceMap = std::map< std::pair<std::size_t, std::size_t>,
                                     std::vector<PointSource> >;
 
    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
    using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView;
    using ElementVolumeVariables = typename GridVariables::GridVolumeVariables::LocalView;
    using VolumeVariables = typename ElementVolumeVariables::VolumeVariables;
 
    using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
 
    static constexpr bool isBox = GridGeometry::discMethod == DiscretizationMethod::box;
    static constexpr bool isStaggered = GridGeometry::discMethod == DiscretizationMethod::staggered;
 
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
    using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
 
public:
    struct Traits
    {
        using Scalar = FVProblem::Scalar;
        using PrimaryVariables = FVProblem::PrimaryVariables;
        using NumEqVector = FVProblem::NumEqVector;
        using BoundaryTypes = FVProblem::BoundaryTypes;
    };
 
    FVProblem(std::shared_ptr<const GridGeometry> gridGeometry, const std::string& paramGroup = "")
    : gridGeometry_(gridGeometry)
    , paramGroup_(paramGroup)
    {
        // set a default name for the problem
        problemName_ = getParamFromGroup<std::string>(paramGroup, "Problem.Name");
    }
 
    const std::string& name() const
    {
        return problemName_;
    }
 
    void setName(const std::string& newName)
    {
        problemName_ = newName;
    }
 
    // \{
 
    BoundaryTypes boundaryTypes(const Element &element,
                                const SubControlVolume &scv) const
    {
        if (!isBox)
            DUNE_THROW(Dune::InvalidStateException,
                       "boundaryTypes(..., scv) called for cell-centered method.");
 
        // forward it to the method which only takes the global coordinate
        return asImp_().boundaryTypesAtPos(scv.dofPosition());
    }
 
    BoundaryTypes boundaryTypes(const Element &element,
                                const SubControlVolumeFace &scvf) const
    {
        if (isBox)
            DUNE_THROW(Dune::InvalidStateException,
                       "boundaryTypes(..., scvf) called for box method.");
 
        // forward it to the method which only takes the global coordinate
        return asImp_().boundaryTypesAtPos(scvf.ipGlobal());
    }
 
    BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
    {
        BoundaryTypes bcTypes;
        bcTypes.setAllDirichlet();
        return bcTypes;
    }
 
    PrimaryVariables dirichlet(const Element &element, const SubControlVolumeFace &scvf) const
    {
        // forward it to the method which only takes the global coordinate
        if (isBox)
        {
            DUNE_THROW(Dune::InvalidStateException, "dirichlet(scvf) called for box method.");
        }
        else
            return asImp_().dirichletAtPos(scvf.ipGlobal());
    }
 
    PrimaryVariables dirichlet(const Element &element, const SubControlVolume &scv) const
    {
        // forward it to the method which only takes the global coordinate
        if (!isBox && !isStaggered)
        {
            DUNE_THROW(Dune::InvalidStateException, "dirichlet(scv) called for other than box or staggered method.");
        }
        else
            return asImp_().dirichletAtPos(scv.dofPosition());
    }
 
    PrimaryVariables dirichletAtPos(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 dirichlet() or a dirichletAtPos() method.");
    }
 
    static constexpr bool enableInternalDirichletConstraints()
    { return false; }
 
    NumEqVector neumann(const Element& element,
                        const FVElementGeometry& fvGeometry,
                        const ElementVolumeVariables& elemVolVars,
                        const ElementFluxVariablesCache& elemFluxVarsCache,
                        const SubControlVolumeFace& scvf) const
    {
        // forward it to the interface with only the global position
        return asImp_().neumannAtPos(scvf.ipGlobal());
    }
 
    NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
    {
        return NumEqVector(0.0);
    }
 
    NumEqVector source(const Element &element,
                       const FVElementGeometry& fvGeometry,
                       const ElementVolumeVariables& elemVolVars,
                       const SubControlVolume &scv) const
    {
        // forward to solution independent, fully-implicit specific interface
        return asImp_().sourceAtPos(scv.center());
    }
 
    NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
    {
        return NumEqVector(0.0);
    }
 
    void addPointSources(std::vector<PointSource>& pointSources) const {}
 
    void pointSource(PointSource& source,
                     const Element &element,
                     const FVElementGeometry& fvGeometry,
                     const ElementVolumeVariables& elemVolVars,
                     const SubControlVolume &scv) const
    {
        // forward to space dependent interface method
        asImp_().pointSourceAtPos(source, source.position());
    }
 
    void pointSourceAtPos(PointSource& pointSource,
                          const GlobalPosition &globalPos) const {}
 
    template<class MatrixBlock>
    void addSourceDerivatives(MatrixBlock& block,
                              const Element& element,
                              const FVElementGeometry& fvGeometry,
                              const VolumeVariables& volVars,
                              const SubControlVolume& scv) const {}
 
    NumEqVector scvPointSources(const Element &element,
                                const FVElementGeometry& fvGeometry,
                                const ElementVolumeVariables& elemVolVars,
                                const SubControlVolume &scv) const
    {
        NumEqVector source(0);
        auto scvIdx = scv.indexInElement();
        auto key = std::make_pair(gridGeometry_->elementMapper().index(element), scvIdx);
        if (pointSourceMap_.count(key))
        {
            // Add the contributions to the dof source values
            // We divide by the volume. In the local residual this will be multiplied with the same
            // factor again. That's because the user specifies absolute values in kg/s.
            const auto volume = scv.volume()*elemVolVars[scv].extrusionFactor();
 
            for (const auto& ps : pointSourceMap_.at(key))
            {
                // we make a copy of the local point source here
                auto pointSource = ps;
 
                // Note: two concepts are implemented here. The PointSource property can be set to a
                // customized point source function achieving variable point sources,
                // see TimeDependentPointSource for an example. The second imitated the standard
                // dumux source interface with solDependentPointSource / pointSourceAtPos, methods
                // that can be overloaded in the actual problem class also achieving variable point sources.
                // The first one is more convenient for simple function like a time dependent source.
                // The second one might be more convenient for e.g. a solution dependent point source.
 
                // we do an update e.g. used for TimeDependentPointSource
                pointSource.update(asImp_(), element, fvGeometry, elemVolVars, scv);
                // call convienience problem interface function
                asImp_().pointSource(pointSource, element, fvGeometry, elemVolVars, scv);
                // at last take care about multiplying with the correct volume
                pointSource /= volume*pointSource.embeddings();
                // add the point source values to the local residual
                source += pointSource.values();
            }
        }
 
        return source;
    }
 
    void computePointSourceMap()
    {
        // clear the given point source maps in case it's not empty
        pointSourceMap_.clear();
 
        // get and apply point sources if any given in the problem
        std::vector<PointSource> sources;
        asImp_().addPointSources(sources);
 
        // if there are point sources compute the DOF to point source map
        if (!sources.empty())
        {
            // calculate point source locations and save them in a map
            PointSourceHelper::computePointSourceMap(*gridGeometry_,
                                                     sources,
                                                     pointSourceMap_);
        }
    }
 
    const PointSourceMap& pointSourceMap() const
    { return pointSourceMap_; }
 
    void applyInitialSolution(SolutionVector& sol) const
    {
        // set the initial values by forwarding to a specialized method
        applyInitialSolutionImpl_(sol, std::integral_constant<bool, isBox>());
    }
 
    template<class Entity>
    PrimaryVariables initial(const Entity& entity) const
    {
        static_assert(int(Entity::codimension) == 0 || int(Entity::codimension) == dim, "Entity must be element or vertex");
        return asImp_().initialAtPos(entity.geometry().center());
    }
 
    PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
    {
        // Throw an exception (there is no reasonable default value
        // for initial values)
        DUNE_THROW(Dune::InvalidStateException,
                   "The problem does not provide "
                   "an initial() or an initialAtPos() method.");
    }
 
    template<class ElementSolution>
    Scalar extrusionFactor(const Element& element,
                           const SubControlVolume& scv,
                           const ElementSolution& elemSol) const
    {
        // forward to generic interface
        return asImp_().extrusionFactorAtPos(scv.center());
    }
 
    Scalar extrusionFactorAtPos(const GlobalPosition &globalPos) const
    {
        // As a default, i.e. if the user's problem does not overload
        // any extrusion factor method, return 1.0
        return 1.0;
    }
 
    // \}
 
    const GridGeometry& gridGeometry() const
    { return *gridGeometry_; }
 
    const std::string& paramGroup() const
    { return paramGroup_; }
 
protected:
    Implementation &asImp_()
    { return *static_cast<Implementation *>(this); }
 
    const Implementation &asImp_() const
    { return *static_cast<const Implementation *>(this); }
 
private:
    void applyInitialSolutionImpl_(SolutionVector& sol, /*isBox=*/std::true_type) const
    {
        const auto numDofs = gridGeometry_->vertexMapper().size();
        const auto numVert = gridGeometry_->gridView().size(dim);
        sol.resize(numDofs);
 
        // if there are more dofs than vertices (enriched nodal dofs), we have to
        // call initial for all dofs at the nodes, coming from all neighboring elements.
        if (numDofs != numVert)
        {
            std::vector<bool> dofVisited(numDofs, false);
            for (const auto& element : elements(gridGeometry_->gridView()))
            {
                for (int i = 0; i < element.subEntities(dim); ++i)
                {
                    const auto dofIdxGlobal = gridGeometry_->vertexMapper().subIndex(element, i, dim);
 
                    // forward to implementation if value at dof is not set yet
                    if (!dofVisited[dofIdxGlobal])
                    {
                        sol[dofIdxGlobal] = asImp_().initial(element.template subEntity<dim>(i));
                        dofVisited[dofIdxGlobal] = true;
                    }
                }
            }
        }
 
        // otherwise we directly loop over the vertices
        else
        {
            for (const auto& vertex : vertices(gridGeometry_->gridView()))
            {
                const auto dofIdxGlobal = gridGeometry_->vertexMapper().index(vertex);
                sol[dofIdxGlobal] = asImp_().initial(vertex);
            }
        }
    }
 
    void applyInitialSolutionImpl_(SolutionVector& sol, /*isBox=*/std::false_type) const
    {
        sol.resize(gridGeometry_->numDofs());
        for (const auto& element : elements(gridGeometry_->gridView()))
        {
            const auto dofIdxGlobal = gridGeometry_->elementMapper().index(element);
            sol[dofIdxGlobal] = asImp_().initial(element);
        }
    }
 
    std::shared_ptr<const GridGeometry> gridGeometry_;
 
    std::string paramGroup_;
 
    std::string problemName_;
 
    PointSourceMap pointSourceMap_;
};
 
} // end namespace Dumux
 
#endif