ENH: cellLimitedGrad gradientLimiters: Added support for run-time selectable gradient limiter function

Minmod is the default limiter function and specified with an explicit name e.g.:

    gradSchemes
    {
        default Gauss linear;
        limited cellLimited Gauss linear 1;
    }

Venkatakrishnan and cubic limiter functions are also provided and may be
specified explicitly e.g.:

    gradSchemes
    {
        default Gauss linear;
        limited cellLimited<Venkatakrishnan> Gauss linear 1;
    }

or

    gradSchemes
    {
        default Gauss linear;
        limited cellLimited<cubic> 1.5 Gauss linear 1;
    }

The standard minmod function is recommended for most applications but if
convergence or stability problems arise it may be beneficial to use one of the
alternatives which smooth the gradient limiting.  The Venkatakrishnan is not
well formulated and allows the limiter to exceed 1 whereas the cubic limiter is
designed to obey all the value and gradient constraints on the limiter function,
see

    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).

The cubic limiter function requires the transition point at which the limiter
function reaches 1 is an input parameter which 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.

VenkatakrishnanGradientLimiter: Updated documentation

cubicGradientLimiter: Documented private data
This commit is contained in:
Henry Weller
2018-02-26 23:14:46 +00:00
committed by Andrew Heather
parent c5125c3299
commit 08079c4be9
6 changed files with 682 additions and 365 deletions

View File

@ -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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "cellLimitedGrad.H"
#include "gaussGrad.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type, class Limiter>
void Foam::fv::cellLimitedGrad<Type, Limiter>::limitGradient
(
const Field<scalar>& limiter,
Field<vector>& gIf
) const
{
gIf *= limiter;
}
template<class Type, class Limiter>
void Foam::fv::cellLimitedGrad<Type, Limiter>::limitGradient
(
const Field<vector>& limiter,
Field<tensor>& 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<class Type, class Limiter>
Foam::tmp
<
Foam::GeometricField
<
typename Foam::outerProduct<Foam::vector, Type>::type,
Foam::fvPatchField,
Foam::volMesh
>
>
Foam::fv::cellLimitedGrad<Type, Limiter>::calcGrad
(
const GeometricField<Type, fvPatchField, volMesh>& vsf,
const word& name
) const
{
const fvMesh& mesh = vsf.mesh();
tmp
<
GeometricField
<typename outerProduct<vector, Type>::type, fvPatchField, volMesh>
> tGrad = basicGradScheme_().calcGrad(vsf, name);
if (k_ < SMALL)
{
return tGrad;
}
GeometricField
<
typename outerProduct<vector, Type>::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<Type> maxVsf(vsf.primitiveField());
Field<Type> 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<Type, fvPatchField, volMesh>::Boundary& bsf =
vsf.boundaryField();
forAll(bsf, patchi)
{
const fvPatchField<Type>& psf = bsf[patchi];
const labelUList& pOwner = mesh.boundary()[patchi].faceCells();
if (psf.coupled())
{
const Field<Type> 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<Type> 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<Type> limiter(vsf.primitiveField().size(), pTraits<Type>::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<Type>::correctBoundaryConditions(vsf, g);
return tGrad;
}
// ************************************************************************* //

View File

@ -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<class Type>
template<class Type, class Limiter>
class cellLimitedGrad
:
public fv::gradScheme<Type>
public fv::gradScheme<Type>,
public Limiter
{
// Private Data
@ -74,6 +74,18 @@ class cellLimitedGrad
// Private Member Functions
void limitGradient
(
const Field<scalar>& limiter,
Field<vector>& gIf
) const;
void limitGradient
(
const Field<vector>& limiter,
Field<tensor>& gIf
) const;
//- Disallow default bitwise copy construct
cellLimitedGrad(const cellLimitedGrad&);
@ -93,6 +105,7 @@ public:
cellLimitedGrad(const fvMesh& mesh, Istream& schemeData)
:
gradScheme<Type>(mesh),
Limiter(schemeData),
basicGradScheme_(fv::gradScheme<Type>::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<scalar>::limitFace
template<class Type, class Limiter>
inline void cellLimitedGrad<Type, Limiter>::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<class Type>
inline void cellLimitedGrad<Type>::limitFace
template<class Type, class Limiter>
inline void cellLimitedGrad<Type, Limiter>::limitFace
(
Type& limiter,
const Type& maxDelta,
const Type& minDelta,
const Type& extrapolate
)
) const
{
for (direction cmpt=0; cmpt<Type::nComponents; cmpt++)
for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; ++cmpt)
{
cellLimitedGrad<scalar>::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<volVectorField> cellLimitedGrad<scalar>::calcGrad
(
const volScalarField& vsf,
const word& name
) const;
template<>
tmp<volTensorField> cellLimitedGrad<vector>::calcGrad
(
const volVectorField& vsf,
const word& name
) const;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace fv
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
#ifdef NoRepository
#include "cellLimitedGrad.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View File

@ -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::volVectorField>
Foam::fv::cellLimitedGrad<Foam::scalar>::calcGrad
(
const volScalarField& vsf,
const word& name
) const
{
const fvMesh& mesh = vsf.mesh();
tmp<volVectorField> tGrad = basicGradScheme_().calcGrad(vsf, name);
if (k_ < SMALL)
{
return tGrad;
#define makeNamedFvLimitedGradTypeScheme(SS, Type, Limiter, Name) \
typedef Foam::fv::SS<Foam::Type, Foam::fv::gradientLimiters::Limiter> \
SS##Type##Limiter##_; \
\
defineTemplateTypeNameAndDebugWithName \
( \
SS##Type##Limiter##_, \
Name, \
0 \
); \
\
namespace Foam \
{ \
namespace fv \
{ \
gradScheme<Type>::addIstreamConstructorToTable \
< \
SS<Type, gradientLimiters::Limiter> \
> 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<scalar>::correctBoundaryConditions(vsf, g);
return tGrad;
}
template<>
Foam::tmp<Foam::volTensorField>
Foam::fv::cellLimitedGrad<Foam::vector>::calcGrad
(
const volVectorField& vsf,
const word& name
) const
{
const fvMesh& mesh = vsf.mesh();
tmp<volTensorField> 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<vector>::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)
// ************************************************************************* //

View File

@ -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 <http://www.gnu.org/licenses/>.
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<Venkatakrishnan> 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
// ************************************************************************* //

View File

@ -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 <http://www.gnu.org/licenses/>.
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<cubic> 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
// ************************************************************************* //

View File

@ -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 <http://www.gnu.org/licenses/>.
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
// ************************************************************************* //