New gradient limiter implementation and rewrite

src/dbns/limiter/MDLimiter.H

@@ -191,8 +191,6 @@ public:
             const GradFieldType& gradPhiIn = gradPhi.internalField();
 
             // Compute limiter values, internal faces
-            Type phiOwnerLimiter, phiNeighbourLimiter;
-
             forAll (owner, faceI)
             {
                 const label& own = owner[faceI];
@@ -201,10 +199,12 @@ public:
                 vector deltaRLeft = faceCentre[faceI] - cellCentre[own];
                 vector deltaRRight = faceCentre[faceI] - cellCentre[nei];
 
+                // Find minimal limiter value in each cell
+
                 // Owner side
                 limitFunction.limiter
                 (
-                    phiOwnerLimiter,
+                    phiLimiterIn[own],
                     cellVolume[own],
                     phiMaxIn[own] - phi_[own],
                     phiMinIn[own] - phi_[own],
@@ -214,19 +214,12 @@ public:
                 // Neighbour side
                 limitFunction.limiter
                 (
-                    phiNeighbourLimiter,
+                    phiLimiterIn[nei],
                     cellVolume[nei],
                     phiMaxIn[nei] - phi_[nei],
                     phiMinIn[nei] - phi_[nei],
                     (deltaRRight & gradPhiIn[nei])                    
                 );
-
-                // Find minimal limiter value in each cell
-                phiLimiterIn[own] =
-                    min(phiLimiterIn[own], phiOwnerLimiter);
-
-                phiLimiterIn[nei] =
-                    min(phiLimiterIn[nei], phiNeighbourLimiter);
             }
 
             // Coupled boundaries
@@ -251,24 +244,19 @@ public:
                     const GradFieldType gradPhiRight =
                         gradPhi.boundaryField()[patchI].patchNeighbourField();
 
-                    Type phiFCLimiter;
+                    // Find minimal limiter value in each cell
-
                     forAll (fc, faceI)
                     {
                         const label& curFC = fc[faceI];
 
                         limitFunction.limiter
                         (
-                            phiFCLimiter,
+                            phiLimiterIn[curFC],
                             cellVolume[curFC],
                             phiMaxIn[curFC] - phi_[curFC],
                             phiMinIn[curFC] - phi_[curFC],
                             (deltaR[faceI] & gradPhiLeft[faceI])
                         );
-
-                        // Find minimal limiter value in each cell
-                        phiLimiterIn[curFC] =
-                            min(phiLimiterIn[curFC], phiFCLimiter);
                     }
                 }
             }

src/finiteVolume/Make/files

@@ -358,6 +358,10 @@ $(limitedGradSchemes)/faceLimitedGrad/faceLimitedGrads.C
 $(limitedGradSchemes)/cellLimitedGrad/cellLimitedGrads.C
 $(limitedGradSchemes)/faceMDLimitedGrad/faceMDLimitedGrads.C
 $(limitedGradSchemes)/cellMDLimitedGrad/cellMDLimitedGrads.C
+$(limitedGradSchemes)/barthJespersenGrad/barthJespersenGrads.C
+$(limitedGradSchemes)/venkatakrishnanGrad/venkatakrishnanGrads.C
+$(limitedGradSchemes)/wangGrad/wangGrads.C
+$(limitedGradSchemes)/michalakGoochGrad/michalakGoochGrads.C
 
 snGradSchemes = finiteVolume/snGradSchemes
 $(snGradSchemes)/snGradScheme/snGradSchemes.C

src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrads.C

@@ -66,6 +66,7 @@ inline void cellLimitedGrad<scalar>::limitFace
     }
 }
 
+
 template<class Type>
 inline void cellLimitedGrad<Type>::limitFace
 (
@@ -229,7 +230,7 @@ tmp<volVectorField> cellLimitedGrad<scalar>::calcGrad
 
     if (fv::debug)
     {
-        Info<< "gradient limiter for: " << vsf.name()
+        Info<< "Explicit scalar gradient limiter for: " << vsf.name()
             << " max = " << gMax(limiter)
             << " min = " << gMin(limiter)
             << " average: " << gAverage(limiter) << endl;
@@ -381,7 +382,7 @@ tmp<volTensorField> cellLimitedGrad<vector>::calcGrad
 
     if (fv::debug)
     {
-        Info<< "gradient limiter for: " << vsf.name()
+        Info<< "Explicit vector gradient limiter for: " << vsf.name()
             << " max = " << gMax(limiter)
             << " min = " << gMin(limiter)
             << " average: " << gAverage(limiter) << endl;
@@ -563,7 +564,7 @@ tmp<BlockLduSystem<vector, vector> > cellLimitedGrad<scalar>::fvmGrad
 
     if (fv::debug)
     {
-        Info<< "gradient limiter for: " << vsf.name()
+        Info<< "Implicit scalar gradient limiter for: " << vsf.name()
             << " max = " << gMax(lfIn)
             << " min = " << gMin(lfIn)
             << " average: " << gAverage(lfIn) << endl;

src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellMDLimitedGrad/cellMDLimitedGrad.H

@@ -141,7 +141,11 @@ inline void cellMDLimitedGrad<scalar>::limitFace
 {
     scalar extrapolate = dcf & g;
 
-    if (extrapolate > maxDelta)
+    if (mag(extrapolate) < SMALL)
+    {
+        // Limiter remains unchanged
+    }
+    else if (extrapolate > maxDelta)
     {
         g = g + dcf*(maxDelta - extrapolate)/magSqr(dcf);
     }

src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/limitedGrad/LimitedGrad.C

@@ -67,16 +67,13 @@ LimitedGrad<Type, GradientLimiter>::limiter
         label own = owner[faceI];
         label nei = neighbour[faceI];
 
-        const Type& vfOwn = vfIn[own];
-        const Type& vfNei = vfIn[nei];
-
         // Owner side
-        maxVf[own] = max(maxVf[own], vfNei);
+        maxVf[own] = Foam::max(maxVf[own], vfIn[nei]);
-        minVf[own] = min(minVf[own], vfNei);
+        minVf[own] = Foam::min(minVf[own], vfIn[nei]);
 
         // Neighbour side
-        maxVf[nei] = max(maxVf[nei], vfOwn);
+        maxVf[nei] = Foam::max(maxVf[nei], vfIn[own]);
-        minVf[nei] = min(minVf[nei], vfOwn);
+        minVf[nei] = Foam::min(minVf[nei], vfIn[own]);
     }
 
     const GeoBoundaryFieldType& bf = vf.boundaryField();
@@ -96,10 +93,9 @@ LimitedGrad<Type, GradientLimiter>::limiter
             forAll (pOwner, pFaceI)
             {
                 label own = pOwner[pFaceI];
-                Type vfNei = psfNei[pFaceI];
 
-                maxVf[own] = max(maxVf[own], vfNei);
+                maxVf[own] = Foam::max(maxVf[own], psfNei[pFaceI]);
-                minVf[own] = min(minVf[own], vfNei);
+                minVf[own] = Foam::min(minVf[own], psfNei[pFaceI]);
             }
         }
         else
@@ -108,18 +104,20 @@ LimitedGrad<Type, GradientLimiter>::limiter
             forAll (pOwner, pFaceI)
             {
                 label own = pOwner[pFaceI];
-                Type vfNei = psf[pFaceI];
+                const Type& vfNei = psf[pFaceI];
 
-                maxVf[own] = max(maxVf[own], vfNei);
+                maxVf[own] = Foam::max(maxVf[own], vfNei);
-                minVf[own] = min(minVf[own], vfNei);
+                minVf[own] = Foam::min(minVf[own], vfNei);
             }
         }
     }
 
     // Subtract the cell value to get differences
-    maxVf -= vf;
+    // Stabilise differences for round-off error, since maxVf must be
-    minVf -= vf;
+    // positive or zero and minVf must be negative or zero
-
+    // HJ, 1/Feb/2016
+    maxVf = Foam::max(maxVf - vf, pTraits<Type>::zero);
+    minVf = Foam::min(minVf - vf, pTraits<Type>::zero);
 
     // Create a limiter
     tmp<GeoFieldType> tlimitField
@@ -135,14 +133,18 @@ LimitedGrad<Type, GradientLimiter>::limiter
                 IOobject::NO_WRITE
             ),
             mesh,
-            dimensioned<Type>("zero", dimless, pTraits<Type>::zero),
+            dimensioned<Type>("zero", dimless, pTraits<Type>::one),
             zeroGradientFvPatchField<Type>::typeName
         )
     );
     GeoFieldType& limitField = tlimitField();
 
     const volVectorField& C = mesh.C();
+    const vectorField& CIn = C.internalField();
+
     const surfaceVectorField& Cf = mesh.Cf();
+    const vectorField& CfIn = Cf.internalField();
+
     const scalarField& vols = mesh.V().field();
 
     Field<Type>& lfIn = limitField.internalField();
@@ -164,7 +166,7 @@ LimitedGrad<Type, GradientLimiter>::limiter
             vols[own],
             maxVf[own],
             minVf[own],
-            (Cf[faceI] - C[own]) & g[own]
+            (CfIn[faceI] - CIn[own]) & g[own]
         );
 
         // Neighbour side
@@ -174,7 +176,7 @@ LimitedGrad<Type, GradientLimiter>::limiter
             vols[nei],
             maxVf[nei],
             minVf[nei],
-            (Cf[faceI] - C[nei]) & g[nei]
+            (CfIn[faceI] - CIn[nei]) & g[nei]
         );
     }
 
@@ -199,7 +201,7 @@ LimitedGrad<Type, GradientLimiter>::limiter
     }
 
     // Update coupled boundaries for patchNeighbourField
-    limitField.boundaryField().evaluateCoupled();
+    limitField.correctBoundaryConditions();
 
     if (fv::debug)
     {
@@ -234,7 +236,7 @@ LimitedGrad<Type, GradientLimiter>::gradientField
     GeoGradFieldType& gradVf = tGrad();
 
     // Apply the limiter
-    GeoFieldType limitField = this->limiter(vf, gradVf);
+    GeoFieldType limitField(this->limiter(vf, gradVf));
 
     gradVf.internalField() =
         scaleRow(gradVf.internalField(), limitField.internalField());
@@ -261,11 +263,15 @@ LimitedGrad<Type, GradientLimiter>::gradientMatrix
     GradMatrixType& bs = tbs();
 
     // Calculate limiter.  Using explicit gradient
-    GeoFieldType limitField = this->limiter
+    GeoFieldType limitField
+    (
+        this->limiter
         (
             vf,
             basicGradScheme_().grad(vf)()
+        )
     );
+
     const Field<Type>& lfIn = limitField.internalField();
 
     typedef typename CoeffField<vector>::linearTypeField
@@ -288,14 +294,14 @@ LimitedGrad<Type, GradientLimiter>::gradientMatrix
     const unallocLabelList& owner = mesh.owner();
     const unallocLabelList& neighbour = mesh.neighbour();
 
-    forAll(u, faceI)
+    forAll (u, faceI)
     {
         u[faceI] = scaleRow(u[faceI], lfIn[owner[faceI]]);
         l[faceI] = scaleRow(l[faceI], lfIn[neighbour[faceI]]);
     }
 
     // Limit diag and source coeffs
-    forAll(d, cellI)
+    forAll (d, cellI)
     {
         d[cellI] = scaleRow(d[cellI], lfIn[cellI]);
 
@@ -303,7 +309,7 @@ LimitedGrad<Type, GradientLimiter>::gradientMatrix
     }
 
     // Limit coupling coeffs
-    forAll(vf.boundaryField(), patchI)
+    forAll (vf.boundaryField(), patchI)
     {
         const fvPatchField<Type>& pf = vf.boundaryField()[patchI];
         const fvPatch& patch = pf.patch();
@@ -321,7 +327,7 @@ LimitedGrad<Type, GradientLimiter>::gradientMatrix
             const Field<Type> lfNei =
                 limitField.boundaryField()[patchI].patchNeighbourField();
 
-            forAll(pf, faceI)
+            forAll (pf, faceI)
             {
                 pcoupleUpper[faceI] =
                     scaleRow(pcoupleUpper[faceI], lfIn[fc[faceI]]);

src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/minModGrad/minModGrad.H

@@ -1,139 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
-   \\    /   O peration     | Version:     3.2
-    \\  /    A nd           | Web:         http://www.foam-extend.org
-     \\/     M anipulation  | For copyright notice see file Copyright
--------------------------------------------------------------------------------
-License
-    This file is part of foam-extend.
-
-    foam-extend 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.
-
-    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::fv::minModGrad
-
-Description
-    MinMod gradient limiter applied to a runTime selected
-     base gradient scheme.
-
-SourceFiles
-    minModGrad.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef minModGrad_H
-#define minModGrad_H
-
-#include "gradScheme.H"
-#include "LimitedGrad.H"
-#include "MinModLimiter.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace fv
-{
-
-/*---------------------------------------------------------------------------*\
-                     Class minModGrad Declaration
-\*---------------------------------------------------------------------------*/
-
-template<class Type>
-class minModGrad
-:
-    public fv::gradScheme<Type>,
-    public LimitedGrad<Type, MinModLimiter>
-{
-    // Private Member Functions
-
-        //- Disallow default bitwise copy construct
-        minModGrad(const minModGrad&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const minModGrad&);
-
-
-public:
-
-    //- RunTime type information
-    TypeName("minMod");
-
-
-    // Constructors
-
-        //- Construct from mesh and schemeData
-        minModGrad(const fvMesh& mesh, Istream& schemeData)
-        :
-            gradScheme<Type>(mesh),
-            LimitedGrad<Type, MinModLimiter>(mesh, schemeData)
-        {}
-
-
-    // Member Functions
-
-        //- Return the gradient of the given field to the gradScheme::grad
-        //  for optional caching
-        virtual tmp
-        <
-            GeometricField
-            <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>
-        > calcGrad
-        (
-            const GeometricField<Type, fvPatchField, volMesh>& vf,
-            const word& name
-        ) const
-        {
-            return LimitedGrad<Type, MinModLimiter>::gradientField
-            (
-                vf,
-                name
-            );
-        }
-
-        //- Return the BlockLduSystem corresponding to the implicit cell
-        //  limited grad discretization.  For block coupled systems.
-        virtual tmp
-        <
-            BlockLduSystem<vector, typename outerProduct<vector, Type>::type>
-        > fvmGrad
-        (
-            const GeometricField<Type, fvPatchField, volMesh>& vf
-        ) const
-        {
-            return LimitedGrad<Type, MinModLimiter>::gradientMatrix
-            (
-                vf
-            );
-        }
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace fv
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/minModGrad/minModGrads.C

@@ -1,89 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
-   \\    /   O peration     | Version:     3.2
-    \\  /    A nd           | Web:         http://www.foam-extend.org
-     \\/     M anipulation  | For copyright notice see file Copyright
--------------------------------------------------------------------------------
-License
-    This file is part of foam-extend.
-
-    foam-extend 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.
-
-    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "minModGrad.H"
-#include "fvMesh.H"
-#include "volMesh.H"
-#include "surfaceMesh.H"
-#include "volFields.H"
-#include "fixedValueFvPatchFields.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace fv
-{
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-makeFvGradScheme(minModGrad)
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-
-template<>
-tmp
-<
-    BlockLduSystem<vector, outerProduct<vector, vector>::type>
->
-minModGrad<vector>::fvmGrad
-(
-    const volVectorField& vf
-) const
-{
-    FatalErrorIn
-    (
-        "tmp<BlockLduSystem> minModGrad<vector>::fvmGrad\n"
-        "(\n"
-        "    GeometricField<vector, fvPatchField, volMesh>&"
-        ")\n"
-    )   << "Implicit gradient operators with cell limiters defined only for "
-        << "scalar."
-        << abort(FatalError);
-
-    typedef outerProduct<vector, vector>::type GradType;
-
-    tmp<BlockLduSystem<vector, GradType> > tbs
-    (
-        new BlockLduSystem<vector, GradType>(vf.mesh())
-    );
-
-    return tbs;
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace fv
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// ************************************************************************* //

src/finiteVolume/finiteVolume/gradientLimiters/BarthJespersenLimiter.H

@@ -28,15 +28,14 @@ Description
     Barth-Jespersen limiter
 
 Author
-    Aleksandar Jemcov, Hrvoje Jasak
+    Aleksandar Jemcov
+    Rewrite by Hrvoje Jasak
 
 \*---------------------------------------------------------------------------*/
 
 #ifndef BarthJespersenLimiter_H
 #define BarthJespersenLimiter_H
 
-#include "vector.H"
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -74,19 +73,16 @@ public:
         {
             if (mag(extrapolate) < SMALL)
             {
-                lim = 1;
+                return;
             }
-            else
+
+            if (extrapolate - deltaOneMax > SMALL)
             {
-                lim = min
+                lim = min(lim, deltaOneMax/extrapolate);
-                (
+            }
-                    max
+            else if (extrapolate - deltaOneMin < -SMALL)
-                    (
+            {
-                        max(deltaOneMax/extrapolate, 0),
+                lim = min(lim, deltaOneMin/extrapolate);
-                        max(deltaOneMin/extrapolate, 0)
-                    ),
-                    1
-                );
             }
         }
 

src/finiteVolume/finiteVolume/gradientLimiters/MinModLimiter.H

@@ -1,126 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
-   \\    /   O peration     | Version:     3.2
-    \\  /    A nd           | Web:         http://www.foam-extend.org
-     \\/     M anipulation  | For copyright notice see file Copyright
--------------------------------------------------------------------------------
-License
-    This file is part of foam-extend.
-
-    foam-extend 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.
-
-    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    MinModLimiter
-
-Description
-    Non-differentiable cell limiter.  The limiter limits the gradient such
-    that the extrapolated face value does not under- or over-shoot the
-    neighbouring cell values
-
-Author
-    Hrvoje Jasak
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef MinModLimiter_H
-#define MinModLimiter_H
-
-#include "vector.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                         Class MinModLimiter Declaration
-\*---------------------------------------------------------------------------*/
-
-class MinModLimiter
-{
-public:
-
-    // Constructors
-
-        //- Construct null
-        MinModLimiter()
-        {}
-
-
-    // Destructor - default
-
-
-    // Member functions
-
-        //- Set scalar limiter value
-        inline void limiter
-        (
-            scalar& lim,
-            const scalar& cellVolume,
-            const scalar& deltaOneMax,
-            const scalar& deltaOneMin,
-            const scalar& extrapolate
-        )
-        {
-            if (mag(extrapolate) < SMALL)
-            {
-                lim = 1;
-            }
-            else
-            {
-                if (extrapolate > deltaOneMax + VSMALL)
-                {
-                    lim = min(1, deltaOneMax/extrapolate);
-                }
-                else if (extrapolate < deltaOneMin - VSMALL)
-                {
-                    lim = min(1, deltaOneMin/extrapolate);
-                }
-            }
-        }
-
-        //- Set Type limiter
-        template<class Type>
-        inline void limiter
-        (
-            Type& lim,
-            const scalar& cellVolume,
-            const Type& deltaOneMax,
-            const Type& deltaOneMin,
-            const Type& extrapolate
-        )
-        {
-            for (direction cmpt = 0; cmpt < Type::nComponents; cmpt++)
-            {
-                limiter
-                (
-                    lim.component(cmpt),
-                    cellVolume,
-                    deltaOneMax.component(cmpt),
-                    deltaOneMin.component(cmpt),
-                    extrapolate.component(cmpt)
-                );
-            }
-        }
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-#endif
-
-// ************************************************************************* //

src/finiteVolume/finiteVolume/gradientLimiters/VenkatakrishnanLimiter.H

@@ -29,15 +29,13 @@ Description
 
 Author
     Aleksandar Jemcov
-    Updated by Hrvoje Jasak
+    Rewrite by Hrvoje Jasak
 
 \*---------------------------------------------------------------------------*/
 
 #ifndef VenkatakrishnanLimiter_H
 #define VenkatakrishnanLimiter_H
 
-#include "vector.H"
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -77,10 +75,16 @@ public:
         {
             scalar epsilonSquare = pow3(k_)*cellVolume;
 
-            if (extrapolate > SMALL)
+            if (mag(extrapolate) < SMALL)
+            {
+                // Limiter remains unchanged
+                return;
+            }
+            else if (extrapolate > VSMALL)
             {
                 lim = max
                 (
+                    0,
                     min
                     (
                         (
@@ -92,21 +96,21 @@ public:
                             extrapolate*
                             (
                                 sqr(deltaOneMax)
-                              + 2.0*sqr(extrapolate)
+                              + 2*sqr(extrapolate)
                               + deltaOneMax*extrapolate
                               + epsilonSquare
                             ),
                             SMALL
                         ),
-                        1
+                        lim
-                    ),
+                    )
-                    0
                 );
             }
-            else if (extrapolate < SMALL)
+            else if (extrapolate < VSMALL)
             {
                 lim = max
                 (
+                    0,
                     min
                     (
                         (
@@ -118,20 +122,16 @@ public:
                             extrapolate*
                             (
                                 sqr(deltaOneMin)
-                              + 2.0*sqr(extrapolate) + deltaOneMin*extrapolate
+                              + 2*sqr(extrapolate)
+                              + deltaOneMin*extrapolate
                               + epsilonSquare
                             ),
                             SMALL
                         ),
-                        1
+                        lim
-                    ),
+                    )
-                    0
                 );
             }
-            else
-            {
-                lim = 1;
-            }
         }
 
         //- Set Type limiter

src/finiteVolume/finiteVolume/gradientLimiters/WangLimiter.H

@@ -28,16 +28,13 @@ Description
     Wang differentiable limiter
 
 Author
-    Aleksandar Jemcov
+    Hrvoje Jasak
-    Updated by Hrvoje Jasak
 
 \*---------------------------------------------------------------------------*/
 
 #ifndef WangLimiter_H
 #define WangLimiter_H
 
-#include "vector.H"
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -78,10 +75,11 @@ public:
             scalar epsilonSquare =
                 sqr(epsilonPrime_*(deltaOneMax - deltaOneMin));
 
-            if (extrapolate > SMALL)
+            if (extrapolate > VSMALL)
             {
                 lim = max
                 (
+                    0,
                     min
                     (
                         (
@@ -93,21 +91,21 @@ public:
                             extrapolate*
                             (
                                 sqr(deltaOneMax)
-                              + 2.0*sqr(extrapolate)
+                              + 2*sqr(extrapolate)
                               + deltaOneMax*extrapolate
                               + epsilonSquare
                             ),
                             SMALL
                         ),
-                        1
+                        lim
-                    ),
+                    )
-                    0
                 );
             }
-            else if (extrapolate < SMALL)
+            else if (extrapolate < VSMALL)
             {
                 lim = max
                 (
+                    0,
                     min
                     (
                         (
@@ -119,19 +117,19 @@ public:
                             extrapolate*
                             (
                                 sqr(deltaOneMin)
-                              + 2.0*sqr(extrapolate) + deltaOneMin*extrapolate
+                              + 2*sqr(extrapolate)
+                              + deltaOneMin*extrapolate
                               + epsilonSquare
                             ),
                             SMALL
                         ),
-                        1
+                        lim
-                    ),
+                    )
-                    0
                 );
             }
             else
             {
-                lim = 1;
+                // Limiter remains unchanged
             }
         }