From 91ee15e3ac8153fa0af79d3ffad0d73f5517a45c Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Fri, 19 Jun 2015 11:52:48 +0100
Subject: [PATCH] rhoCentralFoam: Added experimental LTS support

Select LTS via the ddtScheme:

ddtSchemes
{
    default         localEuler rDeltaT;
}

The LTS algorithm is 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.

For cases with reasonably uniform meshes like the forwardStep tutorial
LTS does not provide much benefit but for cases with large variation in
cell-size like the biconic25-55Run35 tutorial LTS provides significant
speed-up to convergence particularly if started from uniform conditions.
---
 .../rhoCentralFoam/createRDeltaT.H            | 28 ++++++++++++++++++
 .../rhoCentralFoam/rhoCentralFoam.C           | 13 ++++++++-
 .../compressible/rhoCentralFoam/setRDeltaT.H  | 29 +++++++++++++++++++
 3 files changed, 69 insertions(+), 1 deletion(-)
 create mode 100644 applications/solvers/compressible/rhoCentralFoam/createRDeltaT.H
 create mode 100644 applications/solvers/compressible/rhoCentralFoam/setRDeltaT.H

diff --git a/applications/solvers/compressible/rhoCentralFoam/createRDeltaT.H b/applications/solvers/compressible/rhoCentralFoam/createRDeltaT.H
new file mode 100644
index 0000000000..44d181383c
--- /dev/null
+++ b/applications/solvers/compressible/rhoCentralFoam/createRDeltaT.H
@@ -0,0 +1,28 @@
+bool LTS =
+    word(mesh.ddtScheme("ddt(rho)"))
+ == fv::localEulerDdtScheme<scalar>::typeName;
+
+tmp<volScalarField> trDeltaT;
+
+if (LTS)
+{
+    Info<< "Using LTS" << endl;
+
+    trDeltaT = tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "rDeltaT",
+                runTime.timeName(),
+                mesh,
+                IOobject::READ_IF_PRESENT,
+                IOobject::AUTO_WRITE
+            ),
+            mesh,
+            dimensionedScalar("one", dimless/dimTime, 1),
+            zeroGradientFvPatchScalarField::typeName
+        )
+    );
+}
diff --git a/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C b/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C
index 0b5cba1c78..6c0ec3856f 100644
--- a/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C
+++ b/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C
@@ -36,6 +36,8 @@ Description
 #include "zeroGradientFvPatchFields.H"
 #include "fixedRhoFvPatchScalarField.H"
 #include "directionInterpolate.H"
+#include "localEulerDdtScheme.H"
+#include "fvcSmooth.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -46,6 +48,7 @@ int main(int argc, char *argv[])
     #include "createTime.H"
     #include "createMesh.H"
     #include "createFields.H"
+    #include "createRDeltaT.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -134,7 +137,15 @@ int main(int argc, char *argv[])
 
         #include "centralCourantNo.H"
         #include "readTimeControls.H"
-        #include "setDeltaT.H"
+
+        if (LTS)
+        {
+            #include "setRDeltaT.H"
+        }
+        else
+        {
+            #include "setDeltaT.H"
+        }
 
         runTime++;
 
diff --git a/applications/solvers/compressible/rhoCentralFoam/setRDeltaT.H b/applications/solvers/compressible/rhoCentralFoam/setRDeltaT.H
new file mode 100644
index 0000000000..bb7ebdef4f
--- /dev/null
+++ b/applications/solvers/compressible/rhoCentralFoam/setRDeltaT.H
@@ -0,0 +1,29 @@
+{
+    volScalarField& rDeltaT = trDeltaT();
+
+    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(amaxSf)().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;
+}