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 ElementDiscretization,
         class ElementVariables,
         class ElementFluxVariablesCache>
class NavierStokesMomentumFluxContext
{
    using Element = typename ElementDiscretization::Element;
    using SubControlVolumeFace = typename ElementDiscretization::SubControlVolumeFace;
public:

    NavierStokesMomentumFluxContext(
        const Problem& problem,
        const ElementDiscretization& elemDisc,
        const ElementVariables& elemVars,
        const ElementFluxVariablesCache& elemFluxVarsCache,
        const SubControlVolumeFace& scvf
    )
    : problem_(problem)
    , elemDisc_(elemDisc)
    , elemVars_(elemVars)
    , elemFluxVarsCache_(elemFluxVarsCache)
    , scvf_(scvf)
    {}
 
    const Problem& problem() const
    { return problem_; }
 
    const Element& element() const
    { return elemDisc_.element(); }
 
    const SubControlVolumeFace& scvFace() const
    { return scvf_; }
 
    const ElementDiscretization& elemDisc() const
    { return elemDisc_; }
 
    const ElementVariables& elemVars() const
    { return elemVars_; }
 
    const ElementFluxVariablesCache& elemFluxVarsCache() const
    { return elemFluxVarsCache_; }
 
private:
    const Problem& problem_;
    const ElementDiscretization& elemDisc_;
    const ElementVariables& elemVars_;
    const ElementFluxVariablesCache& elemFluxVarsCache_;
    const SubControlVolumeFace& scvf_;
};

template<class Problem,
         class ElementDiscretization,
         class ElementVariables,
         class IpData>
class NavierStokesMomentumFluxFunctionContext
{
    using Element = typename ElementDiscretization::Element;
    using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
 
    static constexpr int dim = ElementDiscretization::GridGeometry::GridView::dimension;
    static constexpr int dimWorld = ElementDiscretization::GridGeometry::GridView::dimensionworld;
 
    using Tensor = Dune::FieldMatrix<typename GlobalPosition::value_type, dim, dimWorld>;
 
public:

    NavierStokesMomentumFluxFunctionContext(
        const Problem& problem,
        const ElementDiscretization& elemDisc,
        const ElementVariables& elemVars,
        const IpData& ipData
    )
    : problem_(problem)
    , elemDisc_(elemDisc)
    , elemVars_(elemVars)
    , ipData_(ipData)
    , velocity_(0.0)
    , gradVelocity_(0.0)
    {
        calculateVelocity_();
        calculateGradVelocity();
    }
 
    const Problem& problem() const
    { return problem_; }
 
    const Element& element() const
    { return elemDisc_.element(); }
 
    const ElementDiscretization& elemDisc() const
    { return elemDisc_; }
 
    const ElementVariables& elemVars() const
    { return elemVars_; }
 
    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(elemDisc_))
            velocity_.axpy(shapeValues[localDof.index()][0], elemVars_[localDof.index()].velocity());
    }
 
    void calculateGradVelocity()
    {
        for (const auto& localDof : localDofs(elemDisc_))
        {
            const auto& vars = elemVars_[localDof.index()];
            for (int dir = 0; dir < dim; ++dir)
                gradVelocity_[dir].axpy(vars.velocity(dir), ipData_.gradN(localDof.index()));
        }
    }
 
    const Problem& problem_;
    const ElementDiscretization& elemDisc_;
    const ElementVariables& elemVars_;
    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& elemDisc = context.elemDisc();
        const auto& elemVars = context.elemVars();
        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(elemDisc))
            v.axpy(shapeValues[localDof.index()][0], elemVars[localDof.index()].velocity());
 
        // get density from the problem
        const Scalar density = context.problem().density(context.element(), context.elemDisc(), fluxVarCache.ipData());
 
        const auto vn = v*scvf.unitOuterNormal();
        const auto& insideVars = elemVars[elemDisc.scv(scvf.insideScvIdx())];
        const auto& outsideVars = elemVars[elemDisc.scv(scvf.outsideScvIdx())];
        const auto upwindVelocity = vn > 0 ? insideVars.velocity() : outsideVars.velocity();
        const auto downwindVelocity = vn > 0 ? outsideVars.velocity() : insideVars.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(elemDisc, scvf) * insideVars.extrusionFactor();
    }

    template<class Context>
    NumEqVector diffusiveMomentumFlux(const Context& context) const
    {
        const auto& element = context.element();
        const auto& elemDisc = context.elemDisc();
        const auto& elemVars = context.elemVars();
        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(elemDisc))
        {
            const auto& vars = elemVars[localDof];
            for (int dir = 0; dir < dim; ++dir)
                gradV[dir].axpy(vars.velocity(dir), fluxVarCache.gradN(localDof.index()));
        }
 
        // get viscosity from the problem
        const auto mu = context.problem().effectiveViscosity(element, elemDisc, 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(elemDisc, scvf) * elemVars[elemDisc.scv(scvf.insideScvIdx())].extrusionFactor();
        return diffusiveFlux;
    }
 
    template<class Context>
    NumEqVector pressureContribution(const Context& context) const
    {
        const auto& element = context.element();
        const auto& elemDisc = context.elemDisc();
        const auto& elemVars = context.elemVars();
        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, elemDisc, 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(elemDisc, scvf)*elemVars[elemDisc.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 ElementDiscretization, class ElementVariables,
             class SubControlVolumeFace, class VelocityVector>
    NumEqVector advectiveMomentumFluxIntegral(const Problem& problem,
                                              const ElementDiscretization& elemDisc,
                                              const ElementVariables& elemVars,
                                              const SubControlVolumeFace& scvf,
                                              const VelocityVector& integratedVelocity) const
    {
        if (!problem.enableInertiaTerms())
            return NumEqVector(0.0);
 
        // get density from the problem
        const Scalar density = problem.density(elemDisc.element(), elemDisc, ipData(elemDisc, scvf));
 
        const auto vn_integral = integratedVelocity*scvf.unitOuterNormal();
        const auto& insideVars = elemVars[elemDisc.scv(scvf.insideScvIdx())];
        const auto& outsideVars = elemVars[elemDisc.scv(scvf.outsideScvIdx())];
        const auto upwindVelocity = vn_integral > 0 ? insideVars.velocity() : outsideVars.velocity();
        const auto downwindVelocity = vn_integral > 0 ? outsideVars.velocity() : insideVars.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& elemDisc = context.elemDisc();
 
        // get viscosity from the problem
        const auto mu = context.problem().effectiveViscosity(element, elemDisc, 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& elemDisc = context.elemDisc();
 
        // The pressure force integrand
        const auto pressure = context.problem().pressure(element, elemDisc, 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