felocalresidual.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_FE_LOCAL_RESIDUAL_HELPER_HH
#define DUMUX_NAVIERSTOKES_MOMENTUM_CVFE_FE_LOCAL_RESIDUAL_HELPER_HH
 
#include <dune/geometry/quadraturerules.hh>
#include <dumux/common/typetraits/localdofs_.hh>
#include <dumux/common/boundaryflag.hh>
#include <dumux/discretization/extrusion.hh>
#include <dumux/discretization/fem/interpolationpointdata.hh>
 
#include "flux.hh"
 
namespace Dumux {
 
template<class Scalar, class NumEqVector, class LocalBasis, class Extrusion>
class NavierStokesMomentumFELocalResidual
{
    using RangeType = typename LocalBasis::Traits::RangeType;
 
public:
    template<class ResidualVector, class Problem, class FVElementGeometry, class ElementVariables>
    static void addStorageTerms(ResidualVector& residual,
                                const Problem& problem,
                                const FVElementGeometry& fvGeometry,
                                const ElementVariables& prevElemVars,
                                const ElementVariables& curElemVars,
                                const Scalar timeStepSize)
    {
        if constexpr (Detail::LocalDofs::hasNonCVLocalDofsInterface<FVElementGeometry>())
        {
            // Make sure we don't iterate over quadrature points if there are no hybrid dofs
            if (nonCVLocalDofs(fvGeometry).empty())
                return;
 
            const auto& localBasis = fvGeometry.feLocalBasis();
            std::vector<RangeType> integralShapeFunctions(localBasis.size(), RangeType(0.0));
 
            // We apply mass lumping such that we only need to calculate the integral of basis functions
            const auto& geometry = fvGeometry.elementGeometry();
            const auto& element = fvGeometry.element();
            using GlobalPosition = typename FVElementGeometry::GridGeometry::GlobalCoordinate;
            using FeIpData = FEInterpolationPointData<GlobalPosition, LocalBasis>;
 
            for (const auto& qpData : CVFE::quadratureRule(fvGeometry, element))
            {
                const auto& ipData = qpData.ipData();
                // Obtain and store shape function values and gradients at the current quad point
                FeIpData feIpData(geometry, ipData.local(), ipData.global(), localBasis);
 
                // get density from the problem
                for (const auto& localDof : nonCVLocalDofs(fvGeometry))
                    integralShapeFunctions[localDof.index()] += qpData.weight() * feIpData.shapeValue(localDof.index());
            }
 
            for (const auto& localDof : nonCVLocalDofs(fvGeometry))
            {
                const auto localDofIdx = localDof.index();
                const auto& data = ipData(fvGeometry, localDof);
                const auto curDensity = problem.density(element, fvGeometry, data, false);
                const auto prevDensity = problem.density(element, fvGeometry, data, true);
                const auto curVelocity = curElemVars[localDofIdx].velocity();
                const auto prevVelocity = prevElemVars[localDofIdx].velocity();
                auto timeDeriv = (curDensity*curVelocity - prevDensity*prevVelocity);
                timeDeriv /= timeStepSize;
 
                // add storage to residual
                for (int eqIdx = 0; eqIdx < NumEqVector::dimension; ++eqIdx)
                    residual[localDofIdx][eqIdx] += integralShapeFunctions[localDofIdx]*timeDeriv[eqIdx];
            }
        }
    }
 
    template<class ResidualVector, class Problem, class FVElementGeometry,
             class ElementVariables, class ElementFluxVariablesCache, class ElementBoundaryTypes>
    [[deprecated("This function is deprecated and will be removed after release 3.11. "
                 "Use addFluxAndSourceTerms(residual, problem, fvGeometry, elemVars, elemBcTypes) instead.")]]
    static void addFluxAndSourceTerms(ResidualVector& residual,
                                      const Problem& problem,
                                      const FVElementGeometry& fvGeometry,
                                      const ElementVariables& elemVars,
                                      const ElementFluxVariablesCache& elemFluxVarsCache,
                                      const ElementBoundaryTypes& elemBcTypes)
    {
        if constexpr (Detail::LocalDofs::hasNonCVLocalDofsInterface<FVElementGeometry>())
        {
            // Make sure we don't iterate over quadrature points if there are no hybrid dofs
            if (nonCVLocalDofs(fvGeometry).empty())
                return;
 
            if (!problem.pointSourceMap().empty())
                DUNE_THROW(Dune::NotImplemented, "Point sources are not implemented for hybrid momentum schemes.");
 
            static const bool enableUnsymmetrizedVelocityGradient
                = getParamFromGroup<bool>(problem.paramGroup(), "FreeFlow.EnableUnsymmetrizedVelocityGradient", false);
 
            const auto& element = fvGeometry.element();
            using FluxVariablesCache = typename ElementFluxVariablesCache::FluxVariablesCache;
            using FluxFunctionContext = NavierStokesMomentumFluxFunctionContext<Problem, FVElementGeometry, ElementVariables, FluxVariablesCache>;
            for (const auto& qpData : CVFE::quadratureRule(fvGeometry, element))
            {
                const auto& ipData = qpData.ipData();
                // Obtain and store shape function values and gradients at the current quad point
                const auto& cache = elemFluxVarsCache[ipData];
                FluxFunctionContext context(problem, fvGeometry, elemVars, cache);
                const auto& v = context.velocity();
                const auto& gradV = context.gradVelocity();
 
                // get viscosity from the problem
                const Scalar mu = problem.effectiveViscosity(element, fvGeometry, ipData);
                // get density from the problem
                const Scalar density = problem.density(element, fvGeometry, ipData);
 
                for (const auto& localDof : nonCVLocalDofs(fvGeometry))
                {
                    const auto localDofIdx = localDof.index();
                    NumEqVector fluxAndSourceTerm(0.0);
                    // add advection term
                    if (problem.enableInertiaTerms())
                        fluxAndSourceTerm -= density*(v*cache.gradN(localDofIdx))*v;
 
                    // add diffusion term
                    fluxAndSourceTerm += enableUnsymmetrizedVelocityGradient ?
                                            mu*mv(gradV, cache.gradN(localDofIdx))
                                            : mu*mv(gradV + getTransposed(gradV), cache.gradN(localDofIdx));
 
                    // add pressure term
                    fluxAndSourceTerm -= problem.pressure(element, fvGeometry, ipData) * cache.gradN(localDofIdx);
 
                    // finally add source and flux boundary terms and add everything to residual
                    auto sourceAtIp = problem.source(fvGeometry, elemVars, ipData);
                    // add gravity term rho*g (note that gravity might be zero in case it's disabled in the problem)
                    sourceAtIp += density * problem.gravity();
 
                    const auto& shapeValues = cache.shapeValues();
                    for (int eqIdx = 0; eqIdx < NumEqVector::dimension; ++eqIdx)
                    {
                        fluxAndSourceTerm[eqIdx] -= shapeValues[localDofIdx] * sourceAtIp[eqIdx];
                        residual[localDofIdx][eqIdx] += qpData.weight()*fluxAndSourceTerm[eqIdx];
                    }
                }
            }
 
            if (elemBcTypes.hasFluxBoundary())
                addBoundaryFluxes(residual, problem, fvGeometry, elemVars, elemFluxVarsCache, elemBcTypes);
        }
    }
 
    template<class ResidualVector, class Problem, class FVElementGeometry,
             class ElementVariables, class ElementFluxVariablesCache, class ElementBoundaryTypes>
    [[deprecated("This function is deprecated and will be removed after release 3.11. "
                 "Use addBoundaryFluxes(residual, problem, fvGeometry, elemVars, elemBcTypes) instead.")]]
    static void addBoundaryFluxes(ResidualVector& residual,
                                  const Problem& problem,
                                  const FVElementGeometry& fvGeometry,
                                  const ElementVariables& elemVars,
                                  const ElementFluxVariablesCache& elemFluxVarsCache,
                                  const ElementBoundaryTypes& elemBcTypes)
    {
        ResidualVector flux(0.0);
 
        const auto& element = fvGeometry.element();
        for (const auto& boundaryFace : boundaryFaces(fvGeometry))
        {
            const auto& bcTypes = elemBcTypes.get(fvGeometry, boundaryFace);
            if (!bcTypes.hasFluxBoundary())
                continue;
 
            problem.addBoundaryFluxIntegrals(flux, fvGeometry, elemVars, elemFluxVarsCache, boundaryFace, bcTypes);
        }
        residual += flux;
    }
 
    template<class ResidualVector, class Problem, class FVElementGeometry,
             class ElementVariables, class ElementBoundaryTypes>
    static void addFluxAndSourceTerms(ResidualVector& residual,
                                      const Problem& problem,
                                      const FVElementGeometry& fvGeometry,
                                      const ElementVariables& elemVars,
                                      const ElementBoundaryTypes& elemBcTypes)
    {
        if constexpr (Detail::LocalDofs::hasNonCVLocalDofsInterface<FVElementGeometry>())
        {
            // Make sure we don't iterate over quadrature points if there are no hybrid dofs
            if (nonCVLocalDofs(fvGeometry).empty())
                return;
 
            if (!problem.pointSourceMap().empty())
                DUNE_THROW(Dune::NotImplemented, "Point sources are not implemented for hybrid momentum schemes.");
 
            static const bool enableUnsymmetrizedVelocityGradient
                = getParamFromGroup<bool>(problem.paramGroup(), "FreeFlow.EnableUnsymmetrizedVelocityGradient", false);
 
            const auto& element = fvGeometry.element();
            using Cache = typename ElementVariables::InterpolationPointData;
            using FluxFunctionContext = NavierStokesMomentumFluxFunctionContext<Problem, FVElementGeometry, ElementVariables, Cache>;
            for (const auto& qpData : CVFE::quadratureRule(fvGeometry, element))
            {
                const auto& ipData = qpData.ipData();
                // Obtain and store shape function values and gradients at the current quad point
                const auto& ipCache = cache(elemVars, ipData);
                FluxFunctionContext context(problem, fvGeometry, elemVars, ipCache);
                const auto& v = context.velocity();
                const auto& gradV = context.gradVelocity();
 
                // get viscosity from the problem
                const Scalar mu = problem.effectiveViscosity(element, fvGeometry, ipData);
                // get density from the problem
                const Scalar density = problem.density(element, fvGeometry, ipData);
 
                for (const auto& localDof : nonCVLocalDofs(fvGeometry))
                {
                    const auto localDofIdx = localDof.index();
                    NumEqVector fluxAndSourceTerm(0.0);
                    // add advection term
                    if (problem.enableInertiaTerms())
                        fluxAndSourceTerm -= density*(v*ipCache.gradN(localDofIdx))*v;
 
                    // add diffusion term
                    fluxAndSourceTerm += enableUnsymmetrizedVelocityGradient ?
                                            mu*mv(gradV, ipCache.gradN(localDofIdx))
                                            : mu*mv(gradV + getTransposed(gradV), ipCache.gradN(localDofIdx));
 
                    // add pressure term
                    fluxAndSourceTerm -= problem.pressure(element, fvGeometry, ipData) * ipCache.gradN(localDofIdx);
 
                    // finally add source and flux boundary term and add everything to residual
                    auto sourceAtIp = problem.source(fvGeometry, elemVars, ipData);
                    // add gravity term rho*g (note that gravity might be zero in case it's disabled in the problem)
                    sourceAtIp += density * problem.gravity();
 
                    const auto& shapeValues = ipCache.shapeValues();
                    for (int eqIdx = 0; eqIdx < NumEqVector::dimension; ++eqIdx)
                    {
                        fluxAndSourceTerm[eqIdx] -= shapeValues[localDofIdx] * sourceAtIp[eqIdx];
                        residual[localDofIdx][eqIdx] += qpData.weight()*fluxAndSourceTerm[eqIdx];
                    }
                }
            }
 
            if (elemBcTypes.hasFluxBoundary())
                addBoundaryFluxes(residual, problem, fvGeometry, elemVars, elemBcTypes);
        }
    }
 
    template<class ResidualVector, class Problem, class FVElementGeometry,
             class ElementVariables, class ElementBoundaryTypes>
    static void addBoundaryFluxes(ResidualVector& residual,
                                  const Problem& problem,
                                  const FVElementGeometry& fvGeometry,
                                  const ElementVariables& elemVars,
                                  const ElementBoundaryTypes& elemBcTypes)
    {
        ResidualVector flux(0.0);
 
        const auto& element = fvGeometry.element();
        for (const auto& boundaryFace : boundaryFaces(fvGeometry))
        {
            const auto& bcTypes = elemBcTypes.get(fvGeometry, boundaryFace);
            if (!bcTypes.hasFluxBoundary())
                continue;
 
            problem.addBoundaryFluxIntegrals(flux, fvGeometry, elemVars, boundaryFace, bcTypes);
        }
        residual += flux;
    }
};
 
} // end namespace Dumux
 
#endif