Commit 86dc9556 authored by Henry Weller's avatar Henry Weller
Browse files

interFoam family of solvers: Improved Crank-Nicolson implementation

Fewer limiter iterations are now required to obtain sufficient boundedness and
restart is more consistent.
parent 093ede5d
......@@ -2,46 +2,57 @@
word alphaScheme("div(phi,alpha)");
word alpharScheme("div(phirb,alpha)");
tmp<fv::ddtScheme<scalar>> ddtAlpha
(
fv::ddtScheme<scalar>::New
(
mesh,
mesh.ddtScheme("ddt(alpha)")
)
);
// Set the off-centering coefficient according to ddt scheme
scalar ocCoeff = 0;
if
(
isType<fv::EulerDdtScheme<scalar>>(ddtAlpha())
|| isType<fv::localEulerDdtScheme<scalar>>(ddtAlpha())
)
{
ocCoeff = 0;
}
else if (isType<fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha()))
{
if (nAlphaSubCycles > 1)
tmp<fv::ddtScheme<scalar>> tddtAlpha
(
fv::ddtScheme<scalar>::New
(
mesh,
mesh.ddtScheme("ddt(alpha)")
)
);
const fv::ddtScheme<scalar>& ddtAlpha = tddtAlpha();
if
(
isType<fv::EulerDdtScheme<scalar>>(ddtAlpha)
|| isType<fv::localEulerDdtScheme<scalar>>(ddtAlpha)
)
{
ocCoeff = 0;
}
else if (isType<fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha))
{
if (nAlphaSubCycles > 1)
{
FatalErrorInFunction
<< "Sub-cycling is not supported "
"with the CrankNicolson ddt scheme"
<< exit(FatalError);
}
if
(
alphaRestart
|| mesh.time().timeIndex() > mesh.time().startTimeIndex() + 1
)
{
ocCoeff =
refCast<const fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha)
.ocCoeff();
}
}
else
{
FatalErrorInFunction
<< "Sub-cycling is not supported "
"with the CrankNicolson ddt scheme"
<< "Only Euler and CrankNicolson ddt schemes are supported"
<< exit(FatalError);
}
ocCoeff =
refCast<const fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha())
.ocCoeff();
}
else
{
FatalErrorInFunction
<< "Only Euler and CrankNicolson ddt schemes are supported"
<< exit(FatalError);
}
// Set the time blending factor, 1 for Euler
scalar cnCoeff = 1.0/(1.0 + ocCoeff);
// Standard face-flux compression coefficient
......@@ -136,8 +147,8 @@
(
fvc::flux
(
phi,
alpha1,
phiCN(),
cnCoeff*alpha1 + (1.0 - cnCoeff)*alpha1.oldTime(),
alphaScheme
)
+ fvc::flux
......@@ -148,13 +159,6 @@
)
);
// Calculate the Crank-Nicolson off-centred alpha flux
if (ocCoeff > 0)
{
talphaPhiUn =
cnCoeff*talphaPhiUn + (1.0 - cnCoeff)*alphaPhi.oldTime();
}
if (MULESCorr)
{
tmp<surfaceScalarField> talphaPhiCorr(talphaPhiUn() - alphaPhi);
......
IOobject alphaPhiHeader
(
"alphaPhi",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
);
const bool alphaRestart = alphaPhiHeader.headerOk();
// MULES flux from previous time-step
surfaceScalarField alphaPhi
(
alphaPhiHeader,
phi*fvc::interpolate(alpha1)
);
// MULES Correction
tmp<surfaceScalarField> talphaPhiCorr0;
......@@ -64,6 +64,7 @@ int main(int argc, char *argv[])
#include "initContinuityErrs.H"
#include "createControl.H"
#include "createFields.H"
#include "createAlphaFluxes.H"
#include "createFvOptions.H"
#include "createUf.H"
#include "createControls.H"
......
......@@ -61,6 +61,7 @@ int main(int argc, char *argv[])
#include "createControl.H"
#include "createTimeControls.H"
#include "createFields.H"
#include "createAlphaFluxes.H"
#include "createFvOptions.H"
volScalarField& p = mixture.p();
......
......@@ -101,21 +101,4 @@ autoPtr<compressible::turbulenceModel> turbulence
Info<< "Creating field kinetic energy K\n" << endl;
volScalarField K("K", 0.5*magSqr(U));
// MULES flux from previous time-step
surfaceScalarField alphaPhi
(
IOobject
(
"alphaPhi",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
phi*fvc::interpolate(alpha1)
);
// MULES Correction
tmp<surfaceScalarField> talphaPhiCorr0;
#include "createMRF.H"
......@@ -121,21 +121,4 @@ if (p_rgh.needReference())
mesh.setFluxRequired(p_rgh.name());
mesh.setFluxRequired(alpha1.name());
// MULES flux from previous time-step
surfaceScalarField alphaPhi
(
IOobject
(
"alphaPhi",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
phi*fvc::interpolate(alpha1)
);
// MULES Correction
tmp<surfaceScalarField> talphaPhiCorr0;
#include "createMRF.H"
......@@ -60,6 +60,7 @@ int main(int argc, char *argv[])
#include "createTimeControls.H"
#include "createDyMControls.H"
#include "createFields.H"
#include "createAlphaFluxes.H"
#include "createFvOptions.H"
volScalarField rAU
......
......@@ -63,6 +63,7 @@ int main(int argc, char *argv[])
#include "createTimeControls.H"
#include "initContinuityErrs.H"
#include "createFields.H"
#include "createAlphaFluxes.H"
#include "createFvOptions.H"
#include "correctPhi.H"
......
......@@ -198,7 +198,7 @@ scalar CrankNicolsonDdtScheme<Type>::coef_
const DDt0Field<GeoField>& ddt0
) const
{
if (mesh().time().timeIndex() - ddt0.startTimeIndex() > 0)
if (mesh().time().timeIndex() > ddt0.startTimeIndex())
{
return 1 + ocCoeff_;
}
......@@ -216,7 +216,7 @@ scalar CrankNicolsonDdtScheme<Type>::coef0_
const DDt0Field<GeoField>& ddt0
) const
{
if (mesh().time().timeIndex() - ddt0.startTimeIndex() > 1)
if (mesh().time().timeIndex() > ddt0.startTimeIndex() + 1)
{
return 1 + ocCoeff_;
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment