freeflow/navierstokes/momentum/problem.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_NAVIERSTOKES_MOMENTUM_PROBLEM_HH
#define DUMUX_NAVIERSTOKES_MOMENTUM_PROBLEM_HH
 
#include <dune/common/exceptions.hh>
#include <dune/common/typetraits.hh>
#include <dumux/common/numeqvector.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/fvproblemwithspatialparams.hh>
#include <dumux/discretization/method.hh>
#include <dumux/freeflow/navierstokes/momentum/boundarytypes.hh>
 
namespace Dumux {
 
template<class TypeTag, class DiscretizationMethod>
class NavierStokesMomentumProblemImpl;
 
template<class TypeTag>
class NavierStokesMomentumProblemImpl<TypeTag, DiscretizationMethods::FCStaggered>
: public FVProblemWithSpatialParams<TypeTag>
{
    using ParentType = FVProblemWithSpatialParams<TypeTag>;
    using Implementation = GetPropType<TypeTag, Properties::Problem>;
 
    using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
    using GridView = typename GridGeometry::GridView;
    using Element = typename GridView::template Codim<0>::Entity;
 
    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
    using GridVolumeVariables = typename GridVariables::GridVolumeVariables;
    using ElementVolumeVariables = typename GridVolumeVariables::LocalView;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
 
    using FVElementGeometry = typename GridGeometry::LocalView;
    using SubControlVolume = typename FVElementGeometry::SubControlVolume;
    using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
    using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
    enum {
        dim = GridView::dimension,
        dimWorld = GridView::dimensionworld
      };
 
    using GlobalPosition = typename SubControlVolumeFace::GlobalPosition;
    using VelocityVector = Dune::FieldVector<Scalar, dimWorld>;
    using GravityVector = Dune::FieldVector<Scalar, dimWorld>;
    using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
    using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
 
    static constexpr bool isCoupled_ = !std::is_empty_v<CouplingManager>;
 
 
public:
    using InitialValues = Dune::FieldVector<Scalar, dimWorld>;
    using Sources = Dune::FieldVector<Scalar, dimWorld>;
    using DirichletValues = Dune::FieldVector<Scalar, dimWorld>;
    using BoundaryFluxes = Dune::FieldVector<Scalar, dimWorld>;
 
    using MomentumFluxType = Dune::FieldVector<Scalar, dimWorld>;
 
    using BoundaryTypes = NavierStokesMomentumBoundaryTypes<ModelTraits::dim()>;
 
    static constexpr bool isMomentumProblem() { return true; }
 
    NavierStokesMomentumProblemImpl(std::shared_ptr<const GridGeometry> gridGeometry,
                                    std::shared_ptr<CouplingManager> couplingManager,
                                    const std::string& paramGroup = "")
    : ParentType(gridGeometry, paramGroup)
    , gravity_(0.0)
    , couplingManager_(couplingManager)
    {
        if (getParamFromGroup<bool>(paramGroup, "Problem.EnableGravity"))
            gravity_[dim-1]  = -9.81;
 
        enableInertiaTerms_ = getParamFromGroup<bool>(paramGroup, "Problem.EnableInertiaTerms");
    }
 
    NavierStokesMomentumProblemImpl(std::shared_ptr<const GridGeometry> gridGeometry,
                                    const std::string& paramGroup = "")
    : NavierStokesMomentumProblemImpl(gridGeometry, {}, paramGroup)
    {}
 
    Sources 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.dofPosition());
    }
 
    Sources sourceAtPos(const GlobalPosition& globalPos) const
    {
        return Sources(0.0);
    }
 
    auto boundaryTypes(const Element& element,
                       const SubControlVolumeFace& scvf) const
    {
        // Forward it to the method which only takes the global coordinate.
        // We evaluate the boundary type at the center of the sub control volume face
        // in order to avoid ambiguities at domain corners.
        return asImp_().boundaryTypesAtPos(scvf.center());
    }
 
    DirichletValues dirichlet(const Element& element, const SubControlVolumeFace& scvf) const
    {
        return asImp_().dirichletAtPos(scvf.ipGlobal());
    }
 
    template<class ElementFluxVariablesCache>
    BoundaryFluxes neumann(const Element& element,
                           const FVElementGeometry& fvGeometry,
                           const ElementVolumeVariables& elemVolVars,
                           const ElementFluxVariablesCache& elemFluxVarsCache,
                           const SubControlVolumeFace& scvf) const
    { return asImp_().neumannAtPos(scvf.ipGlobal()); }
 
    BoundaryFluxes neumannAtPos(const GlobalPosition& globalPos) const
    { return BoundaryFluxes(0.0); } 
 
    const GravityVector& gravity() const
    { return gravity_; }
 
    bool enableInertiaTerms() const
    { return enableInertiaTerms_; }
 
    Scalar pressure(const Element& element,
                    const FVElementGeometry& fvGeometry,
                    const SubControlVolumeFace& scvf) const
    {
        if constexpr (isCoupled_)
            return couplingManager_->pressure(element, fvGeometry, scvf);
        else
            return asImp_().pressureAtPos(scvf.ipGlobal());
    }
 
    Scalar pressureAtPos(const GlobalPosition&) const
    {
        DUNE_THROW(Dune::NotImplemented, "pressureAtPos not implemented");
    }
 
    Scalar referencePressure(const Element& element,
                             const FVElementGeometry& fvGeometry,
                             const SubControlVolumeFace& scvf) const
    { return 0.0; }
 
    Scalar density(const Element& element,
                   const FVElementGeometry& fvGeometry,
                   const SubControlVolumeFace& scvf) const
    {
        if constexpr (isCoupled_)
            return couplingManager_->density(element, fvGeometry, scvf);
        else
            return asImp_().densityAtPos(scvf.ipGlobal());
    }
 
    Scalar density(const Element& element,
                   const SubControlVolume& scv,
                   const bool isPreviousTimeStep = false) const
    {
        if constexpr (isCoupled_)
            return couplingManager_->density(element, scv, isPreviousTimeStep);
        else
            return asImp_().densityAtPos(scv.dofPosition());
    }
 
    auto insideAndOutsideDensity(const Element& element,
                                 const FVElementGeometry& fvGeometry,
                                 const SubControlVolumeFace& scvf,
                                 const bool isPreviousTimeStep = false) const
    {
        if constexpr (isCoupled_)
            return couplingManager_->insideAndOutsideDensity(element, fvGeometry, scvf, isPreviousTimeStep);
        else
        {
            const auto rho = asImp_().densityAtPos(scvf.ipGlobal());
            return std::make_pair(rho, rho);
        }
    }
 
    Scalar densityAtPos(const GlobalPosition&) const
    {
        DUNE_THROW(Dune::NotImplemented, "densityAtPos not implemented");
    }
 
    Scalar effectiveViscosity(const Element& element,
                              const FVElementGeometry& fvGeometry,
                              const SubControlVolumeFace& scvf) const
    {
        if constexpr (isCoupled_)
            return couplingManager_->effectiveViscosity(element, fvGeometry, scvf);
        else
            return asImp_().effectiveViscosityAtPos(scvf.ipGlobal());
    }
 
    Scalar effectiveViscosityAtPos(const GlobalPosition&) const
    {
        DUNE_THROW(Dune::NotImplemented, "effectiveViscosityAtPos not implemented");
    }
 
    template<class SolutionVector>
    void applyInitialSolution(SolutionVector& sol) const
    {
        sol.resize(this->gridGeometry().numDofs());
        std::vector<bool> dofHandled(this->gridGeometry().numDofs(), false);
        auto fvGeometry = localView(this->gridGeometry());
        for (const auto& element : elements(this->gridGeometry().gridView()))
        {
            fvGeometry.bindElement(element);
            for (const auto& scv : scvs(fvGeometry))
            {
                const auto dofIdx = scv.dofIndex();
                if (!dofHandled[dofIdx])
                {
                    dofHandled[dofIdx] = true;
                    sol[dofIdx] = asImp_().initial(scv)[scv.dofAxis()];
                }
            }
        }
    }
 
    InitialValues initial(const SubControlVolume& scv) const
    {
        return asImp_().initialAtPos(scv.dofPosition());
    }
 
    Scalar pseudo3DWallFriction(const Element& element,
                                const FVElementGeometry& fvGeometry,
                                const ElementVolumeVariables& elemVolVars,
                                const SubControlVolume& scv,
                                const Scalar height,
                                const Scalar factor = 8.0) const
    {
        const Scalar velocity = elemVolVars[scv].velocity();
        const auto scvf = scvfs(fvGeometry, scv).begin();
        const Scalar viscosity = effectiveViscosity(element, fvGeometry, *scvf);
        return pseudo3DWallFriction(velocity, viscosity, height, factor);
    }
 
    Scalar pseudo3DWallFriction(const Scalar velocity,
                                const Scalar viscosity,
                                const Scalar height,
                                const Scalar factor = 8.0) const
    {
        static_assert(dim == 2, "Pseudo 3D wall friction may only be used in 2D");
        return -factor * velocity * viscosity / (height*height);
    }
 
    bool onSlipBoundary(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) const
    { return asImp_().onSlipBoundaryAtPos(scvf.center()); }
 
    bool onSlipBoundaryAtPos(const GlobalPosition& pos) const
    { return false; }
 
    Scalar permeability(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) const
    {
        DUNE_THROW(Dune::NotImplemented,
            "When using the Beavers-Joseph-Saffman boundary condition, "
            "the permeability must be returned in the actual problem"
        );
    }
 
    Scalar alphaBJ(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) const
    {
        DUNE_THROW(Dune::NotImplemented,
            "When using the Beavers-Joseph-Saffman boundary condition, "
            "the alpha value must be returned in the actual problem"
        );
    }
 
    Scalar betaBJ(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf, const GlobalPosition& tangentialVector) const
    {
        const Scalar interfacePermeability = interfacePermeability_(fvGeometry, scvf, tangentialVector);
        using std::sqrt;
        return asImp_().alphaBJ(fvGeometry, scvf) / sqrt(interfacePermeability);
    }
 
    VelocityVector porousMediumVelocity(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) const
    {
        return VelocityVector(0.0);
    }
 
    const VelocityVector beaversJosephVelocity(const FVElementGeometry& fvGeometry,
                                               const SubControlVolumeFace& scvf,
                                               const ElementVolumeVariables& elemVolVars,
                                               Scalar tangentialVelocityGradient) const
    {
        assert(scvf.isLateral());
        assert(scvf.boundary());
 
        const auto& scv = fvGeometry.scv(scvf.insideScvIdx());
 
        // create a unit normal vector oriented in positive coordinate direction
        GlobalPosition orientation(0.0);
        orientation[scv.dofAxis()] = 1.0;
 
        // du/dy + dv/dx = beta * (u_boundary-uPM)
        // beta = alpha/sqrt(K)
        const Scalar betaBJ = asImp_().betaBJ(fvGeometry, scvf, orientation);
        const Scalar distanceNormalToBoundary = (scvf.ipGlobal() - scv.dofPosition()).two_norm();
 
        static const bool onlyNormalGradient = getParamFromGroup<bool>(this->paramGroup(), "FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph", false);
        if (onlyNormalGradient)
            tangentialVelocityGradient = 0.0;
 
        const Scalar scalarSlipVelocity = (tangentialVelocityGradient*distanceNormalToBoundary
            + asImp_().porousMediumVelocity(fvGeometry, scvf) * orientation * betaBJ * distanceNormalToBoundary
            + elemVolVars[scv].velocity()) / (betaBJ*distanceNormalToBoundary + 1.0);
 
        orientation *= scalarSlipVelocity;
        return orientation;
    }
 
    const CouplingManager& couplingManager() const
    {
        if constexpr (isCoupled_)
            return *couplingManager_;
        else
            DUNE_THROW(Dune::InvalidStateException,
                "Accessing coupling manager of an uncoupled problem is not possible."
            );
    }
 
private:
    template<class Scvf>
    Scalar interfacePermeability_(const FVElementGeometry& fvGeometry, const Scvf& scvf, const GlobalPosition& tangentialVector) const
    {
        const auto& K = asImp_().permeability(fvGeometry, scvf);
 
        // use t*K*t for permeability tensors
        if constexpr (Dune::IsNumber<std::decay_t<decltype(K)>>::value)
            return K;
        else
            return vtmv(tangentialVector, K, tangentialVector);
    }
 
    Implementation& asImp_()
    { return *static_cast<Implementation *>(this); }
 
    const Implementation& asImp_() const
    { return *static_cast<const Implementation *>(this); }
 
    GravityVector gravity_;
    bool enableInertiaTerms_;
    std::shared_ptr<CouplingManager> couplingManager_;
};
 
template<class TypeTag, class DM>
class NavierStokesMomentumProblemImpl<TypeTag, DiscretizationMethods::CVFE<DM>>
: public FVProblemWithSpatialParams<TypeTag>
{
    using ParentType = FVProblemWithSpatialParams<TypeTag>;
    using Implementation = GetPropType<TypeTag, Properties::Problem>;
 
    using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
    using GridView = typename GridGeometry::GridView;
    using Element = typename GridView::template Codim<0>::Entity;
 
    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
    using GridVolumeVariables = typename GridVariables::GridVolumeVariables;
    using ElementVolumeVariables = typename GridVolumeVariables::LocalView;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
 
    using FVElementGeometry = typename GridGeometry::LocalView;
    using SubControlVolume = typename FVElementGeometry::SubControlVolume;
    using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
    using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
    static constexpr int dim = GridView::dimension;
    static constexpr int dimWorld = GridView::dimensionworld;
 
    using GlobalPosition = typename SubControlVolumeFace::GlobalPosition;
    using VelocityVector = Dune::FieldVector<Scalar, dimWorld>;
    using GravityVector = Dune::FieldVector<Scalar, dimWorld>;
    using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
    using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
 
    static constexpr bool isVertexCentered
        = GridGeometry::discMethod == DiscretizationMethods::box
        || GridGeometry::discMethod == DiscretizationMethods::pq1bubble;
 
public:
    using InitialValues = Dune::FieldVector<Scalar, dimWorld>;
    using Sources = Dune::FieldVector<Scalar, dimWorld>;
    using DirichletValues = Dune::FieldVector<Scalar, dimWorld>;
    using BoundaryFluxes = Dune::FieldVector<Scalar, dimWorld>;
 
    using BoundaryTypes = NavierStokesMomentumBoundaryTypes<ModelTraits::dim()>;
 
    static constexpr bool isMomentumProblem() { return true; }
 
    NavierStokesMomentumProblemImpl(std::shared_ptr<const GridGeometry> gridGeometry,
                                    std::shared_ptr<CouplingManager> couplingManager,
                                    const std::string& paramGroup = "")
    : ParentType(gridGeometry, paramGroup)
    , gravity_(0.0)
    , couplingManager_(couplingManager)
    {
        if (getParamFromGroup<bool>(paramGroup, "Problem.EnableGravity"))
            gravity_[dim-1]  = -9.81;
 
        enableInertiaTerms_ = getParamFromGroup<bool>(paramGroup, "Problem.EnableInertiaTerms");
    }
 
    NavierStokesMomentumProblemImpl(std::shared_ptr<const GridGeometry> gridGeometry,
                                    const std::string& paramGroup = "")
    : NavierStokesMomentumProblemImpl(gridGeometry, {}, paramGroup)
    {}
 
    template<class ElementVolumeVariables>
    Sources 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());
    }
 
    Sources sourceAtPos(const GlobalPosition &globalPos) const
    {
        return Sources(0.0);
    }
 
    BoundaryTypes boundaryTypes(const Element& element,
                                const SubControlVolume& scv) const
    {
        if (!isVertexCentered)
            DUNE_THROW(Dune::InvalidStateException,
                       "boundaryTypes(..., scv) called for for a non vertex-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 (isVertexCentered)
            DUNE_THROW(Dune::InvalidStateException,
                       "boundaryTypes(..., scvf) called for a vertex-centered method.");
 
        // forward it to the method which only takes the global coordinate
        return asImp_().boundaryTypesAtPos(scvf.ipGlobal());
    }
 
    DirichletValues dirichlet(const Element& element, const SubControlVolume& scv) const
    {
        // forward it to the method which only takes the global coordinate
        if (!isVertexCentered)
            DUNE_THROW(Dune::InvalidStateException, "dirichlet(scv) called for a non vertex-centered method.");
        else
            return asImp_().dirichletAtPos(scv.dofPosition());
    }
 
    DirichletValues dirichlet(const Element& element, const SubControlVolumeFace& scvf) const
    {
        // forward it to the method which only takes the global coordinate
        if (isVertexCentered)
            DUNE_THROW(Dune::InvalidStateException, "dirichlet(scvf) called for vertex-centered method.");
        else
            return asImp_().dirichletAtPos(scvf.ipGlobal());
    }
 
    const GravityVector& gravity() const
    { return gravity_; }
 
    bool enableInertiaTerms() const
    { return enableInertiaTerms_; }
 
    Scalar referencePressure() const
    { return 0.0; }
 
    Scalar pressure(const Element& element,
                    const FVElementGeometry& fvGeometry,
                    const SubControlVolumeFace& scvf) const
    {
        if constexpr (std::is_empty_v<CouplingManager>)
            return asImp_().pressureAtPos(scvf.ipGlobal());
        else
            return couplingManager_->pressure(element, fvGeometry, scvf);
    }
 
    Scalar pressure(const Element& element,
                    const FVElementGeometry& fvGeometry,
                    const SubControlVolume& scv,
                    const bool isPreviousTimeStep = false) const
    {
        if constexpr (std::is_empty_v<CouplingManager>)
            return asImp_().pressureAtPos(scv.dofPosition());
        else
            return couplingManager_->pressure(element, fvGeometry, scv, isPreviousTimeStep);
    }
 
    Scalar pressureAtPos(const GlobalPosition&) const
    {
        DUNE_THROW(Dune::NotImplemented, "pressureAtPos not implemented");
    }
 
    Scalar density(const Element& element,
                   const FVElementGeometry& fvGeometry,
                   const SubControlVolumeFace& scvf) const
    {
        if constexpr (std::is_empty_v<CouplingManager>)
            return asImp_().densityAtPos(scvf.ipGlobal());
        else
            return couplingManager_->density(element, fvGeometry, scvf);
    }
 
    Scalar density(const Element& element,
                   const FVElementGeometry& fvGeometry,
                   const SubControlVolume& scv,
                   const bool isPreviousTimeStep = false) const
    {
        if constexpr (std::is_empty_v<CouplingManager>)
            return asImp_().densityAtPos(scv.dofPosition());
        else
            return couplingManager_->density(element, fvGeometry, scv, isPreviousTimeStep);
    }
 
 
    Scalar densityAtPos(const GlobalPosition&) const
    {
        DUNE_THROW(Dune::NotImplemented, "densityAtPos not implemented");
    }
 
    Scalar effectiveViscosity(const Element& element,
                              const FVElementGeometry& fvGeometry,
                              const SubControlVolumeFace& scvf) const
    {
        if constexpr (std::is_empty_v<CouplingManager>)
            return asImp_().effectiveViscosityAtPos(scvf.ipGlobal());
        else
            return couplingManager_->effectiveViscosity(element, fvGeometry, scvf);
    }
 
    Scalar effectiveViscosity(const Element& element,
                              const FVElementGeometry& fvGeometry,
                              const SubControlVolume& scv) const
    {
        if constexpr (std::is_empty_v<CouplingManager>)
            return asImp_().effectiveViscosityAtPos(scv.dofPosition());
        else
            return couplingManager_->effectiveViscosity(element, fvGeometry, scv);
    }
 
    Scalar effectiveViscosityAtPos(const GlobalPosition&) const
    {
        DUNE_THROW(Dune::NotImplemented, "effectiveViscosityAtPos not implemented");
    }
 
    template<class SolutionVector>
    void applyInitialSolution(SolutionVector& sol) const
    {
        static_assert(GridGeometry::discMethod == DiscretizationMethods::CVFE<DM>{});
        sol.resize(this->gridGeometry().numDofs());
        std::vector<bool> dofHandled(this->gridGeometry().numDofs(), false);
        auto fvGeometry = localView(this->gridGeometry());
        for (const auto& element : elements(this->gridGeometry().gridView()))
        {
            fvGeometry.bindElement(element);
            for (const auto& scv : scvs(fvGeometry))
            {
                const auto dofIdx = scv.dofIndex();
                if (!dofHandled[dofIdx])
                {
                    dofHandled[dofIdx] = true;
                    sol[dofIdx] = asImp_().initial(scv);
                }
            }
        }
    }
 
    InitialValues initial(const SubControlVolume& scv) const
    {
        static_assert(GridGeometry::discMethod == DiscretizationMethods::CVFE<DM>{});
        return asImp_().initialAtPos(scv.dofPosition());
    }
 
private:
    Implementation &asImp_()
    { return *static_cast<Implementation *>(this); }
 
    const Implementation &asImp_() const
    { return *static_cast<const Implementation *>(this); }
 
    GravityVector gravity_;
    bool enableInertiaTerms_;
    std::shared_ptr<CouplingManager> couplingManager_ = {};
};
 
template<class TypeTag>
using NavierStokesMomentumProblem = NavierStokesMomentumProblemImpl<
    TypeTag, typename GetPropType<TypeTag, Properties::GridGeometry>::DiscretizationMethod
>;
 
} // end namespace Dumux
 
#endif