Commit e00ae44a authored by Henry Weller's avatar Henry Weller
Browse files

interFoam family: Added run-time selectable LTS support

LTS is selected by the ddt scheme e.g. in the
tutorials/multiphase/interFoam/ras/DTCHull case:

ddtSchemes
{
    default         localEuler rDeltaT;
}

LTSInterFoam is no longer needed now that interFoam includes LTS
support.
parent bf95b5c2
......@@ -15,7 +15,7 @@
(
1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT),
fvc::surfaceSum(amaxSf)().dimensionedInternalField()
/(2*maxCo*mesh.V())
/((2*maxCo)*mesh.V())
);
// Update tho boundary values of the reciprocal time-step
......
EXE_INC = \
-I$(LIB_SRC)/transportModels/compressible/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude
LIB_LIBS = \
-lcompressibleTransportModels \
-lfluidThermophysicalModels \
-lspecie \
-lfiniteVolume
-lfiniteVolume \
-lmeshTools
......@@ -6,11 +6,7 @@
{
fvScalarMatrix alpha1Eqn
(
#ifdef LTSSOLVE
fv::localEulerDdtScheme<scalar>(mesh, rDeltaT.name()).fvmDdt(alpha1)
#else
fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
#endif
+ fv::gaussConvectionScheme<scalar>
(
mesh,
......@@ -34,16 +30,6 @@
{
Info<< "Applying the previous iteration correction flux" << endl;
#ifdef LTSSOLVE
MULES::LTScorrect
(
alpha1,
phiAlpha,
tphiAlphaCorr0(),
mixture.alphaMax(),
0
);
#else
MULES::correct
(
alpha1,
......@@ -52,7 +38,6 @@
mixture.alphaMax(),
0
);
#endif
phiAlpha += tphiAlphaCorr0();
}
......@@ -84,16 +69,6 @@
tmp<surfaceScalarField> tphiAlphaCorr(tphiAlphaUn() - phiAlpha);
volScalarField alpha10("alpha10", alpha1);
#ifdef LTSSOLVE
MULES::LTScorrect
(
alpha1,
tphiAlphaUn(),
tphiAlphaCorr(),
mixture.alphaMax(),
0
);
#else
MULES::correct
(
alpha1,
......@@ -102,7 +77,6 @@
mixture.alphaMax(),
0
);
#endif
// Under-relax the correction for all but the 1st corrector
if (aCorr == 0)
......@@ -119,16 +93,6 @@
{
phiAlpha = tphiAlphaUn;
#ifdef LTSSOLVE
MULES::explicitLTSSolve
(
alpha1,
phi,
phiAlpha,
mixture.alphaMax(),
0
);
#else
MULES::explicitSolve
(
alpha1,
......@@ -137,7 +101,6 @@
mixture.alphaMax(),
0
);
#endif
}
}
......
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 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/>.
Application
interFoam
Description
Solver for 2 incompressible, isothermal immiscible fluids using a VOF
(volume of fluid) phase-fraction based interface capturing approach.
The momentum and other fluid properties are of the "mixture" and a single
momentum equation is solved.
Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
For a two-fluid approach see twoPhaseEulerFoam.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "CMULES.H"
#include "EulerDdtScheme.H"
#include "localEulerDdtScheme.H"
#include "CrankNicolsonDdtScheme.H"
#include "subCycle.H"
#include "immiscibleIncompressibleTwoPhaseMixture.H"
#include "turbulentTransportModel.H"
#include "fvcSmooth.H"
#include "pimpleControl.H"
#include "fvIOoptionList.H"
#include "CorrectPhi.H"
#include "fixedFluxPressureFvPatchScalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
pimpleControl pimple(mesh);
#include "initContinuityErrs.H"
#include "createFields.H"
#include "createMRF.H"
#include "createFvOptions.H"
#include "correctPhi.H"
#include "CourantNo.H"
#include "createRDeltaT.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.run())
{
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
#include "setrDeltaT.H"
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
#include "alphaControls.H"
#define LTSSOLVE
#include "alphaEqnSubCycle.H"
#undef LTSSOLVE
mixture.correct();
turbulence->correct();
#include "UEqn.H"
// --- Pressure corrector loop
while (pimple.correct())
{
#include "pEqn.H"
}
}
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //
EXE_INC = \
-I.. \
-I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/fvOptions/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude
EXE_LIBS = \
-limmiscibleIncompressibleTwoPhaseMixture \
-lturbulenceModels \
-lincompressibleTurbulenceModels \
-lfiniteVolume \
-lmeshTools \
-lfvOptions \
-lsampling
scalar maxDeltaT
(
pimple.dict().lookupOrDefault<scalar>("maxDeltaT", GREAT)
);
volScalarField rDeltaT
(
IOobject
(
"rDeltaT",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT),
zeroGradientFvPatchScalarField::typeName
);
volScalarField rSubDeltaT
(
IOobject
(
"rSubDeltaT",
runTime.timeName(),
mesh
),
mesh,
1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT)
);
......@@ -77,11 +77,15 @@
{
fvScalarMatrix alpha1Eqn
(
#ifdef LTSSOLVE
fv::localEulerDdtScheme<scalar>(mesh, rDeltaT.name()).fvmDdt(alpha1)
#else
fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
#endif
(
LTS
? fv::localEulerDdtScheme<scalar>
(
mesh,
trSubDeltaT().name()
).fvmDdt(alpha1)
: fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
)
+ fv::gaussConvectionScheme<scalar>
(
mesh,
......@@ -104,11 +108,7 @@
if (alphaApplyPrevCorr && tphiAlphaCorr0.valid())
{
Info<< "Applying the previous iteration compression flux" << endl;
#ifdef LTSSOLVE
MULES::LTScorrect(alpha1, phiAlpha, tphiAlphaCorr0(), 1, 0);
#else
MULES::correct(alpha1, phiAlpha, tphiAlphaCorr0(), 1, 0);
#endif
phiAlpha += tphiAlphaCorr0();
}
......@@ -154,11 +154,7 @@
tmp<surfaceScalarField> tphiAlphaCorr(tphiAlphaUn() - phiAlpha);
volScalarField alpha10("alpha10", alpha1);
#ifdef LTSSOLVE
MULES::LTScorrect(alpha1, tphiAlphaUn(), tphiAlphaCorr(), 1, 0);
#else
MULES::correct(alpha1, tphiAlphaUn(), tphiAlphaCorr(), 1, 0);
#endif
// Under-relax the correction for all but the 1st corrector
if (aCorr == 0)
......@@ -175,11 +171,7 @@
{
phiAlpha = tphiAlphaUn;
#ifdef LTSSOLVE
MULES::explicitLTSSolve(alpha1, phi, phiAlpha, 1, 0);
#else
MULES::explicitSolve(alpha1, phiCN, phiAlpha, 1, 0);
#endif
}
alpha2 = 1.0 - alpha1;
......
bool LTS =
word(mesh.ddtScheme("default"))
== fv::localEulerDdtScheme<scalar>::typeName;
tmp<volScalarField> trDeltaT;
tmp<volScalarField> trSubDeltaT;
if (LTS)
{
scalar maxDeltaT
(
pimple.dict().lookupOrDefault<scalar>("maxDeltaT", GREAT)
);
trDeltaT = tmp<volScalarField>
(
new volScalarField
(
IOobject
(
"rDeltaT",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT),
zeroGradientFvPatchScalarField::typeName
)
);
trSubDeltaT = tmp<volScalarField>
(
new volScalarField
(
IOobject
(
"rSubDeltaT",
runTime.timeName(),
mesh
),
mesh,
1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT)
)
);
}
......@@ -45,6 +45,8 @@ Description
#include "fvIOoptionList.H"
#include "CorrectPhi.H"
#include "fixedFluxPressureFvPatchScalarField.H"
#include "localEulerDdtScheme.H"
#include "fvcSmooth.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -60,7 +62,6 @@ int main(int argc, char *argv[])
#include "createFields.H"
#include "createMRF.H"
#include "createFvOptions.H"
#include "readTimeControls.H"
volScalarField rAU
(
......@@ -78,8 +79,14 @@ int main(int argc, char *argv[])
#include "correctPhi.H"
#include "createUf.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"
#include "createRDeltaT.H"
if (!LTS)
{
#include "readTimeControls.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
......@@ -87,10 +94,17 @@ int main(int argc, char *argv[])
while (runTime.run())
{
#include "readControls.H"
#include "alphaCourantNo.H"
#include "CourantNo.H"
#include "setDeltaT.H"
if (LTS)
{
#include "setRDeltaT.H"
}
else
{
#include "CourantNo.H"
#include "alphaCourantNo.H"
#include "setDeltaT.H"
}
runTime++;
......
......@@ -49,6 +49,8 @@ Description
#include "fvIOoptionList.H"
#include "CorrectPhi.H"
#include "fixedFluxPressureFvPatchScalarField.H"
#include "localEulerDdtScheme.H"
#include "fvcSmooth.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -64,10 +66,15 @@ int main(int argc, char *argv[])
#include "createFields.H"
#include "createMRF.H"
#include "createFvOptions.H"
#include "readTimeControls.H"
#include "correctPhi.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"
#include "createRDeltaT.H"
if (!LTS)
{
#include "readTimeControls.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -76,9 +83,17 @@ int main(int argc, char *argv[])
while (runTime.run())
{
#include "readTimeControls.H"
#include "CourantNo.H"
#include "alphaCourantNo.H"
#include "setDeltaT.H"
if (LTS)
{
#include "setRDeltaT.H"
}
else
{
#include "CourantNo.H"
#include "alphaCourantNo.H"
#include "setDeltaT.H"
}
runTime++;
......
......@@ -72,19 +72,41 @@
phiAlpha1 -= phiAlpha1BD;
// Calculate the limiter for alpha1
MULES::limiter
(
allLambda,
1.0/runTime.deltaT().value(),
geometricOneField(),
alpha1,
phiAlpha1BD,
phiAlpha1,
zeroField(),
zeroField(),
1,
0
);
if (LTS)
{
const volScalarField& rDeltaT =
mesh.objectRegistry::lookupObject<volScalarField>("rSubDeltaT");
MULES::limiter
(
allLambda,
rDeltaT,
geometricOneField(),
alpha1,
phiAlpha1BD,
phiAlpha1,
zeroField(),
zeroField(),
1,
0
);
}
else
{
MULES::limiter
(
allLambda,
1.0/runTime.deltaT().value(),
geometricOneField(),
alpha1,
phiAlpha1BD,
phiAlpha1,
zeroField(),
zeroField(),
1,
0
);
}
// Create the complete flux for alpha2
surfaceScalarField phiAlpha2
......@@ -113,19 +135,41 @@
phiAlpha2 -= phiAlpha2BD;
// Further limit the limiter for alpha2
MULES::limiter
(
allLambda,
1.0/runTime.deltaT().value(),
geometricOneField(),
alpha2,
phiAlpha2BD,
phiAlpha2,
zeroField(),
zeroField(),
1,
0
);
if (LTS)
{
const volScalarField& rDeltaT =
mesh.objectRegistry::lookupObject<volScalarField>("rSubDeltaT");
MULES::limiter
(
allLambda,
rDeltaT,
geometricOneField(),
alpha2,
phiAlpha2BD,
phiAlpha2,
zeroField(),
zeroField(),
1,
0
);
}
else
{
MULES::limiter
(
allLambda,
1.0/runTime.deltaT().value(),
geometricOneField(),
alpha2,
phiAlpha2BD,
phiAlpha2,
zeroField(),
zeroField(),
1,
0
);
}
// Construct the limited fluxes
phiAlpha1 = phiAlpha1BD + lambda*phiAlpha1;
......
......@@ -39,6 +39,8 @@ Description
#include "fvIOoptionList.H"
#include "CorrectPhi.H"