From 1d16db3fa0c69c15d688d2e5da975959967f7300 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Fri, 26 Jun 2015 18:52:23 +0100
Subject: [PATCH] reactingTwoPhaseEulerFoam: Added experimental run-time
 selectable LTS support

Select LTS via the ddtScheme:

    ddtSchemes
    {
        default         localEuler rDeltaT;
    }

The LTS algorithm is currently controlled with the standard settings in
controlDict, e.g.:

    maxCo           0.5;
    maxDeltaT       2e-8;

with the addition of the optional rDeltaT smoothing coefficient:

    rDeltaTSmoothingCoeff 0.02;

which defaults to 0.02.

    ddtSchemes
    {
        default         localEuler rDeltaT;
    }
---
 .../reactingTwoPhaseEulerFoam/createRDeltaT.H | 47 +++++++++++++++++++
 .../reactingTwoPhaseEulerFoam.C               | 25 ++++++++--
 .../reactingTwoPhaseEulerFoam/setRDeltaT.H    | 35 ++++++++++++++
 3 files changed, 102 insertions(+), 5 deletions(-)
 create mode 100644 applications/solvers/multiphase/reactingTwoPhaseEulerFoam/createRDeltaT.H
 create mode 100644 applications/solvers/multiphase/reactingTwoPhaseEulerFoam/setRDeltaT.H

diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/createRDeltaT.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/createRDeltaT.H
new file mode 100644
index 00000000000..b45dff137e4
--- /dev/null
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/createRDeltaT.H
@@ -0,0 +1,47 @@
+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)
+        )
+    );
+}
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C
index 27826f12a19..b0434b45d67 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C
@@ -38,6 +38,8 @@ Description
 #include "PhaseCompressibleTurbulenceModel.H"
 #include "fixedFluxPressureFvPatchScalarField.H"
 #include "pimpleControl.H"
+#include "localEulerDdtScheme.H"
+#include "fvcSmooth.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -50,11 +52,16 @@ int main(int argc, char *argv[])
 
     pimpleControl pimple(mesh);
 
+    #include "createRDeltaT.H"
     #include "createFields.H"
     #include "initContinuityErrs.H"
-    #include "readTimeControls.H"
-    #include "CourantNos.H"
-    #include "setInitialDeltaT.H"
+
+    if (!LTS)
+    {
+        #include "readTimeControls.H"
+        #include "CourantNo.H"
+        #include "setInitialDeltaT.H"
+    }
 
     Switch faceMomentum
     (
@@ -78,8 +85,16 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
         #include "readTimeControls.H"
-        #include "CourantNos.H"
-        #include "setDeltaT.H"
+
+        if (LTS)
+        {
+            #include "setRDeltaT.H"
+        }
+        else
+        {
+            #include "CourantNos.H"
+            #include "setDeltaT.H"
+        }
 
         runTime++;
         Info<< "Time = " << runTime.timeName() << nl << endl;
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/setRDeltaT.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/setRDeltaT.H
new file mode 100644
index 00000000000..99e08c3f565
--- /dev/null
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/setRDeltaT.H
@@ -0,0 +1,35 @@
+{
+    volScalarField& rDeltaT = trDeltaT();
+    volScalarField& rSubDeltaT = trSubDeltaT();
+
+    scalar rDeltaTSmoothingCoeff
+    (
+        runTime.controlDict().lookupOrDefault<scalar>
+        (
+            "rDeltaTSmoothingCoeff",
+            0.02
+        )
+    );
+
+    // Set the reciprocal time-step from the local Courant number
+    rDeltaT.dimensionedInternalField() = max
+    (
+        1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT),
+        fvc::surfaceSum(max(mag(phi1), mag(phi2)))().dimensionedInternalField()
+       /((2*maxCo)*mesh.V())
+    );
+
+    // Update tho boundary values of the reciprocal time-step
+    rDeltaT.correctBoundaryConditions();
+
+    fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
+
+    Info<< "Flow time scale min/max = "
+        << gMin(1/rDeltaT.internalField())
+        << ", " << gMax(1/rDeltaT.internalField()) << endl;
+
+    const dictionary& alphaControls = mesh.solverDict(alpha1.name());
+    label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
+
+    rSubDeltaT = rDeltaT*nAlphaSubCycles;
+}
-- 
GitLab