flux/upwindscheme.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_DISCRETIZATION_UPWINDSCHEME_HH
#define DUMUX_DISCRETIZATION_UPWINDSCHEME_HH
 
#include <dumux/common/parameters.hh>
#include <dumux/discretization/method.hh>
 
namespace Dumux {
 
template<class GridGeometry, class DiscretizationMethod>
class UpwindSchemeImpl;
 
template<class GridGeometry>
using UpwindScheme = UpwindSchemeImpl<GridGeometry, typename GridGeometry::DiscretizationMethod>;
 
namespace Detail {
 
template<class ElemVolVars, class SubControlVolumeFace, class UpwindTermFunction, class Scalar>
Scalar upwindSchemeMultiplier(const ElemVolVars& elemVolVars,
                              const SubControlVolumeFace& scvf,
                              const UpwindTermFunction& upwindTerm,
                              Scalar flux, int phaseIdx)
{
    // TODO: pass this from outside?
    static const Scalar upwindWeight = getParamFromGroup<Scalar>(elemVolVars.gridVolVars().problem().paramGroup(), "Flux.UpwindWeight");
 
    const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
    const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
 
    using std::signbit;
    if (signbit(flux)) // if sign of flux is negative
        return (upwindWeight*upwindTerm(outsideVolVars)
                + (1.0 - upwindWeight)*upwindTerm(insideVolVars));
    else
        return (upwindWeight*upwindTerm(insideVolVars)
                + (1.0 - upwindWeight)*upwindTerm(outsideVolVars));
}
 
} // end namespace Detail
 
template<class GridGeometry>
class UpwindSchemeImpl<GridGeometry, DiscretizationMethods::Box>
{
public:
    template<class FluxVariables, class UpwindTermFunction, class Scalar>
    static Scalar apply(const FluxVariables& fluxVars,
                        const UpwindTermFunction& upwindTerm,
                        Scalar flux, int phaseIdx)
    {
        return apply(fluxVars.elemVolVars(), fluxVars.scvFace(), upwindTerm, flux, phaseIdx);
    }
 
    template<class ElemVolVars, class SubControlVolumeFace, class UpwindTermFunction, class Scalar>
    static Scalar apply(const ElemVolVars& elemVolVars,
                        const SubControlVolumeFace& scvf,
                        const UpwindTermFunction& upwindTerm,
                        Scalar flux, int phaseIdx)
    {
        return flux*multiplier(elemVolVars, scvf, upwindTerm, flux, phaseIdx);
    }
 
    template<class ElemVolVars, class SubControlVolumeFace, class UpwindTermFunction, class Scalar>
    static Scalar multiplier(const ElemVolVars& elemVolVars,
                             const SubControlVolumeFace& scvf,
                             const UpwindTermFunction& upwindTerm,
                             Scalar flux, int phaseIdx)
    {
        return Detail::upwindSchemeMultiplier(elemVolVars, scvf, upwindTerm, flux, phaseIdx);
    }
};
 
template<class GridGeometry>
class UpwindSchemeImpl<GridGeometry, DiscretizationMethods::CCTpfa>
{
    using GridView = typename GridGeometry::GridView;
    static constexpr int dim = GridView::dimension;
    static constexpr int dimWorld = GridView::dimensionworld;
 
public:
    template<class ElemVolVars, class SubControlVolumeFace, class UpwindTermFunction, class Scalar>
    static Scalar multiplier(const ElemVolVars& elemVolVars,
                             const SubControlVolumeFace& scvf,
                             const UpwindTermFunction& upwindTerm,
                             Scalar flux, int phaseIdx)
    {
        static_assert(dim == dimWorld, "Multiplier cannot be computed on surface grids using cell-centered scheme!");
        return Detail::upwindSchemeMultiplier(elemVolVars, scvf, upwindTerm, flux, phaseIdx);
    }
 
    template<class ElemVolVars, class SubControlVolumeFace, class UpwindTermFunction, class Scalar>
    static Scalar apply(const ElemVolVars& elemVolVars,
                        const SubControlVolumeFace& scvf,
                        const UpwindTermFunction& upwindTerm,
                        Scalar flux, int phaseIdx)
    {
        static_assert(dim == dimWorld, "This upwind scheme cannot be used on surface grids using cell-centered schemes!");
        return flux*multiplier(elemVolVars, scvf, upwindTerm, flux, phaseIdx);
    }
 
    // For surface and network grids (dim < dimWorld) we have to do a special upwind scheme
    template<class FluxVariables, class UpwindTermFunction, class Scalar>
    static Scalar apply(const FluxVariables& fluxVars,
                        const UpwindTermFunction& upwindTerm,
                        Scalar flux, int phaseIdx)
    {
        const auto& scvf = fluxVars.scvFace();
        const auto& elemVolVars = fluxVars.elemVolVars();
 
        // standard scheme
        if constexpr (dim == dimWorld)
            return apply(elemVolVars, scvf, upwindTerm, flux, phaseIdx);
 
        // on non-branching points the standard scheme works
        if (scvf.numOutsideScvs() == 1)
            return flux*Detail::upwindSchemeMultiplier(elemVolVars, scvf, upwindTerm, flux, phaseIdx);
 
        // handle branching points with a more complicated upwind scheme
        // we compute a flux-weighted average of all inflowing branches
        const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
 
        Scalar branchingPointUpwindTerm = 0.0;
        Scalar sumUpwindFluxes = 0.0;
 
        // if the inside flux is positive (outflow) do fully upwind and return flux
        using std::signbit;
        if (!signbit(flux))
            return upwindTerm(insideVolVars)*flux;
        else
            sumUpwindFluxes += flux;
 
        for (unsigned int i = 0; i < scvf.numOutsideScvs(); ++i)
        {
            // compute the outside flux
            const auto& fvGeometry = fluxVars.fvGeometry();
            const auto outsideScvIdx = scvf.outsideScvIdx(i);
            const auto outsideElement = fvGeometry.gridGeometry().element(outsideScvIdx);
            const auto& flippedScvf = fvGeometry.flipScvf(scvf.index(), i);
 
            using AdvectionType = typename FluxVariables::AdvectionType;
            const auto outsideFlux = AdvectionType::flux(fluxVars.problem(),
                                                         outsideElement,
                                                         fvGeometry,
                                                         elemVolVars,
                                                         flippedScvf,
                                                         phaseIdx,
                                                         fluxVars.elemFluxVarsCache());
 
            if (!signbit(outsideFlux))
                branchingPointUpwindTerm += upwindTerm(elemVolVars[outsideScvIdx])*outsideFlux;
            else
                sumUpwindFluxes += outsideFlux;
        }
 
        // the flux might be zero
        if (sumUpwindFluxes != 0.0)
            branchingPointUpwindTerm /= -sumUpwindFluxes;
        else
            branchingPointUpwindTerm = 0.0;
 
        // upwind scheme (always do fully upwind at branching points)
        // a weighting here would lead to an error since the derivation is based on a fully upwind scheme
        // TODO How to implement a weight of e.g. 0.5
        if (signbit(flux))
            return flux*branchingPointUpwindTerm;
        else
            return flux*upwindTerm(insideVolVars);
    }
};
 
template<class GridGeometry>
class UpwindSchemeImpl<GridGeometry, DiscretizationMethods::CCMpfa>
: public UpwindSchemeImpl<GridGeometry, DiscretizationMethods::CCTpfa> {};
 
} // end namespace Dumux
 
#endif