freeflow/navierstokes/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_PROBLEM_HH
#define DUMUX_NAVIERSTOKES_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/common/staggeredfvproblem.hh>
#include <dumux/discretization/method.hh>
#include <dumux/freeflow/navierstokes/momentum/boundarytypes.hh>
 
namespace Dumux {
 
template<class TypeTag, class DiscretizationMethod> struct NavierStokesParentProblemImpl;
 
// compatibility with old-style Navier-Stokes models
template<class TypeTag>
struct NavierStokesParentProblemImpl<TypeTag, DiscretizationMethods::Staggered>
{
    using type = StaggeredFVProblem<TypeTag>;
};
 
template<class TypeTag>
using NavierStokesParentProblem = typename NavierStokesParentProblemImpl<
    TypeTag, typename GetPropType<TypeTag, Properties::GridGeometry>::DiscretizationMethod
>::type;
 
template<class TypeTag, class DiscretizationMethod>
class NavierStokesProblemImpl;
 
template<class TypeTag>
class NavierStokesProblemImpl<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; }
 
    NavierStokesProblemImpl(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");
    }
 
    NavierStokesProblemImpl(std::shared_ptr<const GridGeometry> gridGeometry,
                            const std::string& paramGroup = "")
    : NavierStokesProblemImpl(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 NavierStokesProblemImpl<TypeTag, DiscretizationMethods::CCTpfa>
: 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 FVElementGeometry = typename GridGeometry::LocalView;
    using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
    using Element = typename GridView::template Codim<0>::Entity;
    using GlobalPosition = typename SubControlVolumeFace::GlobalPosition;
    using VelocityVector = GlobalPosition;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    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, ModelTraits::numEq()>;
    using Sources = Dune::FieldVector<Scalar, ModelTraits::numEq()>;
    using DirichletValues = Dune::FieldVector<Scalar, ModelTraits::numEq()>;
    using BoundaryFluxes = Dune::FieldVector<Scalar, ModelTraits::numEq()>;
 
    using BoundaryTypes = Dumux::BoundaryTypes<ModelTraits::numEq()>;
 
    static constexpr bool isMomentumProblem() { return false; }
 
    NavierStokesProblemImpl(std::shared_ptr<const GridGeometry> gridGeometry,
                            std::shared_ptr<CouplingManager> couplingManager,
                            const std::string& paramGroup = "")
    : ParentType(gridGeometry, paramGroup)
    , couplingManager_(couplingManager)
    {}
 
    NavierStokesProblemImpl(std::shared_ptr<const GridGeometry> gridGeometry,
                            const std::string& paramGroup = "")
    : NavierStokesProblemImpl(gridGeometry, {}, paramGroup)
    {}
 
    VelocityVector faceVelocity(const Element& element,
                                const FVElementGeometry& fvGeometry,
                                const SubControlVolumeFace& scvf) const
    {
        if constexpr (isCoupled_)
            return couplingManager_->faceVelocity(element, scvf);
        else
            return asImp_().velocityAtPos(scvf.ipGlobal());
    }
 
    VelocityVector velocityAtPos(const GlobalPosition&) const
    {
        DUNE_THROW(Dune::NotImplemented, "velocityAtPos not implemented");
    }
 
    [[deprecated("temperature should now be defined in the spatial params with temperature(globalPos)")]]
    Scalar temperatureAtPos(const GlobalPosition& globalPos, int deprecationHelper = 0) const
    { return asImp_().temperature(); }
 
    [[deprecated("temperature should now be defined in the spatial params with temperature(globalPos)")]]
    Scalar temperature(int deprecationHelper = 0) const
    {
        DUNE_THROW(Dune::NotImplemented, "temperature() method not implemented by the actual problem");
    }
 
    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:
    Implementation &asImp_()
    { return *static_cast<Implementation *>(this); }
 
    const Implementation &asImp_() const
    { return *static_cast<const Implementation *>(this); }
 
    std::shared_ptr<CouplingManager> couplingManager_;
};
 
template<class TypeTag>
class NavierStokesProblemImpl<TypeTag, DiscretizationMethods::Staggered>
: public NavierStokesParentProblem<TypeTag>
{
    using ParentType = NavierStokesParentProblem<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 GridFaceVariables = typename GridVariables::GridFaceVariables;
    using ElementFaceVariables = typename GridFaceVariables::LocalView;
    using FaceVariables = typename GridFaceVariables::FaceVariables;
    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 PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
    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>;
 
public:
    NavierStokesProblemImpl(std::shared_ptr<const GridGeometry> gridGeometry, const std::string& paramGroup = "")
    : ParentType(gridGeometry, paramGroup)
    , gravity_(0.0)
    {
        if (getParamFromGroup<bool>(paramGroup, "Problem.EnableGravity"))
            gravity_[dim-1]  = -9.81;
 
        enableInertiaTerms_ = getParamFromGroup<bool>(paramGroup, "Problem.EnableInertiaTerms");
    }
 
    [[deprecated("temperature should now be defined in the spatial params with temperature(globalPos)")]]
    Scalar temperatureAtPos(const GlobalPosition &globalPos, int deprecationHelper = 0) const
    { return asImp_().temperature(); }
 
    [[deprecated("temperature should now be defined in the spatial params with temperature(globalPos)")]]
    Scalar temperature(int deprecationHelper = 0) const
    { DUNE_THROW(Dune::NotImplemented, "temperature() method not implemented by the actual problem"); }
 
    const GravityVector& gravity() const
    { return gravity_; }
 
    bool enableInertiaTerms() const
    { return enableInertiaTerms_; }
 
    template <class SolutionVector, class G = GridGeometry>
    typename std::enable_if<G::discMethod == DiscretizationMethods::staggered, void>::type
    applyInitialFaceSolution(SolutionVector& sol,
                             const SubControlVolumeFace& scvf,
                             const PrimaryVariables& initSol) const
    {
        sol[GridGeometry::faceIdx()][scvf.dofIndex()][0] = initSol[Indices::velocity(scvf.directionIndex())];
    }
 
    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);
    }
 
    template <class ElementVolumeVariables, class ElementFaceVariables, class G = GridGeometry>
    typename std::enable_if<G::discMethod == DiscretizationMethods::staggered, Scalar>::type
    pseudo3DWallFriction(const SubControlVolumeFace& scvf,
                         const ElementVolumeVariables& elemVolVars,
                         const ElementFaceVariables& elemFaceVars,
                         const Scalar height,
                         const Scalar factor = 8.0) const
    {
        const Scalar velocity = elemFaceVars[scvf].velocitySelf();
        const Scalar viscosity = elemVolVars[scvf.insideScvIdx()].effectiveViscosity();
        return pseudo3DWallFriction(velocity, viscosity, height, factor);
    }
 
    Scalar permeability(const Element& element, 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 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 Element& element, const SubControlVolumeFace& scvf, const GlobalPosition& tangentialVector) const
    {
        const Scalar interfacePermeability = interfacePermeability_(element, scvf, tangentialVector);
        using std::sqrt;
        return asImp_().alphaBJ(scvf) / sqrt(interfacePermeability);
    }
 
    VelocityVector porousMediumVelocity(const Element& element, const SubControlVolumeFace& scvf) const
    {
        return VelocityVector(0.0);
    }
 
    const Scalar beaversJosephVelocity(const Element& element,
                                       const SubControlVolume& scv,
                                       const SubControlVolumeFace& ownScvf,
                                       const SubControlVolumeFace& faceOnPorousBoundary,
                                       const Scalar velocitySelf,
                                       const Scalar tangentialVelocityGradient) const
    {
        // create a unit normal vector oriented in positive coordinate direction
        GlobalPosition orientation = ownScvf.unitOuterNormal();
        orientation[ownScvf.directionIndex()] = 1.0;
 
        // du/dy + dv/dx = alpha/sqrt(K) * (u_boundary-uPM)
        // beta = alpha/sqrt(K)
        const Scalar betaBJ = asImp_().betaBJ(element, faceOnPorousBoundary, orientation);
        const Scalar distanceNormalToBoundary = (faceOnPorousBoundary.center() - scv.center()).two_norm();
 
        return (tangentialVelocityGradient*distanceNormalToBoundary
              + asImp_().porousMediumVelocity(element, faceOnPorousBoundary) * orientation * betaBJ * distanceNormalToBoundary
              + velocitySelf) / (betaBJ*distanceNormalToBoundary + 1.0);
    }
 
private:
    Scalar interfacePermeability_(const Element& element, const SubControlVolumeFace& scvf, const GlobalPosition& tangentialVector) const
    {
        const auto& K = asImp_().permeability(element, 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_;
};
 
 
template<class TypeTag>
using NavierStokesProblem = NavierStokesProblemImpl<
    TypeTag, typename GetPropType<TypeTag, Properties::GridGeometry>::DiscretizationMethod
>;
 
} // end namespace Dumux
 
#endif