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:
//
// SPDX-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
// SPDX-License-Identifier: GPL-3.0-or-later
//
#ifndef DUMUX_NAVIERSTOKES_MOMENTUM_PROBLEM_HH
#define DUMUX_NAVIERSTOKES_MOMENTUM_PROBLEM_HH
 
#include <dune/common/exceptions.hh>
#include <dune/common/fvector.hh>
#include <dune/common/typetraits.hh>
#include <dumux/common/numeqvector.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/concepts/variables_.hh>
#include <dumux/discretization/cvfe/quadraturerules.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);
    }
 
    Sources scvPointSources(const Element& element,
                            const FVElementGeometry& fvGeometry,
                            const ElementVolumeVariables& elemVolVars,
                            const SubControlVolume& scv) const
    {
        if (!this->pointSourceMap().empty())
            DUNE_THROW(Dune::NotImplemented, "scvPointSources not implemented");
 
        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 effectiveViscosity(const Element& element,
                              const FVElementGeometry& fvGeometry,
                              const SubControlVolume& scv) const
    {
        if constexpr (isCoupled_)
            return couplingManager_->effectiveViscosity(element, fvGeometry, scv);
        else
            return asImp_().effectiveViscosityAtPos(scv.dofPosition());
    }
 
    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; }
 
    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); }
 
    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>;
 
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
    {
        // 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
    {
        DUNE_THROW(Dune::InvalidStateException, "boundaryTypes(..., scvf) called for a CVFE method.");
    }
 
    DirichletValues dirichlet(const Element& element, const SubControlVolume& scv) const
    {
        // forward it to the method which only takes the global coordinate
        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
        DUNE_THROW(Dune::InvalidStateException, "dirichlet(scvf) called for CVFE method.");
    }
 
    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);
    }
 
    template<class IpData>
    Scalar pressure(const Element& element,
                    const FVElementGeometry& fvGeometry,
                    const IpData& ipData,
                    const bool isPreviousTimeStep = false) const
    {
        if constexpr (std::is_empty_v<CouplingManager>)
            return asImp_().pressureAtPos(ipData.global());
        else
            return couplingManager_->pressure(element, fvGeometry, ipData, 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);
    }
 
    template<class IpData>
    Scalar density(const Element& element,
                   const FVElementGeometry& fvGeometry,
                   const IpData& ipData,
                   const bool isPreviousTimeStep = false) const
    {
        if constexpr (std::is_empty_v<CouplingManager>)
            return asImp_().densityAtPos(ipData.global());
        else
            return couplingManager_->density(element, fvGeometry, ipData, 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 bool isPreviousTimeStep = false) const
    {
        if constexpr (std::is_empty_v<CouplingManager>)
            return asImp_().effectiveViscosityAtPos(scv.dofPosition());
        else
            return couplingManager_->effectiveViscosity(element, fvGeometry, scv, isPreviousTimeStep);
    }
 
    template<class IpData>
    Scalar effectiveViscosity(const Element& element,
                              const FVElementGeometry& fvGeometry,
                              const IpData& ipData,
                              const bool isPreviousTimeStep = false) const
    {
        if constexpr (std::is_empty_v<CouplingManager>)
            return asImp_().effectiveViscosityAtPos(ipData.global());
        else
            return couplingManager_->effectiveViscosity(element, fvGeometry, ipData, isPreviousTimeStep);
    }
 
    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>
class CVFENavierStokesMomentumProblem
: 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 GridVariablesCache = Concept::GridVariablesCache_t<GridVariables>;
    using ElementVariables = typename GridVariablesCache::LocalView;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
 
    using FVElementGeometry = typename GridGeometry::LocalView;
    using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
    using GlobalPosition = typename SubControlVolumeFace::GlobalPosition;
    using LocalPosition = typename Element::Geometry::LocalCoordinate;
 
    static constexpr int dim = GridView::dimension;
    static constexpr int dimWorld = GridView::dimensionworld;
 
    using GravityVector = Dune::FieldVector<Scalar, dimWorld>;
    using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
    using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
 
public:
    using InitialValues = Dune::FieldVector<Scalar, ModelTraits::numEq()>;
    using Sources = Dune::FieldVector<Scalar, ModelTraits::numEq()>;
    using DirichletValues = Dune::FieldVector<Scalar, dimWorld>;
    using BoundaryFluxes = Dune::FieldVector<Scalar, ModelTraits::numEq()>;
 
    using BoundaryTypes = NavierStokesMomentumBoundaryTypes<ModelTraits::dim()>;
 
    static constexpr bool isMomentumProblem() { return true; }
    static constexpr bool providesIntegralInterface() { return true; }
 
    CVFENavierStokesMomentumProblem(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");
    }
 
    CVFENavierStokesMomentumProblem(std::shared_ptr<const GridGeometry> gridGeometry,
                                    const std::string& paramGroup = "")
    : CVFENavierStokesMomentumProblem(gridGeometry, {}, paramGroup)
    {}
 
    template<class IpData>
    Sources source(const FVElementGeometry& fvGeometry,
                   const ElementVariables& elemVars,
                   const IpData& ipData) const
    {
        return asImp_().sourceAtPos(ipData.global());
    }
 
    Sources sourceAtPos(const GlobalPosition &globalPos) const
    {
        return Sources(0.0);
    }
 
    template<class Intersection>
    BoundaryTypes boundaryTypes(const FVElementGeometry& fvGeometry,
                                const Intersection& intersection) const
    {
        // forward it to the method which only takes the global coordinate
        return asImp_().boundaryTypesAtPos(intersection.geometry().center());
    }
 
    template<class FaceIpData>
    DirichletValues dirichlet(const FVElementGeometry& fvGeometry,
                              const FaceIpData& faceIpData) const
    {
        // forward it to the method which only takes the global coordinate
        return asImp_().dirichletAtPos(faceIpData.global());
    }
 
    BoundaryFluxes boundaryFluxIntegral(const FVElementGeometry& fvGeometry,
                                        const ElementVariables& elemVars,
                                        const SubControlVolumeFace& scvf) const
    {
        BoundaryFluxes flux(0.0);
        for (const auto& qpData : CVFE::quadratureRule(fvGeometry, scvf))
            flux += qpData.weight() * asImp_().boundaryFlux(fvGeometry, elemVars, qpData.ipData());
 
        return flux * elemVars[fvGeometry.scv(scvf.insideScvIdx())].extrusionFactor();
    }
 
    template<class ElementFluxVariablesCache>
    [[deprecated("This function is deprecated and will be removed after release 3.11. "
                 "Use boundaryFluxIntegral without elemFluxVarsCache instead.")]]
    BoundaryFluxes boundaryFluxIntegral(const FVElementGeometry& fvGeometry,
                                        const ElementVariables& elemVars,
                                        const ElementFluxVariablesCache& elemFluxVarsCache,
                                        const SubControlVolumeFace& scvf) const
    {
        BoundaryFluxes flux(0.0);
        for (const auto& qpData : CVFE::quadratureRule(fvGeometry, scvf))
            flux += qpData.weight() * asImp_().boundaryFlux(fvGeometry, elemVars, elemFluxVarsCache, qpData.ipData());
 
        return flux * elemVars[fvGeometry.scv(scvf.insideScvIdx())].extrusionFactor();
    }
 
    template<class ResidualVector, class BoundaryTypes>
    void addBoundaryFluxIntegrals(ResidualVector& residual,
                                  const FVElementGeometry& fvGeometry,
                                  const ElementVariables& elemVars,
                                  const typename FVElementGeometry::GridGeometry::GridView::Intersection& intersection,
                                  const BoundaryTypes& bcTypes) const
    {
        // quadrature rule for intersections (dim-1)
        for (const auto& qpData : CVFE::quadratureRule(fvGeometry, intersection))
        {
            const auto& ipData = qpData.ipData();
            const auto& ipCache = cache(elemVars, ipData);
            for (const auto& localDof : nonCVLocalDofs(fvGeometry))
            {
                const BoundaryFluxes& boundaryFlux = qpData.weight()*asImp_().boundaryFlux(fvGeometry, elemVars, ipData);
                const auto& shapeValues = ipCache.shapeValues();
                // only add fluxes to equations for which Neumann is set
                for (int eqIdx = 0; eqIdx < BoundaryFluxes::dimension; ++eqIdx)
                    if (bcTypes.isNeumann(eqIdx))
                        residual[localDof.index()][eqIdx] += shapeValues[localDof.index()] * boundaryFlux[eqIdx];
            }
        }
    }
 
    template<class ResidualVector, class ElementFluxVariablesCache, class BoundaryTypes>
    [[deprecated("This function is deprecated and will be removed after release 3.11. "
                 "Use addBoundaryFluxIntegrals without elemFluxVarsCache instead.")]]
    void addBoundaryFluxIntegrals(ResidualVector& residual,
                                  const FVElementGeometry& fvGeometry,
                                  const ElementVariables& elemVars,
                                  const ElementFluxVariablesCache& elemFluxVarsCache,
                                  const typename FVElementGeometry::GridGeometry::GridView::Intersection& intersection,
                                  const BoundaryTypes& bcTypes) const
    {
        // quadrature rule for intersections (dim-1)
        for (const auto& qpData : CVFE::quadratureRule(fvGeometry, intersection))
        {
            const auto& ipData = qpData.ipData();
            const auto& cache = elemFluxVarsCache[ipData];
            for (const auto& localDof : nonCVLocalDofs(fvGeometry))
            {
                const BoundaryFluxes& boundaryFlux = qpData.weight()*asImp_().boundaryFlux(fvGeometry, elemVars, elemFluxVarsCache, ipData);
                const auto& shapeValues = cache.shapeValues();
                // only add fluxes to equations for which Neumann is set
                for (int eqIdx = 0; eqIdx < BoundaryFluxes::dimension; ++eqIdx)
                    if (bcTypes.isNeumann(eqIdx))
                        residual[localDof.index()][eqIdx] += shapeValues[localDof.index()] * boundaryFlux[eqIdx];
            }
        }
    }
 
    template<class FaceIpData>
    BoundaryFluxes boundaryFlux(const FVElementGeometry& fvGeometry,
                                const ElementVariables& elemVars,
                                const FaceIpData& faceIpData) const
    {
        return asImp_().boundaryFluxAtPos(faceIpData.global());
    }
 
    template<class ElementFluxVariablesCache, class FaceIpData>
    [[deprecated("This function is deprecated and will be removed after release 3.11. "
                 "Use boundaryFlux without elemFluxVarsCache instead.")]]
    BoundaryFluxes boundaryFlux(const FVElementGeometry& fvGeometry,
                                const ElementVariables& elemVars,
                                const ElementFluxVariablesCache& elemFluxVarsCache,
                                const FaceIpData& faceIpData) const
    {
        return asImp_().boundaryFluxAtPos(faceIpData.global());
    }
 
    BoundaryFluxes boundaryFluxAtPos(const GlobalPosition& globalPos) const
    { return BoundaryFluxes(0.0); } 
 
    const GravityVector& gravity() const
    { return gravity_; }
 
    bool enableInertiaTerms() const
    { return enableInertiaTerms_; }
 
    Scalar referencePressure() const
    { return 0.0; }
 
    template<class IpData>
    Scalar pressure(const Element& element,
                    const FVElementGeometry& fvGeometry,
                    const IpData& ipData,
                    const bool isPreviousTimeStep = false) const
    {
        if constexpr (std::is_empty_v<CouplingManager>)
            return asImp_().pressureAtPos(ipData.global());
        else
            return couplingManager_->pressure(element, fvGeometry, ipData, isPreviousTimeStep);
    }
 
    Scalar pressureAtPos(const GlobalPosition&) const
    {
        DUNE_THROW(Dune::NotImplemented, "pressureAtPos not implemented");
    }
 
    template<class IpData>
    Scalar density(const Element& element,
                   const FVElementGeometry& fvGeometry,
                   const IpData& ipData,
                   const bool isPreviousTimeStep = false) const
    {
        if constexpr (std::is_empty_v<CouplingManager>)
            return asImp_().densityAtPos(ipData.global());
        else
            return couplingManager_->density(element, fvGeometry, ipData, isPreviousTimeStep);
    }
 
    Scalar densityAtPos(const GlobalPosition&) const
    {
        DUNE_THROW(Dune::NotImplemented, "densityAtPos not implemented");
    }
 
 
    template<class IpData>
    Scalar effectiveViscosity(const Element& element,
                              const FVElementGeometry& fvGeometry,
                              const IpData& ipData,
                              const bool isPreviousTimeStep = false) const
    {
        if constexpr (std::is_empty_v<CouplingManager>)
            return asImp_().effectiveViscosityAtPos(ipData.global());
        else
            return couplingManager_->effectiveViscosity(element, fvGeometry, ipData, isPreviousTimeStep);
    }
 
    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& localDof : localDofs(fvGeometry))
            {
                const auto dofIdx = localDof.dofIndex();
                if (!dofHandled[dofIdx])
                {
                    dofHandled[dofIdx] = true;
                    sol[dofIdx] = asImp_().initial(fvGeometry, ipData(fvGeometry, localDof));
                }
            }
        }
    }
 
    template<class IpData>
    InitialValues initial(const FVElementGeometry& fvGeometry, const IpData& ipData) const
    {
        return asImp_().initialAtPos(ipData.global());
    }
 
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