flux.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_CVFE_FLUXVARIABLES_HH
#define DUMUX_NAVIERSTOKES_MOMENTUM_CVFE_FLUXVARIABLES_HH
 
#include <dune/common/fmatrix.hh>
 
#include <dumux/common/math.hh>
#include <dumux/common/exceptions.hh>
#include <dumux/common/parameters.hh>
 
#include <dumux/discretization/extrusion.hh>
 
namespace Dumux {
 
template<class Problem,
         class FVElementGeometry,
         class ElementVolumeVariables,
         class ElementFluxVariablesCache>
class NavierStokesMomentumFluxContext
{
    using Element = typename FVElementGeometry::Element;
    using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
public:
 
    NavierStokesMomentumFluxContext(
        const Problem& problem,
        const FVElementGeometry& fvGeometry,
        const ElementVolumeVariables& elemVolVars,
        const ElementFluxVariablesCache& elemFluxVarsCache,
        const SubControlVolumeFace& scvf
    )
    : problem_(problem)
    , fvGeometry_(fvGeometry)
    , elemVolVars_(elemVolVars)
    , elemFluxVarsCache_(elemFluxVarsCache)
    , scvf_(scvf)
    {}
 
    const Problem& problem() const
    { return problem_; }
 
    const Element& element() const
    { return fvGeometry_.element(); }
 
    const SubControlVolumeFace& scvFace() const
    { return scvf_; }
 
    const FVElementGeometry& fvGeometry() const
    { return fvGeometry_; }
 
    const ElementVolumeVariables& elemVolVars() const
    { return elemVolVars_; }
 
    const ElementFluxVariablesCache& elemFluxVarsCache() const
    { return elemFluxVarsCache_; }
 
private:
    const Problem& problem_;
    const FVElementGeometry& fvGeometry_;
    const ElementVolumeVariables& elemVolVars_;
    const ElementFluxVariablesCache& elemFluxVarsCache_;
    const SubControlVolumeFace& scvf_;
};
 
template<class Problem,
         class FVElementGeometry,
         class ElementVolumeVariables,
         class IpData>
class NavierStokesMomentumFluxFunctionContext
{
    using Element = typename FVElementGeometry::Element;
    using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
    using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
 
    static constexpr int dim = FVElementGeometry::GridGeometry::GridView::dimension;
    static constexpr int dimWorld = FVElementGeometry::GridGeometry::GridView::dimensionworld;
 
    using Tensor = Dune::FieldMatrix<typename GlobalPosition::value_type, dim, dimWorld>;
 
public:
 
    NavierStokesMomentumFluxFunctionContext(
        const Problem& problem,
        const FVElementGeometry& fvGeometry,
        const ElementVolumeVariables& elemVolVars,
        const IpData& ipData
    )
    : problem_(problem)
    , fvGeometry_(fvGeometry)
    , elemVolVars_(elemVolVars)
    , ipData_(ipData)
    , velocity_(0.0)
    , gradVelocity_(0.0)
    {
        calculateVelocity_();
        calculateGradVelocity();
    }
 
    const Problem& problem() const
    { return problem_; }
 
    const Element& element() const
    { return fvGeometry_.element(); }
 
    const FVElementGeometry& fvGeometry() const
    { return fvGeometry_; }
 
    const ElementVolumeVariables& elemVolVars() const
    { return elemVolVars_; }
 
    const GlobalPosition& velocity() const
    { return velocity_; }
 
    const Tensor& gradVelocity() const
    { return gradVelocity_; }
 
private:
    void calculateVelocity_()
    {
        const auto& shapeValues = ipData_.shapeValues();
        for (const auto& localDof : localDofs(fvGeometry_))
            velocity_.axpy(shapeValues[localDof.index()][0], elemVolVars_[localDof.index()].velocity());
    }
 
    void calculateGradVelocity()
    {
        for (const auto& localDof : localDofs(fvGeometry_))
        {
            const auto& volVars = elemVolVars_[localDof.index()];
            for (int dir = 0; dir < dim; ++dir)
                gradVelocity_[dir].axpy(volVars.velocity(dir), ipData_.gradN(localDof.index()));
        }
    }
 
    const Problem& problem_;
    const FVElementGeometry& fvGeometry_;
    const ElementVolumeVariables& elemVolVars_;
    const IpData& ipData_;
    GlobalPosition velocity_;
    Tensor gradVelocity_;
};
 
template<class GridGeometry, class NumEqVector>
class NavierStokesMomentumFluxCVFE
{
    using GridView = typename GridGeometry::GridView;
    using Element = typename GridView::template Codim<0>::Entity;
    using Scalar = typename NumEqVector::value_type;
 
    using Extrusion = Extrusion_t<GridGeometry>;
 
    static constexpr int dim = GridView::dimension;
    static constexpr int dimWorld = GridView::dimensionworld;
 
    using Tensor = Dune::FieldMatrix<Scalar, dim, dimWorld>;
    static_assert(NumEqVector::dimension == dimWorld, "Wrong dimension of velocity vector");
 
public:
    template<class Context>
    NumEqVector advectiveMomentumFlux(const Context& context) const
    {
        if (!context.problem().enableInertiaTerms())
            return NumEqVector(0.0);
 
        const auto& fvGeometry = context.fvGeometry();
        const auto& elemVolVars = context.elemVolVars();
        const auto& scvf = context.scvFace();
        const auto& fluxVarCache = context.elemFluxVarsCache()[scvf];
        const auto& shapeValues = fluxVarCache.shapeValues();
 
        // interpolate velocity at scvf
        NumEqVector v(0.0);
        for (const auto& localDof : localDofs(fvGeometry))
            v.axpy(shapeValues[localDof.index()][0], elemVolVars[localDof.index()].velocity());
 
        // get density from the problem
        const Scalar density = context.problem().density(context.element(), context.fvGeometry(), fluxVarCache.ipData());
 
        const auto vn = v*scvf.unitOuterNormal();
        const auto& insideVolVars = elemVolVars[fvGeometry.scv(scvf.insideScvIdx())];
        const auto& outsideVolVars = elemVolVars[fvGeometry.scv(scvf.outsideScvIdx())];
        const auto upwindVelocity = vn > 0 ? insideVolVars.velocity() : outsideVolVars.velocity();
        const auto downwindVelocity = vn > 0 ? outsideVolVars.velocity() : insideVolVars.velocity();
        static const auto upwindWeight = getParamFromGroup<Scalar>(context.problem().paramGroup(), "Flux.UpwindWeight");
        const auto advectiveTermIntegrand = density*vn * (upwindWeight * upwindVelocity + (1.0-upwindWeight)*downwindVelocity);
 
        return advectiveTermIntegrand * Extrusion::area(fvGeometry, scvf) * insideVolVars.extrusionFactor();
    }
 
    template<class Context>
    NumEqVector diffusiveMomentumFlux(const Context& context) const
    {
        const auto& element = context.element();
        const auto& fvGeometry = context.fvGeometry();
        const auto& elemVolVars = context.elemVolVars();
        const auto& scvf = context.scvFace();
        const auto& fluxVarCache = context.elemFluxVarsCache()[scvf];
 
        // interpolate velocity gradient at scvf
        Tensor gradV(0.0);
        for (const auto& localDof : localDofs(fvGeometry))
        {
            const auto& volVars = elemVolVars[localDof];
            for (int dir = 0; dir < dim; ++dir)
                gradV[dir].axpy(volVars.velocity(dir), fluxVarCache.gradN(localDof.index()));
        }
 
        // get viscosity from the problem
        const auto mu = context.problem().effectiveViscosity(element, fvGeometry, fluxVarCache.ipData());
 
        static const bool enableUnsymmetrizedVelocityGradient
            = getParamFromGroup<bool>(context.problem().paramGroup(), "FreeFlow.EnableUnsymmetrizedVelocityGradient", false);
 
        // compute -mu*gradV*n*dA
        NumEqVector diffusiveFlux = enableUnsymmetrizedVelocityGradient ?
                mv(gradV, scvf.unitOuterNormal())
                : mv(gradV + getTransposed(gradV),scvf.unitOuterNormal());
 
        diffusiveFlux *= -mu;
 
        static const bool enableDilatationTerm = getParamFromGroup<bool>(context.problem().paramGroup(), "FreeFlow.EnableDilatationTerm", false);
        if (enableDilatationTerm)
            diffusiveFlux += 2.0/3.0 * mu * trace(gradV) * scvf.unitOuterNormal();
 
        diffusiveFlux *= Extrusion::area(fvGeometry, scvf) * elemVolVars[fvGeometry.scv(scvf.insideScvIdx())].extrusionFactor();
        return diffusiveFlux;
    }
 
    template<class Context>
    NumEqVector pressureContribution(const Context& context) const
    {
        const auto& element = context.element();
        const auto& fvGeometry = context.fvGeometry();
        const auto& elemVolVars = context.elemVolVars();
        const auto& scvf = context.scvFace();
        const auto& fluxVarCache = context.elemFluxVarsCache()[scvf];
 
        // The pressure force needs to take the extruded scvf area into account
        const auto pressure = context.problem().pressure(element, fvGeometry, fluxVarCache.ipData());
 
        // The pressure contribution calculated above might have a much larger numerical value compared to the viscous or inertial forces.
        // This may lead to numerical inaccuracies due to loss of significance (cancellation) for the final residual value.
        // In the end, we are only interested in a pressure gradient between the two relevant faces so we can
        // subtract a constant reference value from the actual pressure contribution.
        const auto referencePressure = context.problem().referencePressure();
 
        NumEqVector pn(scvf.unitOuterNormal());
        pn *= (pressure-referencePressure)*Extrusion::area(fvGeometry, scvf)*elemVolVars[fvGeometry.scv(scvf.insideScvIdx())].extrusionFactor();
 
        return pn;
    }
};
 
template<class GridGeometry, class NumEqVector>
class NavierStokesMomentumFluxFunctionCVFE
{
    using GridView = typename GridGeometry::GridView;
    using Element = typename GridView::template Codim<0>::Entity;
    using Scalar = typename NumEqVector::value_type;
 
    using Extrusion = Extrusion_t<GridGeometry>;
 
    static constexpr int dim = GridView::dimension;
    static constexpr int dimWorld = GridView::dimensionworld;
 
    using Tensor = Dune::FieldMatrix<Scalar, dim, dimWorld>;
    static_assert(NumEqVector::dimension == dimWorld, "Wrong dimension of velocity vector");
 
public:
    template<class Problem, class FVElementGeometry, class ElementVariables,
             class ElementFluxVariablesCache, class SubControlVolumeFace, class VelocityVector>
    NumEqVector advectiveMomentumFluxIntegral(const Problem& problem,
                                              const FVElementGeometry& fvGeometry,
                                              const ElementVariables& elemVars,
                                              const ElementFluxVariablesCache& elemFluxVarsCache,
                                              const SubControlVolumeFace& scvf,
                                              const VelocityVector& integratedVelocity) const
    {
        if (!problem.enableInertiaTerms())
            return NumEqVector(0.0);
 
        // get density from the problem
        const Scalar density = problem.density(fvGeometry.element(), fvGeometry, ipData(fvGeometry, scvf));
 
        const auto vn_integral = integratedVelocity*scvf.unitOuterNormal();
        const auto& insideVolVars = elemVars[fvGeometry.scv(scvf.insideScvIdx())];
        const auto& outsideVolVars = elemVars[fvGeometry.scv(scvf.outsideScvIdx())];
        const auto upwindVelocity = vn_integral > 0 ? insideVolVars.velocity() : outsideVolVars.velocity();
        const auto downwindVelocity = vn_integral > 0 ? outsideVolVars.velocity() : insideVolVars.velocity();
        static const auto upwindWeight = getParamFromGroup<Scalar>(problem.paramGroup(), "Flux.UpwindWeight");
        const auto advectiveFlux = density*vn_integral * (upwindWeight * upwindVelocity + (1.0-upwindWeight)*downwindVelocity);
 
        return advectiveFlux;
    }
 
    template<class Context, class IpData>
    NumEqVector diffusiveMomentumFluxIntegrand(const Context& context, const IpData& ipData) const
    {
        const auto& element = context.element();
        const auto& fvGeometry = context.fvGeometry();
 
        // get viscosity from the problem
        const auto mu = context.problem().effectiveViscosity(element, fvGeometry, ipData);
 
        static const bool enableUnsymmetrizedVelocityGradient
            = getParamFromGroup<bool>(context.problem().paramGroup(), "FreeFlow.EnableUnsymmetrizedVelocityGradient", false);
 
        // compute -mu*gradV*n*dA
        const auto& gradV = context.gradVelocity();
        NumEqVector diffusiveFluxIntegrand = enableUnsymmetrizedVelocityGradient ?
                mv(gradV, ipData.unitOuterNormal())
                : mv(gradV + getTransposed(gradV), ipData.unitOuterNormal());
 
        diffusiveFluxIntegrand *= -mu;
 
        static const bool enableDilatationTerm = getParamFromGroup<bool>(context.problem().paramGroup(), "FreeFlow.EnableDilatationTerm", false);
        if (enableDilatationTerm)
            diffusiveFluxIntegrand += 2.0/3.0 * mu * trace(gradV) * ipData.unitOuterNormal();
 
        return diffusiveFluxIntegrand;
    }
 
    template<class Context, class IpData>
    NumEqVector pressureFluxIntegrand(const Context& context, const IpData& ipData) const
    {
        const auto& element = context.element();
        const auto& fvGeometry = context.fvGeometry();
 
        // The pressure force integrand
        const auto pressure = context.problem().pressure(element, fvGeometry, ipData);
 
        // The pressure contribution calculated above might have a much larger numerical value compared to the viscous or inertial forces.
        // This may lead to numerical inaccuracies due to loss of significance (cancellation) for the final residual value.
        // Subtracting a constant reference value from the actual pressure contribution helps improving cancellation issues.
        const auto referencePressure = context.problem().referencePressure();
 
        NumEqVector pn(ipData.unitOuterNormal());
        pn *= (pressure-referencePressure);
 
        return pn;
    }
};
 
} // end namespace Dumux
 
#endif