From f316980e296272067c3af471f3136be49725eb16 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Wed, 22 Mar 2017 22:27:47 +0000
Subject: [PATCH] CrankNicolsonDdtScheme: Added option "ramp" function to
 smoothly transition from Euler

    Off-centering is specified via the mandatory coefficient \c ocCoeff in the
    range [0,1] following the scheme name e.g.
    \verbatim
    ddtSchemes
    {
        default         CrankNicolson 0.9;
    }
    \endverbatim
    or with an optional "ramp" function to transition from the Euler scheme to
    Crank-Nicolson over a initial period to avoid start-up problems, e.g.
    \verbatim
    ddtSchemes
    {
        default         CrankNicolson
        ocCoeff
        {
            type scale;
            scale linearRamp;
            duration 0.01;
            value 0.9;
        };
    }
    \endverbatim

Note this functionality is experimental and the specification and implementation
may change if issues arise.
---
 .../CrankNicolsonDdtScheme.C                  |  8 +--
 .../CrankNicolsonDdtScheme.H                  | 57 +++++++++++++++----
 2 files changed, 49 insertions(+), 16 deletions(-)

diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C
index dd12b2552..341cf88b0 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C
@@ -200,7 +200,7 @@ scalar CrankNicolsonDdtScheme<Type>::coef_
 {
     if (mesh().time().timeIndex() > ddt0.startTimeIndex())
     {
-        return 1 + ocCoeff_;
+        return 1 + ocCoeff();
     }
     else
     {
@@ -218,7 +218,7 @@ scalar CrankNicolsonDdtScheme<Type>::coef0_
 {
     if (mesh().time().timeIndex() > ddt0.startTimeIndex() + 1)
     {
-        return 1 + ocCoeff_;
+        return 1 + ocCoeff();
     }
     else
     {
@@ -256,9 +256,9 @@ tmp<GeoField> CrankNicolsonDdtScheme<Type>::offCentre_
     const GeoField& ddt0
 ) const
 {
-    if (ocCoeff_ < 1)
+    if (ocCoeff() < 1)
     {
-        return ocCoeff_*ddt0;
+        return ocCoeff()*ddt0;
     }
     else
     {
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.H b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.H
index 78c343c3b..4b002fd45 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.H
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.H
@@ -32,13 +32,28 @@ Description
     geometries and it is necessary to "off-centre" the scheme to stabilize it
     while retaining greater temporal accuracy than the first-order
     Euler-implicit scheme.  Off-centering is specified via the mandatory
-    coefficient in the range [0,1] following the scheme name e.g.
+    coefficient \c ocCoeff in the range [0,1] following the scheme name e.g.
     \verbatim
     ddtSchemes
     {
         default         CrankNicolson 0.9;
     }
     \endverbatim
+    or with an optional "ramp" function to transition from the Euler scheme to
+    Crank-Nicolson over a initial period to avoid start-up problems, e.g.
+    \verbatim
+    ddtSchemes
+    {
+        default         CrankNicolson
+        ocCoeff
+        {
+            type scale;
+            scale linearRamp;
+            duration 0.01;
+            value 0.9;
+        };
+    }
+    \endverbatim
     With a coefficient of 1 the scheme is fully centred and second-order,
     with a coefficient of 0 the scheme is equivalent to Euler-implicit.
     A coefficient of 0.9 has been found to be suitable for a range of cases for
@@ -76,6 +91,7 @@ SourceFiles
 #define CrankNicolsonDdtScheme_H
 
 #include "ddtScheme.H"
+#include "Constant.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -138,7 +154,7 @@ class CrankNicolsonDdtScheme
 
 
         //- Off-centering coefficient, 1 -> CN, less than one blends with EI
-        scalar ocCoeff_;
+        autoPtr<Function1<scalar>> ocCoeff_;
 
 
     // Private Member Functions
@@ -197,7 +213,7 @@ public:
         CrankNicolsonDdtScheme(const fvMesh& mesh)
         :
             ddtScheme<Type>(mesh),
-            ocCoeff_(1.0)
+            ocCoeff_(new Function1Types::Constant<scalar>("ocCoeff", 1))
         {
             // Ensure the old-old-time cell volumes are available
             // for moving meshes
@@ -210,17 +226,34 @@ public:
         //- Construct from mesh and Istream
         CrankNicolsonDdtScheme(const fvMesh& mesh, Istream& is)
         :
-            ddtScheme<Type>(mesh, is),
-            ocCoeff_(readScalar(is))
+            ddtScheme<Type>(mesh, is)
         {
-            if (ocCoeff_ < 0 || ocCoeff_ > 1)
+            token firstToken(is);
+
+            if (firstToken.isNumber())
             {
-                FatalIOErrorInFunction
+                const scalar ocCoeff = firstToken.scalarToken();
+                if (ocCoeff < 0 || ocCoeff > 1)
+                {
+                    FatalIOErrorInFunction
+                    (
+                        is
+                    )   << "Off-centreing coefficient = " << ocCoeff
+                        << " should be >= 0 and <= 1"
+                        << exit(FatalIOError);
+                }
+
+                ocCoeff_ = new Function1Types::Constant<scalar>
                 (
-                    is
-                )   << "Off-centreing coefficient = " << ocCoeff_
-                    << " should be >= 0 and <= 1"
-                    << exit(FatalIOError);
+                    "ocCoeff",
+                    ocCoeff
+                );
+            }
+            else
+            {
+                is.putBack(firstToken);
+                dictionary dict(is);
+                ocCoeff_ = Function1<scalar>::New("ocCoeff", dict);
             }
 
             // Ensure the old-old-time cell volumes are available
@@ -243,7 +276,7 @@ public:
         //- Return the off-centreing coefficient
         scalar ocCoeff() const
         {
-            return ocCoeff_;
+            return ocCoeff_->value(mesh().time().value());
         }
 
         tmp<GeometricField<Type, fvPatchField, volMesh>> fvcDdt