ENH: Added relaxed nonOrthogonal laplacian corrector scheme

This commit is contained in:
mattijs
2019-09-25 12:40:23 +01:00
committed by Andrew Heather
parent 391c1eaa40
commit 86da88cdfd
4 changed files with 627 additions and 0 deletions

View File

@ -467,6 +467,7 @@ $(convectionSchemes)/boundedConvectionScheme/boundedConvectionSchemes.C
laplacianSchemes = finiteVolume/laplacianSchemes
$(laplacianSchemes)/laplacianScheme/laplacianSchemes.C
$(laplacianSchemes)/gaussLaplacianScheme/gaussLaplacianSchemes.C
$(laplacianSchemes)/relaxedNonOrthoGaussLaplacianScheme/relaxedNonOrthoGaussLaplacianSchemes.C
finiteVolume/fvc/fvcFlux.C
finiteVolume/fvc/fvcMeshPhi.C

View File

@ -0,0 +1,276 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2019 OpenCFD Ltd.
-------------------------------------------------------------------------------
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 "relaxedNonOrthoGaussLaplacianScheme.H"
#include "surfaceInterpolate.H"
#include "fvcDiv.H"
#include "fvcGrad.H"
#include "fvMatrices.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace fv
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type, class GType>
tmp<fvMatrix<Type>>
relaxedNonOrthoGaussLaplacianScheme<Type, GType>::fvmLaplacianUncorrected
(
const surfaceScalarField& gammaMagSf,
const surfaceScalarField& deltaCoeffs,
const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
tmp<fvMatrix<Type>> tfvm
(
new fvMatrix<Type>
(
vf,
deltaCoeffs.dimensions()*gammaMagSf.dimensions()*vf.dimensions()
)
);
fvMatrix<Type>& fvm = tfvm.ref();
fvm.upper() = deltaCoeffs.primitiveField()*gammaMagSf.primitiveField();
fvm.negSumDiag();
forAll(vf.boundaryField(), patchi)
{
const fvPatchField<Type>& pvf = vf.boundaryField()[patchi];
const fvsPatchScalarField& pGamma = gammaMagSf.boundaryField()[patchi];
const fvsPatchScalarField& pDeltaCoeffs =
deltaCoeffs.boundaryField()[patchi];
if (pvf.coupled())
{
fvm.internalCoeffs()[patchi] =
pGamma*pvf.gradientInternalCoeffs(pDeltaCoeffs);
fvm.boundaryCoeffs()[patchi] =
-pGamma*pvf.gradientBoundaryCoeffs(pDeltaCoeffs);
}
else
{
fvm.internalCoeffs()[patchi] = pGamma*pvf.gradientInternalCoeffs();
fvm.boundaryCoeffs()[patchi] = -pGamma*pvf.gradientBoundaryCoeffs();
}
}
return tfvm;
}
template<class Type, class GType>
tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
relaxedNonOrthoGaussLaplacianScheme<Type, GType>::gammaSnGradCorr
(
const surfaceVectorField& SfGammaCorr,
const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
const fvMesh& mesh = this->mesh();
tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> tgammaSnGradCorr
(
new GeometricField<Type, fvsPatchField, surfaceMesh>
(
IOobject
(
"gammaSnGradCorr("+vf.name()+')',
vf.instance(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
SfGammaCorr.dimensions()
*vf.dimensions()*mesh.deltaCoeffs().dimensions()
)
);
for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
{
tgammaSnGradCorr.ref().replace
(
cmpt,
fvc::dotInterpolate(SfGammaCorr, fvc::grad(vf.component(cmpt)))
);
}
return tgammaSnGradCorr;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type, class GType>
tmp<GeometricField<Type, fvPatchField, volMesh>>
relaxedNonOrthoGaussLaplacianScheme<Type, GType>::fvcLaplacian
(
const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
const fvMesh& mesh = this->mesh();
tmp<GeometricField<Type, fvPatchField, volMesh>> tLaplacian
(
fvc::div(this->tsnGradScheme_().snGrad(vf)*mesh.magSf())
);
tLaplacian.ref().rename("laplacian(" + vf.name() + ')');
return tLaplacian;
}
template<class Type, class GType>
tmp<fvMatrix<Type>>
relaxedNonOrthoGaussLaplacianScheme<Type, GType>::fvmLaplacian
(
const GeometricField<GType, fvsPatchField, surfaceMesh>& gamma,
const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
const fvMesh& mesh = this->mesh();
typedef GeometricField<Type, fvsPatchField, surfaceMesh> SType;
const surfaceVectorField Sn(mesh.Sf()/mesh.magSf());
const surfaceVectorField SfGamma(mesh.Sf() & gamma);
const GeometricField<scalar, fvsPatchField, surfaceMesh> SfGammaSn
(
SfGamma & Sn
);
const surfaceVectorField SfGammaCorr(SfGamma - SfGammaSn*Sn);
tmp<fvMatrix<Type>> tfvm = fvmLaplacianUncorrected
(
SfGammaSn,
this->tsnGradScheme_().deltaCoeffs(vf),
vf
);
fvMatrix<Type>& fvm = tfvm.ref();
tmp<SType> tfaceFluxCorrection = gammaSnGradCorr(SfGammaCorr, vf);
if (this->tsnGradScheme_().corrected())
{
tfaceFluxCorrection.ref() +=
SfGammaSn*this->tsnGradScheme_().correction(vf);
}
const word corrName(tfaceFluxCorrection().name());
tmp<SType> trelaxedCorrection(new SType(tfaceFluxCorrection()));
const word oldName(corrName + "_0");
const scalar relax(vf.mesh().equationRelaxationFactor(oldName));
const objectRegistry& obr = vf.db();
if (obr.foundObject<SType>(oldName))
{
SType& oldCorrection = obr.lookupObjectRef<SType>(oldName);
trelaxedCorrection.ref() *= relax;
trelaxedCorrection.ref() += (1.0-relax)*oldCorrection;
oldCorrection = tfaceFluxCorrection;
}
else
{
SType* s = new SType(oldName, tfaceFluxCorrection);
s->store();
}
fvm.source() -=
mesh.V()
*fvc::div
(
trelaxedCorrection()
)().primitiveField();
if (mesh.fluxRequired(vf.name()))
{
fvm.faceFluxCorrectionPtr() = trelaxedCorrection.ptr();
}
return tfvm;
}
template<class Type, class GType>
tmp<GeometricField<Type, fvPatchField, volMesh>>
relaxedNonOrthoGaussLaplacianScheme<Type, GType>::fvcLaplacian
(
const GeometricField<GType, fvsPatchField, surfaceMesh>& gamma,
const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
const fvMesh& mesh = this->mesh();
const surfaceVectorField Sn(mesh.Sf()/mesh.magSf());
const surfaceVectorField SfGamma(mesh.Sf() & gamma);
const GeometricField<scalar, fvsPatchField, surfaceMesh> SfGammaSn
(
SfGamma & Sn
);
const surfaceVectorField SfGammaCorr(SfGamma - SfGammaSn*Sn);
tmp<GeometricField<Type, fvPatchField, volMesh>> tLaplacian
(
fvc::div
(
SfGammaSn*this->tsnGradScheme_().snGrad(vf)
+ gammaSnGradCorr(SfGammaCorr, vf)
)
);
tLaplacian.ref().rename
(
"laplacian(" + gamma.name() + ',' + vf.name() + ')'
);
return tLaplacian;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace fv
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View File

@ -0,0 +1,210 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2019 OpenCFD Ltd.
-------------------------------------------------------------------------------
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::relaxedNonOrthoGaussLaplacianScheme
Description
Basic second-order laplacian using face-gradients and Gauss' theorem.
Usage
Minimal example by using \c system/fvSchemes:
\verbatim
laplacianSchemes
{
laplacian(<term>) relaxedNonOrthoGauss <other options>;
}
\endverbatim
and by using \c system/fvSolution:
\verbatim
relaxationFactors
{
equations
{
<term> <relaxation factor>;
}
}
\endverbatim
SourceFiles
relaxedNonOrthoGaussLaplacianScheme.C
\*---------------------------------------------------------------------------*/
#ifndef relaxedNonOrthoGaussLaplacianScheme_H
#define relaxedNonOrthoGaussLaplacianScheme_H
#include "laplacianScheme.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace fv
{
/*---------------------------------------------------------------------------*\
Class relaxedNonOrthoGaussLaplacianScheme Declaration
\*---------------------------------------------------------------------------*/
template<class Type, class GType>
class relaxedNonOrthoGaussLaplacianScheme
:
public fv::laplacianScheme<Type, GType>
{
// Private Member Functions
tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> gammaSnGradCorr
(
const surfaceVectorField& SfGammaCorr,
const GeometricField<Type, fvPatchField, volMesh>&
);
//- No copy construct
relaxedNonOrthoGaussLaplacianScheme
(
const relaxedNonOrthoGaussLaplacianScheme&
) = delete;
//- No copy assignment
void operator=(const relaxedNonOrthoGaussLaplacianScheme&) = delete;
public:
//- Runtime type information
TypeName("relaxedNonOrthoGauss");
// Constructors
//- Construct null
relaxedNonOrthoGaussLaplacianScheme(const fvMesh& mesh)
:
laplacianScheme<Type, GType>(mesh)
{}
//- Construct from Istream
relaxedNonOrthoGaussLaplacianScheme(const fvMesh& mesh, Istream& is)
:
laplacianScheme<Type, GType>(mesh, is)
{}
//- Construct from mesh, interpolation and snGradScheme schemes
relaxedNonOrthoGaussLaplacianScheme
(
const fvMesh& mesh,
const tmp<surfaceInterpolationScheme<GType>>& igs,
const tmp<snGradScheme<Type>>& sngs
)
:
laplacianScheme<Type, GType>(mesh, igs, sngs)
{}
//- Destructor
virtual ~relaxedNonOrthoGaussLaplacianScheme() = default;
// Member Functions
static tmp<fvMatrix<Type>> fvmLaplacianUncorrected
(
const surfaceScalarField& gammaMagSf,
const surfaceScalarField& deltaCoeffs,
const GeometricField<Type, fvPatchField, volMesh>&
);
tmp<GeometricField<Type, fvPatchField, volMesh>> fvcLaplacian
(
const GeometricField<Type, fvPatchField, volMesh>&
);
tmp<fvMatrix<Type>> fvmLaplacian
(
const GeometricField<GType, fvsPatchField, surfaceMesh>&,
const GeometricField<Type, fvPatchField, volMesh>&
);
tmp<GeometricField<Type, fvPatchField, volMesh>> fvcLaplacian
(
const GeometricField<GType, fvsPatchField, surfaceMesh>&,
const GeometricField<Type, fvPatchField, volMesh>&
);
};
// Use macros to emulate partial-specialisation of the Laplacian functions
// for scalar diffusivity gamma
#define defineFvmLaplacianScalarGamma(Type) \
\
template<> \
tmp<fvMatrix<Type>> \
relaxedNonOrthoGaussLaplacianScheme<Type, scalar>::fvmLaplacian \
( \
const GeometricField<scalar, fvsPatchField, surfaceMesh>&, \
const GeometricField<Type, fvPatchField, volMesh>& \
); \
\
template<> \
tmp<GeometricField<Type, fvPatchField, volMesh>> \
relaxedNonOrthoGaussLaplacianScheme<Type, scalar>::fvcLaplacian \
( \
const GeometricField<scalar, fvsPatchField, surfaceMesh>&, \
const GeometricField<Type, fvPatchField, volMesh>& \
);
defineFvmLaplacianScalarGamma(scalar);
defineFvmLaplacianScalarGamma(vector);
defineFvmLaplacianScalarGamma(sphericalTensor);
defineFvmLaplacianScalarGamma(symmTensor);
defineFvmLaplacianScalarGamma(tensor);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace fv
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "relaxedNonOrthoGaussLaplacianScheme.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,140 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2019 OpenCFD Ltd.
-------------------------------------------------------------------------------
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 "relaxedNonOrthoGaussLaplacianScheme.H"
#include "fvMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
makeFvLaplacianScheme(relaxedNonOrthoGaussLaplacianScheme)
#define declareFvmLaplacianScalarGamma(Type) \
\
template<> \
Foam::tmp<Foam::fvMatrix<Foam::Type>> \
Foam::fv::relaxedNonOrthoGaussLaplacianScheme<Foam::Type, Foam::scalar>:: \
fvmLaplacian \
( \
const GeometricField<scalar, fvsPatchField, surfaceMesh>& gamma, \
const GeometricField<Type, fvPatchField, volMesh>& vf \
) \
{ \
const fvMesh& mesh = this->mesh(); \
\
typedef GeometricField<Type, fvsPatchField, surfaceMesh> SType; \
\
GeometricField<scalar, fvsPatchField, surfaceMesh> gammaMagSf \
( \
gamma*mesh.magSf() \
); \
\
tmp<fvMatrix<Type>> tfvm = fvmLaplacianUncorrected \
( \
gammaMagSf, \
this->tsnGradScheme_().deltaCoeffs(vf), \
vf \
); \
fvMatrix<Type>& fvm = tfvm.ref(); \
\
if (this->tsnGradScheme_().corrected()) \
{ \
tmp<SType> tCorr(this->tsnGradScheme_().correction(vf)); \
const word corrName(tCorr().name()); \
tmp<SType> tfaceFluxCorrection(gammaMagSf*tCorr); \
\
tmp<SType> trelaxedCorrection(new SType(tfaceFluxCorrection())); \
\
const word oldName(corrName + "_0"); \
const scalar relax(vf.mesh().equationRelaxationFactor(corrName)); \
const objectRegistry& obr = vf.db(); \
if (obr.foundObject<SType>(oldName)) \
{ \
SType& oldCorrection = obr.lookupObjectRef<SType>(oldName); \
trelaxedCorrection.ref() *= relax; \
trelaxedCorrection.ref() += (1.0-relax)*oldCorrection; \
\
oldCorrection = trelaxedCorrection(); \
} \
else \
{ \
SType* s = new SType(oldName, tfaceFluxCorrection); \
s->store(); \
} \
\
tmp<Field<Type>> tcorr \
( \
mesh.V() \
*fvc::div \
( \
trelaxedCorrection() \
)().primitiveField() \
); \
\
fvm.source() -= tcorr(); \
\
if (mesh.fluxRequired(vf.name())) \
{ \
fvm.faceFluxCorrectionPtr() = trelaxedCorrection.ptr(); \
} \
} \
\
return tfvm; \
} \
\
\
template<> \
Foam::tmp<Foam::GeometricField<Foam::Type, Foam::fvPatchField, Foam::volMesh>> \
Foam::fv::relaxedNonOrthoGaussLaplacianScheme<Foam::Type, Foam::scalar>::fvcLaplacian \
( \
const GeometricField<scalar, fvsPatchField, surfaceMesh>& gamma, \
const GeometricField<Type, fvPatchField, volMesh>& vf \
) \
{ \
const fvMesh& mesh = this->mesh(); \
\
tmp<GeometricField<Type, fvPatchField, volMesh>> tLaplacian \
( \
fvc::div(gamma*this->tsnGradScheme_().snGrad(vf)*mesh.magSf()) \
); \
\
tLaplacian.ref().rename \
( \
"laplacian(" + gamma.name() + ',' + vf.name() + ')' \
); \
\
return tLaplacian; \
}
declareFvmLaplacianScalarGamma(scalar);
declareFvmLaplacianScalarGamma(vector);
declareFvmLaplacianScalarGamma(sphericalTensor);
declareFvmLaplacianScalarGamma(symmTensor);
declareFvmLaplacianScalarGamma(tensor);
// ************************************************************************* //