[openfoammaofnepo]

Dear All,

In the openfoam limiter function, there is a scalar called k_, does anybody know what does it stand for? For example, we can have a look at the limiter for Gamma differencing:

    GammaLimiter(Istream& is)
    :
        k_(readScalar(is))
    {
        if (k_ < 0 || k_ > 1)
        {
            FatalIOErrorIn("GammaLimiter(Istream& is)", is)
                << "coefficient = " << k_
                << " should be >= 0 and <= 1"
                << exit(FatalIOError);
        }

        // Rescale k_ to be >= 0 and <= 0.5 (TVD conformant)
        // and avoid the /0 when k_ = 0
        k_ = max(k_/2.0, SMALL);
    }

The above code lines are from:

OpenFOAM/OpenFOAM-<version>/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/Gamma/Gamma.H

The second question is: we need the calculate the gradiant ratio for the limiter function (always the input for the limiter). In TVD scheme, we call it r and in NVD we call it phict. These two quantities are defined in:

OpenFOAM/OpenFOAM-<version>/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/NVDTVD.H

When calculating r or phict, we need to get the difference of phi in Cell N and P (actually upwind and downwind cells). Take r for example. The original code is:

        scalar r
        (
            const scalar faceFlux,
            const scalar phiP,
            const scalar phiN,
            const vector& gradcP,
            const vector& gradcN,
            const vector& d
        ) const
        {
            scalar gradf = phiN - phiP;

            scalar gradcf;

            if (faceFlux > 0)
            {
                gradcf = d & gradcP;
            }
            else
            {
                gradcf = d & gradcN;
            }

            if (mag(gradcf) >= 1000*mag(gradf))
            {
                return 2*1000*sign(gradcf)*sign(gradf) - 1;
            }
            else
            {
                return 2*(gradcf/gradf) - 1;
            }
        }

In my opinion, it should be like:

        scalar r
        (
            const scalar faceFlux,
            const scalar phiP,
            const scalar phiN,
            const vector& gradcP,
            const vector& gradcN,
            const vector& d
        ) const
        {
            

            scalar gradcf;

            if (faceFlux > 0)
            {
                gradcf = d & gradcP;
                scalar gradf = phiN - phiP;
            }
            else
            {
                gradcf = d & gradcN;
                scalar gradf = phiP - phiN;
            }

            if (mag(gradcf) >= 1000*mag(gradf))
            {
                return 2*1000*sign(gradcf)*sign(gradf) - 1;
            }
            else
            {
                return 2*(gradcf/gradf) - 1;
            }
        }

About this issue, I am not sure this is a bug or I make some mistakes here. In my opinion, gradf depends on the direction of flow, i.e. upwind or downwind, but as you know owner and neighbor cells only stands for the geometrical position, instead of any flow direction information. Does anybody have some comments here?

OFFO