diff --git a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrad.C b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrad.C new file mode 100644 index 0000000000..bbd61fbad6 --- /dev/null +++ b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrad.C @@ -0,0 +1,230 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2018 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM 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. + + OpenFOAM 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 OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "cellLimitedGrad.H" +#include "gaussGrad.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template +void Foam::fv::cellLimitedGrad::limitGradient +( + const Field& limiter, + Field& gIf +) const +{ + gIf *= limiter; +} + + +template +void Foam::fv::cellLimitedGrad::limitGradient +( + const Field& limiter, + Field& gIf +) const +{ + forAll(gIf, celli) + { + gIf[celli] = tensor + ( + cmptMultiply(limiter[celli], gIf[celli].x()), + cmptMultiply(limiter[celli], gIf[celli].y()), + cmptMultiply(limiter[celli], gIf[celli].z()) + ); + } +} + + +template +Foam::tmp +< + Foam::GeometricField + < + typename Foam::outerProduct::type, + Foam::fvPatchField, + Foam::volMesh + > +> +Foam::fv::cellLimitedGrad::calcGrad +( + const GeometricField& vsf, + const word& name +) const +{ + const fvMesh& mesh = vsf.mesh(); + + tmp + < + GeometricField + ::type, fvPatchField, volMesh> + > tGrad = basicGradScheme_().calcGrad(vsf, name); + + if (k_ < SMALL) + { + return tGrad; + } + + GeometricField + < + typename outerProduct::type, + fvPatchField, + volMesh + >& g = tGrad.ref(); + + const labelUList& owner = mesh.owner(); + const labelUList& neighbour = mesh.neighbour(); + + const volVectorField& C = mesh.C(); + const surfaceVectorField& Cf = mesh.Cf(); + + Field maxVsf(vsf.primitiveField()); + Field minVsf(vsf.primitiveField()); + + forAll(owner, facei) + { + label own = owner[facei]; + label nei = neighbour[facei]; + + const Type& vsfOwn = vsf[own]; + const Type& vsfNei = vsf[nei]; + + maxVsf[own] = max(maxVsf[own], vsfNei); + minVsf[own] = min(minVsf[own], vsfNei); + + maxVsf[nei] = max(maxVsf[nei], vsfOwn); + minVsf[nei] = min(minVsf[nei], vsfOwn); + } + + + const typename GeometricField::Boundary& bsf = + vsf.boundaryField(); + + forAll(bsf, patchi) + { + const fvPatchField& psf = bsf[patchi]; + const labelUList& pOwner = mesh.boundary()[patchi].faceCells(); + + if (psf.coupled()) + { + const Field psfNei(psf.patchNeighbourField()); + + forAll(pOwner, pFacei) + { + label own = pOwner[pFacei]; + const Type& vsfNei = psfNei[pFacei]; + + maxVsf[own] = max(maxVsf[own], vsfNei); + minVsf[own] = min(minVsf[own], vsfNei); + } + } + else + { + forAll(pOwner, pFacei) + { + label own = pOwner[pFacei]; + const Type& vsfNei = psf[pFacei]; + + maxVsf[own] = max(maxVsf[own], vsfNei); + minVsf[own] = min(minVsf[own], vsfNei); + } + } + } + + maxVsf -= vsf; + minVsf -= vsf; + + if (k_ < 1.0) + { + const Field maxMinVsf((1.0/k_ - 1.0)*(maxVsf - minVsf)); + maxVsf += maxMinVsf; + minVsf -= maxMinVsf; + } + + + // Create limiter initialized to 1 + // Note: the limiter is not permitted to be > 1 + Field limiter(vsf.primitiveField().size(), pTraits::one); + + forAll(owner, facei) + { + label own = owner[facei]; + label nei = neighbour[facei]; + + // owner side + limitFace + ( + limiter[own], + maxVsf[own], + minVsf[own], + (Cf[facei] - C[own]) & g[own] + ); + + // neighbour side + limitFace + ( + limiter[nei], + maxVsf[nei], + minVsf[nei], + (Cf[facei] - C[nei]) & g[nei] + ); + } + + forAll(bsf, patchi) + { + const labelUList& pOwner = mesh.boundary()[patchi].faceCells(); + const vectorField& pCf = Cf.boundaryField()[patchi]; + + forAll(pOwner, pFacei) + { + label own = pOwner[pFacei]; + + limitFace + ( + limiter[own], + maxVsf[own], + minVsf[own], + ((pCf[pFacei] - C[own]) & g[own]) + ); + } + } + + if (fv::debug) + { + Info<< "gradient limiter for: " << vsf.name() + << " max = " << gMax(limiter) + << " min = " << gMin(limiter) + << " average: " << gAverage(limiter) << endl; + } + + limitGradient(limiter, g); + g.correctBoundaryConditions(); + gaussGrad::correctBoundaryConditions(vsf, g); + + return tGrad; +} + + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrad.H b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrad.H index 8b31b5d97b..462c3e63b5 100644 --- a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrad.H +++ b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrad.H @@ -44,14 +44,13 @@ SourceFiles #define cellLimitedGrad_H #include "gradScheme.H" +#include "Field.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // namespace Foam { -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - namespace fv { @@ -59,10 +58,11 @@ namespace fv Class cellLimitedGrad Declaration \*---------------------------------------------------------------------------*/ -template +template class cellLimitedGrad : - public fv::gradScheme + public fv::gradScheme, + public Limiter { // Private Data @@ -74,6 +74,18 @@ class cellLimitedGrad // Private Member Functions + void limitGradient + ( + const Field& limiter, + Field& gIf + ) const; + + void limitGradient + ( + const Field& limiter, + Field& gIf + ) const; + //- Disallow default bitwise copy construct cellLimitedGrad(const cellLimitedGrad&); @@ -93,6 +105,7 @@ public: cellLimitedGrad(const fvMesh& mesh, Istream& schemeData) : gradScheme(mesh), + Limiter(schemeData), basicGradScheme_(fv::gradScheme::New(mesh, schemeData)), k_(readScalar(schemeData)) { @@ -110,13 +123,21 @@ public: // Member Functions - static inline void limitFace + inline void limitFaceCmpt + ( + scalar& limiter, + const scalar maxDelta, + const scalar minDelta, + const scalar extrapolate + ) const; + + inline void limitFace ( Type& limiter, const Type& maxDelta, const Type& minDelta, const Type& extrapolate - ); + ) const; //- Return the gradient of the given field to the gradScheme::grad // for optional caching @@ -134,73 +155,67 @@ public: // * * * * * * * * * * * * Inline Member Function * * * * * * * * * * * * * // -template<> -inline void cellLimitedGrad::limitFace +template +inline void cellLimitedGrad::limitFaceCmpt ( scalar& limiter, - const scalar& maxDelta, - const scalar& minDelta, - const scalar& extrapolate -) + const scalar maxDelta, + const scalar minDelta, + const scalar extrapolate +) const { - if (extrapolate > maxDelta + VSMALL) + scalar r = 1; + + if (extrapolate > SMALL) { - limiter = min(limiter, maxDelta/extrapolate); + r = maxDelta/extrapolate; } - else if (extrapolate < minDelta - VSMALL) + else if (extrapolate < -SMALL) { - limiter = min(limiter, minDelta/extrapolate); + r = minDelta/extrapolate; } + else + { + return; + } + + limiter = min(limiter, Limiter::limiter(r)); } -template -inline void cellLimitedGrad::limitFace +template +inline void cellLimitedGrad::limitFace ( Type& limiter, const Type& maxDelta, const Type& minDelta, const Type& extrapolate -) +) const { - for (direction cmpt=0; cmpt::nComponents; ++cmpt) { - cellLimitedGrad::limitFace + limitFaceCmpt ( - limiter.component(cmpt), - maxDelta.component(cmpt), - minDelta.component(cmpt), - extrapolate.component(cmpt) + setComponent(limiter, cmpt), + component(maxDelta, cmpt), + component(minDelta, cmpt), + component(extrapolate, cmpt) ); } } -// * * * * * * * * Template Member Function Specialisations * * * * * * * * // - -template<> -tmp cellLimitedGrad::calcGrad -( - const volScalarField& vsf, - const word& name -) const; - - -template<> -tmp cellLimitedGrad::calcGrad -( - const volVectorField& vsf, - const word& name -) const; - - // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace fv +} // End namespace Foam + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -} // End namespace Foam +#ifdef NoRepository + #include "cellLimitedGrad.C" +#endif // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrads.C b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrads.C index 7a13f98077..356beaaeff 100644 --- a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrads.C +++ b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrads.C @@ -24,335 +24,49 @@ License \*---------------------------------------------------------------------------*/ #include "cellLimitedGrad.H" -#include "gaussGrad.H" -#include "fvMesh.H" -#include "volMesh.H" -#include "surfaceMesh.H" -#include "volFields.H" -#include "fixedValueFvPatchFields.H" +#include "minmodGradientLimiter.H" +#include "VenkatakrishnanGradientLimiter.H" +#include "cubicGradientLimiter.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -makeFvGradScheme(cellLimitedGrad) - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -template<> -Foam::tmp -Foam::fv::cellLimitedGrad::calcGrad -( - const volScalarField& vsf, - const word& name -) const -{ - const fvMesh& mesh = vsf.mesh(); - - tmp tGrad = basicGradScheme_().calcGrad(vsf, name); - - if (k_ < SMALL) - { - return tGrad; +#define makeNamedFvLimitedGradTypeScheme(SS, Type, Limiter, Name) \ + typedef Foam::fv::SS \ + SS##Type##Limiter##_; \ + \ + defineTemplateTypeNameAndDebugWithName \ + ( \ + SS##Type##Limiter##_, \ + Name, \ + 0 \ + ); \ + \ + namespace Foam \ + { \ + namespace fv \ + { \ + gradScheme::addIstreamConstructorToTable \ + < \ + SS \ + > add##SS##Type##Limiter##IstreamConstructorToTable_; \ + } \ } - volVectorField& g = tGrad.ref(); +#define makeFvLimitedGradTypeScheme(SS, Type, Limiter) \ + makeNamedFvLimitedGradTypeScheme(SS##Grad, Type, Limiter, #SS"<"#Limiter">") - const labelUList& owner = mesh.owner(); - const labelUList& neighbour = mesh.neighbour(); +#define makeFvLimitedGradScheme(SS, Limiter) \ + \ + makeFvLimitedGradTypeScheme(SS, scalar, Limiter) \ + makeFvLimitedGradTypeScheme(SS, vector, Limiter) - const volVectorField& C = mesh.C(); - const surfaceVectorField& Cf = mesh.Cf(); - scalarField maxVsf(vsf.primitiveField()); - scalarField minVsf(vsf.primitiveField()); - - forAll(owner, facei) - { - label own = owner[facei]; - label nei = neighbour[facei]; - - scalar vsfOwn = vsf[own]; - scalar vsfNei = vsf[nei]; - - maxVsf[own] = max(maxVsf[own], vsfNei); - minVsf[own] = min(minVsf[own], vsfNei); - - maxVsf[nei] = max(maxVsf[nei], vsfOwn); - minVsf[nei] = min(minVsf[nei], vsfOwn); - } - - - const volScalarField::Boundary& bsf = vsf.boundaryField(); - - forAll(bsf, patchi) - { - const fvPatchScalarField& psf = bsf[patchi]; - - const labelUList& pOwner = mesh.boundary()[patchi].faceCells(); - - if (psf.coupled()) - { - const scalarField psfNei(psf.patchNeighbourField()); - - forAll(pOwner, pFacei) - { - label own = pOwner[pFacei]; - scalar vsfNei = psfNei[pFacei]; - - maxVsf[own] = max(maxVsf[own], vsfNei); - minVsf[own] = min(minVsf[own], vsfNei); - } - } - else - { - forAll(pOwner, pFacei) - { - label own = pOwner[pFacei]; - scalar vsfNei = psf[pFacei]; - - maxVsf[own] = max(maxVsf[own], vsfNei); - minVsf[own] = min(minVsf[own], vsfNei); - } - } - } - - maxVsf -= vsf; - minVsf -= vsf; - - if (k_ < 1.0) - { - const scalarField maxMinVsf((1.0/k_ - 1.0)*(maxVsf - minVsf)); - maxVsf += maxMinVsf; - minVsf -= maxMinVsf; - - //maxVsf *= 1.0/k_; - //minVsf *= 1.0/k_; - } - - - // create limiter - scalarField limiter(vsf.primitiveField().size(), 1.0); - - forAll(owner, facei) - { - label own = owner[facei]; - label nei = neighbour[facei]; - - // owner side - limitFace - ( - limiter[own], - maxVsf[own], - minVsf[own], - (Cf[facei] - C[own]) & g[own] - ); - - // neighbour side - limitFace - ( - limiter[nei], - maxVsf[nei], - minVsf[nei], - (Cf[facei] - C[nei]) & g[nei] - ); - } - - forAll(bsf, patchi) - { - const labelUList& pOwner = mesh.boundary()[patchi].faceCells(); - const vectorField& pCf = Cf.boundaryField()[patchi]; - - forAll(pOwner, pFacei) - { - label own = pOwner[pFacei]; - - limitFace - ( - limiter[own], - maxVsf[own], - minVsf[own], - (pCf[pFacei] - C[own]) & g[own] - ); - } - } - - if (fv::debug) - { - Info<< "gradient limiter for: " << vsf.name() - << " max = " << gMax(limiter) - << " min = " << gMin(limiter) - << " average: " << gAverage(limiter) << endl; - } - - g.primitiveFieldRef() *= limiter; - g.correctBoundaryConditions(); - gaussGrad::correctBoundaryConditions(vsf, g); - - return tGrad; -} - - -template<> -Foam::tmp -Foam::fv::cellLimitedGrad::calcGrad -( - const volVectorField& vsf, - const word& name -) const -{ - const fvMesh& mesh = vsf.mesh(); - - tmp tGrad = basicGradScheme_().calcGrad(vsf, name); - - if (k_ < SMALL) - { - return tGrad; - } - - volTensorField& g = tGrad.ref(); - - const labelUList& owner = mesh.owner(); - const labelUList& neighbour = mesh.neighbour(); - - const volVectorField& C = mesh.C(); - const surfaceVectorField& Cf = mesh.Cf(); - - vectorField maxVsf(vsf.primitiveField()); - vectorField minVsf(vsf.primitiveField()); - - forAll(owner, facei) - { - label own = owner[facei]; - label nei = neighbour[facei]; - - const vector& vsfOwn = vsf[own]; - const vector& vsfNei = vsf[nei]; - - maxVsf[own] = max(maxVsf[own], vsfNei); - minVsf[own] = min(minVsf[own], vsfNei); - - maxVsf[nei] = max(maxVsf[nei], vsfOwn); - minVsf[nei] = min(minVsf[nei], vsfOwn); - } - - - const volVectorField::Boundary& bsf = vsf.boundaryField(); - - forAll(bsf, patchi) - { - const fvPatchVectorField& psf = bsf[patchi]; - const labelUList& pOwner = mesh.boundary()[patchi].faceCells(); - - if (psf.coupled()) - { - const vectorField psfNei(psf.patchNeighbourField()); - - forAll(pOwner, pFacei) - { - label own = pOwner[pFacei]; - const vector& vsfNei = psfNei[pFacei]; - - maxVsf[own] = max(maxVsf[own], vsfNei); - minVsf[own] = min(minVsf[own], vsfNei); - } - } - else - { - forAll(pOwner, pFacei) - { - label own = pOwner[pFacei]; - const vector& vsfNei = psf[pFacei]; - - maxVsf[own] = max(maxVsf[own], vsfNei); - minVsf[own] = min(minVsf[own], vsfNei); - } - } - } - - maxVsf -= vsf; - minVsf -= vsf; - - if (k_ < 1.0) - { - const vectorField maxMinVsf((1.0/k_ - 1.0)*(maxVsf - minVsf)); - maxVsf += maxMinVsf; - minVsf -= maxMinVsf; - - //maxVsf *= 1.0/k_; - //minVsf *= 1.0/k_; - } - - - // create limiter - vectorField limiter(vsf.primitiveField().size(), vector::one); - - forAll(owner, facei) - { - label own = owner[facei]; - label nei = neighbour[facei]; - - // owner side - limitFace - ( - limiter[own], - maxVsf[own], - minVsf[own], - (Cf[facei] - C[own]) & g[own] - ); - - // neighbour side - limitFace - ( - limiter[nei], - maxVsf[nei], - minVsf[nei], - (Cf[facei] - C[nei]) & g[nei] - ); - } - - forAll(bsf, patchi) - { - const labelUList& pOwner = mesh.boundary()[patchi].faceCells(); - const vectorField& pCf = Cf.boundaryField()[patchi]; - - forAll(pOwner, pFacei) - { - label own = pOwner[pFacei]; - - limitFace - ( - limiter[own], - maxVsf[own], - minVsf[own], - ((pCf[pFacei] - C[own]) & g[own]) - ); - } - } - - if (fv::debug) - { - Info<< "gradient limiter for: " << vsf.name() - << " max = " << gMax(limiter) - << " min = " << gMin(limiter) - << " average: " << gAverage(limiter) << endl; - } - - tensorField& gIf = g.primitiveFieldRef(); - - forAll(gIf, celli) - { - gIf[celli] = tensor - ( - cmptMultiply(limiter[celli], gIf[celli].x()), - cmptMultiply(limiter[celli], gIf[celli].y()), - cmptMultiply(limiter[celli], gIf[celli].z()) - ); - } - - g.correctBoundaryConditions(); - gaussGrad::correctBoundaryConditions(vsf, g); - - return tGrad; -} +// Default limiter in minmod specified without the limiter name +// for backward compatibility +makeNamedFvLimitedGradTypeScheme(cellLimitedGrad, scalar, minmod, "cellLimited") +makeNamedFvLimitedGradTypeScheme(cellLimitedGrad, vector, minmod, "cellLimited") +makeFvLimitedGradScheme(cellLimited, Venkatakrishnan) +makeFvLimitedGradScheme(cellLimited, cubic) // ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/gradientLimiters/VenkatakrishnanGradientLimiter.H b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/gradientLimiters/VenkatakrishnanGradientLimiter.H new file mode 100644 index 0000000000..ece67cb60d --- /dev/null +++ b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/gradientLimiters/VenkatakrishnanGradientLimiter.H @@ -0,0 +1,113 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2018 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM 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. + + OpenFOAM 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 OpenFOAM. If not, see . + +Class + Foam::fv::gradientLimiters::Venkatakrishnan + +Description + Venkatakrishnan gradient limiter + + to be used with the Foam::fv::cellLimitedGrad limited gradient. + + Reference: + \verbatim + Venkatakrishnan, V. (1993). + On the accuracy of limiters and convergence to steady state solutions. + In 31st Aerospace Sciences Meeting (p. 880). + \endverbatim + + Example: + \verbatim + gradSchemes + { + default Gauss linear; + limited cellLimited Gauss linear 1; + } + \endverbatim + + Note: this limiter formally allows the limiter function to slightly exceed 1 + which is generally not a good idea and can cause unboundedness. To avoid + this problem the limiter function is clipped to 1 which is likely to be + better behaved but is then not differentiable and so no longer conforms to + the basic principles of this kind of limiter function. All these issues are + resolved in a consistent and differentiable manner by the + Foam::fv::gradientLimiters::cubic limiter which is recommended in + preference to the Foam::fv::gradientLimiters::Venkatakrishnan limiter. + +See also + Foam::fv::cellLimitedGrad + Foam::fv::gradientLimiters::cubic + +\*---------------------------------------------------------------------------*/ + +#ifndef VenkatakrishnanGradientLimiter_H +#define VenkatakrishnanGradientLimiter_H + +#include "Istream.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +namespace fv +{ + +namespace gradientLimiters +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +class Venkatakrishnan +{ + +public: + + // Constructors + + Venkatakrishnan(Istream&) + {} + + + // Member Functions + + inline scalar limiter(const scalar r) const + { + return (sqr(r) + 2*r)/(sqr(r) + r + 2); + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace gradientLimiters + +} // End namespace fv + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/gradientLimiters/cubicGradientLimiter.H b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/gradientLimiters/cubicGradientLimiter.H new file mode 100644 index 0000000000..9270d1b6d0 --- /dev/null +++ b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/gradientLimiters/cubicGradientLimiter.H @@ -0,0 +1,146 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2018 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM 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. + + OpenFOAM 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 OpenFOAM. If not, see . + +Class + Foam::fv::gradientLimiters::minmod + +Description + Cubic gradient limiter + + to be used with the Foam::fv::cellLimitedGrad limited gradient. This + limiter function is similar to Foam::fv::gradientLimiters::Venkatakrishnan + but is a fit to obey the value and gradient constraints and avoids the + problem of the limiter exceeding 1 present in the Venkatakrishnan function. + + The transition point at which the limiter function reaches 1 is an input + parameter and should be set to a value between 1 and 2 although values + larger than 2 are physical but likely to significantly reduce the accuracy + of the scheme. + + Reference: + \verbatim + Michalak, K., & Ollivier-Gooch, C. (2008). + Limiters for unstructured higher-order accurate solutions + of the Euler equations. + In 46th AIAA Aerospace Sciences Meeting and Exhibit (p. 776). + \endverbatim + + Example: + \verbatim + gradSchemes + { + default Gauss linear; + limited cellLimited 1.5 Gauss linear 1; + } + \endverbatim + +See also + Foam::fv::cellLimitedGrad + Foam::fv::gradientLimiters::Venkatakrishnan + +\*---------------------------------------------------------------------------*/ + +#ifndef cubicGradientLimiter_H +#define cubicGradientLimiter_H + +#include "Istream.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +namespace fv +{ + +namespace gradientLimiters +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +class cubic +{ + // Private Data + + //- Limiter transition point at which the limiter function -> 1 + // Must be > 1 + const scalar rt_; + + //- Coefficient of the r^3 term (evaluated from rt) + const scalar a_; + + // - Coefficient of the r^2 term (evaluated from rt) + const scalar b_; + + +public: + + // Constructors + + cubic(Istream& schemeData) + : + rt_(readScalar(schemeData)), + a_(2.0/sqr(rt_) - 2.0/pow3(rt_)), + b_(-(3.0/2.0)*a_*rt_) + { + if (rt_ < 1) + { + FatalIOErrorInFunction + ( + schemeData + ) << "coefficient = " << rt_ + << " should be > 1" + << exit(FatalIOError); + } + } + + + // Member Functions + + inline scalar limiter(const scalar r) const + { + if (r < rt_) + { + return ((a_*r + b_)*r + 1)*r; + } + else + { + return 1; + } + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace gradientLimiters + +} // End namespace fv + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/gradientLimiters/minmodGradientLimiter.H b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/gradientLimiters/minmodGradientLimiter.H new file mode 100644 index 0000000000..33e8984476 --- /dev/null +++ b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/gradientLimiters/minmodGradientLimiter.H @@ -0,0 +1,99 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2018 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM 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. + + OpenFOAM 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 OpenFOAM. If not, see . + +Class + Foam::fv::gradientLimiters::minmod + +Description + Minmod gradient limiter + + This is the default gradient limiter which clips the gradient to remove cell + to face extrapolation unboundedness. It corresponds to the original + cellLimitedGrad implementation before the addition of run-time select-able + limiter functions. + + Example: + \verbatim + gradSchemes + { + default Gauss linear; + limited cellLimited Gauss linear 1; + } + \endverbatim + +See also + Foam::fv::cellLimitedGrad + +\*---------------------------------------------------------------------------*/ + +#ifndef minmodGradientLimiter_H +#define minmodGradientLimiter_H + +#include "Istream.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +namespace fv +{ + +namespace gradientLimiters +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +class minmod +{ + +public: + + // Constructors + + minmod(Istream&) + {} + + + // Member Functions + + inline scalar limiter(const scalar r) const + { + return min(r, 1); + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace gradientLimiters + +} // End namespace fv + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* //