Commit bab72d51 authored by Johan Roenby's avatar Johan Roenby Committed by Andrew Heather
Browse files

ENH: multiple updates to interIsoFoam related code

- Updated tutorial headers
- Added copyright note to isoAdvector src
- Removed outcommented code lines in interIsoFoam solver
- Removed all LTS from interIsoFoam since this is not currently supported
- Confirmed that discInConstantFlow gives identical results with N subCylces and time step N*dt
- Confirmed that this also holds when nOuterCorrectors > 1.
parent 83b8032d
const dictionary& alphaControls = mesh.solverDict(alpha1.name());
//label nAlphaCorr(readLabel(alphaControls.lookup("nAlphaCorr")));
label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
/*
bool MULESCorr(alphaControls.lookupOrDefault("MULESCorr", false));
// Apply the compression correction from the previous iteration
// Improves efficiency for steady-simulations but can only be applied
// once the alpha field is reasonably steady, i.e. fully developed
bool alphaApplyPrevCorr
(
alphaControls.lookupOrDefault("alphaApplyPrevCorr", false)
);
// Isotropic compression coefficient
scalar icAlpha
(
alphaControls.lookupOrDefault<scalar>("icAlpha", 0)
);
// Shear compression coefficient
scalar scAlpha
(
alphaControls.lookupOrDefault<scalar>("scAlpha", 0)
);
*/
{
// Note for AMR implementation:
// At this point we have just entered the new time step,
// the mesh has been refined and the alpha, phi and U contain
// the field values at the beginning of this time step mapped
// to the new mesh.
// The above assumes that we are in firstIter() of the outer
// corrector loop. If we are in any subsequent iter of this loop
// the alpha1, U and phi will be overwritten with the new time step
// values but still on the same mesh.
if (pimple.firstIter())
{
// Note for AMR implementation:
// At this point we have just entered the new time step,
// the mesh has been refined and the alpha, phi and U contain
// the field values at the beginning of this time step mapped
// to the new mesh.
// The above assumes that we are in firstIter() of the outer
// corrector loop. If we are in any subsequent iter of this loop
// the alpha1, U and phi will be overwritten with the new time step
// values but still on the same mesh.
if (pimple.firstIter())
{
// Note: This assumes moveMeshOuterCorrectors is false
alpha10 = alpha1;
U0 = U;
phi0 = phi;
advector.advect();
#include "rhofs.H"
rhoPhi = advector.getRhoPhi(rho1f, rho2f);
}
else
{
alpha1 = alpha10;
// Temporarily setting U and phi to average of old and new value
// Note: Illegal additions if mesh is topoChanging
// (e.g. if moveMeshOuterCorrectors and AMR)
U = 0.5*U0 + 0.5*U;
phi = 0.5*phi0 + 0.5*phi;
isoAdvection advector(alpha1, phi, U);
advector.advect();
#include "rhofs.H"
rhoPhi = advector.getRhoPhi(rho1f, rho2f);
// Resetting U and phi to the new value
U = 2.0*U - U0;
phi = 2.0*phi - phi0;
}
alpha2 = 1.0 - alpha1;
mixture.correct();
// Note: This assumes moveMeshOuterCorrectors is false
// Saving field values before advection in first PIMPLE iteration
alpha10 = alpha1;
U0 = U;
phi0 = phi;
}
else
{
// Resetting alpha1 to value before advection in first PIMPLE iteration
alpha1 = alpha10;
// Temporarily setting U and phi to average of old and new value
// Note: Illegal additions if mesh is topoChanging
// (e.g. if moveMeshOuterCorrectors and AMR)
U = 0.5*U0 + 0.5*U;
phi = 0.5*phi0 + 0.5*phi;
}
advector.advect();
#include "rhofs.H"
rhoPhi = advector.getRhoPhi(rho1f, rho2f);
if (!pimple.firstIter())
{
// Resetting U and phi to the new value
U = 2.0*U - U0;
phi = 2.0*phi - phi0;
}
Info<< "Phase-1 volume fraction = "
<< alpha1.weightedAverage(mesh.Vsc()).value()
<< " Min(" << alpha1.name() << ") = " << min(alpha1).value()
<< " Max(" << alpha1.name() << ") = " << max(alpha1).value()
<< endl;
alpha2 = 1.0 - alpha1;
mixture.correct();
}
Info<< "Phase-1 volume fraction = "
<< alpha1.weightedAverage(mesh.Vsc()).value()
<< " Min(" << alpha1.name() << ") = " << min(alpha1).value()
<< " Max(" << alpha1.name() << ") = " << max(alpha1).value()
<< endl;
......@@ -13,14 +13,6 @@ if (nAlphaSubCycles > 1)
dimensionedScalar(rhoPhi.dimensions(), Zero)
);
tmp<volScalarField> trSubDeltaT;
if (LTS)
{
trSubDeltaT =
fv::localEulerDdt::localRSubDeltaT(mesh, nAlphaSubCycles);
}
for
(
subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
......
......@@ -121,23 +121,6 @@ if (p_rgh.needReference())
mesh.setFluxRequired(p_rgh.name());
mesh.setFluxRequired(alpha1.name());
/*
// MULES compressed flux is registered in case scalarTransport FO needs it.
surfaceScalarField alphaPhiUn
(
IOobject
(
"alphaPhiUn",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar(phi.dimensions(), Zero)
);
*/
#include "createMRF.H"
#include "createFvOptions.H"
......@@ -155,7 +138,6 @@ surfaceScalarField phi0
phi
);
volVectorField U0
(
IOobject
......@@ -169,7 +151,6 @@ volVectorField U0
U
);
volScalarField alpha10
(
IOobject
......@@ -183,5 +164,4 @@ volScalarField alpha10
alpha1
);
isoAdvection advector(alpha1, phi, U);
......@@ -75,17 +75,13 @@ int main(int argc, char *argv[])
#include "initContinuityErrs.H"
#include "createDyMControls.H"
#include "createFields.H"
// #include "createAlphaFluxes.H"
#include "initCorrectPhi.H"
#include "createUfIfPresent.H"
turbulence->validate();
if (!LTS)
{
#include "CourantNo.H"
#include "setInitialDeltaT.H"
}
#include "CourantNo.H"
#include "setInitialDeltaT.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
......@@ -93,17 +89,9 @@ int main(int argc, char *argv[])
while (runTime.run())
{
#include "readDyMControls.H"
if (LTS)
{
#include "setRDeltaT.H"
}
else
{
#include "CourantNo.H"
#include "alphaCourantNo.H"
#include "setDeltaT.H"
}
#include "CourantNo.H"
#include "alphaCourantNo.H"
#include "setDeltaT.H"
runTime++;
......@@ -118,12 +106,6 @@ int main(int argc, char *argv[])
if (mesh.changing())
{
// Do not apply previous time-step mesh compression flux
// if the mesh topology changed
if (mesh.topoChanging())
{
// talphaPhi1Corr0.clear();
}
gh = (g & mesh.C()) - ghRef;
ghf = (g & mesh.Cf()) - ghRef;
......
{
volScalarField& rDeltaT = trDeltaT.ref();
const dictionary& pimpleDict = pimple.dict();
scalar maxCo
(
pimpleDict.lookupOrDefault<scalar>("maxCo", 0.9)
);
scalar maxAlphaCo
(
pimpleDict.lookupOrDefault<scalar>("maxAlphaCo", 0.2)
);
scalar rDeltaTSmoothingCoeff
(
pimpleDict.lookupOrDefault<scalar>("rDeltaTSmoothingCoeff", 0.1)
);
label nAlphaSpreadIter
(
pimpleDict.lookupOrDefault<label>("nAlphaSpreadIter", 1)
);
scalar alphaSpreadDiff
(
pimpleDict.lookupOrDefault<scalar>("alphaSpreadDiff", 0.2)
);
scalar alphaSpreadMax
(
pimpleDict.lookupOrDefault<scalar>("alphaSpreadMax", 0.99)
);
scalar alphaSpreadMin
(
pimpleDict.lookupOrDefault<scalar>("alphaSpreadMin", 0.01)
);
label nAlphaSweepIter
(
pimpleDict.lookupOrDefault<label>("nAlphaSweepIter", 5)
);
scalar rDeltaTDampingCoeff
(
pimpleDict.lookupOrDefault<scalar>("rDeltaTDampingCoeff", 1.0)
);
scalar maxDeltaT
(
pimpleDict.lookupOrDefault<scalar>("maxDeltaT", GREAT)
);
volScalarField rDeltaT0("rDeltaT0", rDeltaT);
// Set the reciprocal time-step from the local Courant number
rDeltaT.ref() = max
(
1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT),
fvc::surfaceSum(mag(rhoPhi))()()
/((2*maxCo)*mesh.V()*rho())
);
if (maxAlphaCo < maxCo)
{
// Further limit the reciprocal time-step
// in the vicinity of the interface
volScalarField alpha1Bar(fvc::average(alpha1));
rDeltaT.ref() = max
(
rDeltaT(),
pos0(alpha1Bar() - alphaSpreadMin)
*pos0(alphaSpreadMax - alpha1Bar())
*fvc::surfaceSum(mag(phi))()()
/((2*maxAlphaCo)*mesh.V())
);
}
// Update tho boundary values of the reciprocal time-step
rDeltaT.correctBoundaryConditions();
Info<< "Flow time scale min/max = "
<< gMin(1/rDeltaT.primitiveField())
<< ", " << gMax(1/rDeltaT.primitiveField()) << endl;
if (rDeltaTSmoothingCoeff < 1.0)
{
fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
}
if (nAlphaSpreadIter > 0)
{
fvc::spread
(
rDeltaT,
alpha1,
nAlphaSpreadIter,
alphaSpreadDiff,
alphaSpreadMax,
alphaSpreadMin
);
}
if (nAlphaSweepIter > 0)
{
fvc::sweep(rDeltaT, alpha1, nAlphaSweepIter, alphaSpreadDiff);
}
Info<< "Smoothed flow time scale min/max = "
<< gMin(1/rDeltaT.primitiveField())
<< ", " << gMax(1/rDeltaT.primitiveField()) << endl;
// Limit rate of change of time scale
// - reduce as much as required
// - only increase at a fraction of old time scale
if
(
rDeltaTDampingCoeff < 1.0
&& runTime.timeIndex() > runTime.startTimeIndex() + 1
)
{
rDeltaT = max
(
rDeltaT,
(scalar(1.0) - rDeltaTDampingCoeff)*rDeltaT0
);
Info<< "Damped flow time scale min/max = "
<< gMin(1/rDeltaT.primitiveField())
<< ", " << gMax(1/rDeltaT.primitiveField()) << endl;
}
}
......@@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
isoAdvector | Copyright (C) 2016-2017 DHI
Modified work | Copyright (C) 2018 Johan Roenby
-------------------------------------------------------------------------------
License
......
......@@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
isoAdvector | Copyright (C) 2016-2017 DHI
Modified work | Copyright (C) 2018 Johan Roenby
-------------------------------------------------------------------------------
License
......
......@@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
isoAdvector | Copyright (C) 2016-2017 DHI
Modified work | Copyright (C) 2018 Johan Roenby
-------------------------------------------------------------------------------
License
......
......@@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
isoAdvector | Copyright (C) 2016-2017 DHI
Modified work | Copyright (C) 2018 Johan Roenby
-------------------------------------------------------------------------------
License
......
......@@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
isoAdvector | Copyright (C) 2016-2017 DHI
Modified work | Copyright (C) 2018 Johan Roenby
-------------------------------------------------------------------------------
License
......
......@@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
isoAdvector | Copyright (C) 2016-2017 DHI
Modified work | Copyright (C) 2018 Johan Roenby
-------------------------------------------------------------------------------
License
......
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 5 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\ / O peration | Version: plus |
| \\ / A nd | Web: www.OpenFOAM.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
......
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 5 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\ / O peration | Version: plus |
| \\ / A nd | Web: www.OpenFOAM.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
......
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 5 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\ / O peration | Version: plus |
| \\ / A nd | Web: www.OpenFOAM.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
......
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 5 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\ / O peration | Version: plus |
| \\ / A nd | Web: www.OpenFOAM.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
......
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 5 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\ / O peration | Version: plus |
| \\ / A nd | Web: www.OpenFOAM.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
......
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 5 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\ / O peration | Version: plus |
| \\ / A nd | Web: www.OpenFOAM.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
......
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 5 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\ / O peration | Version: plus |
| \\ / A nd | Web: www.OpenFOAM.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
......
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