INT: ddtPhiCorr - reinstated v1712 behaviour and provided experimental

version on a switch.  See #867

By default the code will use the same form as previous versions

To use the experimental version integrated from openfoam.org commit
da787200 set the info switch in the controlDict:

    InfoSwitches
    {
        experimentalDdtCorr 1;
    }

src/finiteVolume/Make/files

@@ -358,6 +358,7 @@ finiteVolume/fv/fv.C
 finiteVolume/fvSchemes/fvSchemes.C
 
 ddtSchemes = finiteVolume/ddtSchemes
+$(ddtSchemes)/ddtScheme/ddtSchemeBase.C
 $(ddtSchemes)/ddtScheme/ddtSchemes.C
 $(ddtSchemes)/steadyStateDdtScheme/steadyStateDdtSchemes.C
 $(ddtSchemes)/EulerDdtScheme/EulerDdtSchemes.C

src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.C

@@ -28,6 +28,7 @@ License
 #include "surfaceInterpolate.H"
 #include "fvMatrix.H"
 #include "cyclicAMIFvPatch.H"
+#include "registerSwitch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -146,6 +147,11 @@ tmp<surfaceScalarField> ddtScheme<Type>::fvcDdtPhiCoeff
     const fluxFieldType& phiCorr
 )
 {
+    if (fv::debug)
+    {
+        InfoInFunction << "Using standard version" << endl;
+    }
+
     tmp<surfaceScalarField> tddtCouplingCoeff
     (
         new surfaceScalarField
@@ -166,18 +172,84 @@ tmp<surfaceScalarField> ddtScheme<Type>::fvcDdtPhiCoeff
     if (ddtPhiCoeff_ < 0)
     {
         // v1712 and earlier
-        // ddtCouplingCoeff -= min
+        ddtCouplingCoeff -= min
-        // (
+        (
-        //     mag(phiCorr)
+            mag(phiCorr)
-        //    /(mag(phi) + dimensionedScalar("small", phi.dimensions(), SMALL)),
+           /(mag(phi) + dimensionedScalar("small", phi.dimensions(), SMALL)),
-        //     scalar(1)
+            scalar(1)
-        // );
+        );
+    }
+    else
+    {
+        ddtCouplingCoeff =
+            dimensionedScalar("ddtPhiCoeff", dimless, ddtPhiCoeff_);
+    }
 
+    surfaceScalarField::Boundary& ccbf = ddtCouplingCoeff.boundaryFieldRef();
+
+    forAll(U.boundaryField(), patchi)
+    {
+        if
+        (
+            U.boundaryField()[patchi].fixesValue()
+         || isA<cyclicAMIFvPatch>(mesh().boundary()[patchi])
+        )
+        {
+            ccbf[patchi] = 0.0;
+        }
+    }
+
+    if (debug > 1)
+    {
+        InfoInFunction
+            << "ddtCouplingCoeff mean max min = "
+            << gAverage(ddtCouplingCoeff.primitiveField())
+            << " " << gMax(ddtCouplingCoeff.primitiveField())
+            << " " << gMin(ddtCouplingCoeff.primitiveField())
+            << endl;
+    }
+
+    return tddtCouplingCoeff;
+}
+
+
+template<class Type>
+tmp<surfaceScalarField> ddtScheme<Type>::fvcDdtPhiCoeffExperimental
+(
+    const GeometricField<Type, fvPatchField, volMesh>& U,
+    const fluxFieldType& phi,
+    const fluxFieldType& phiCorr
+)
+{
+    if (fv::debug)
+    {
+        InfoInFunction << "Using experimental version" << endl;
+    }
+
+    tmp<surfaceScalarField> tddtCouplingCoeff
+    (
+        new surfaceScalarField
+        (
+            IOobject
+            (
+                "ddtCouplingCoeff",
+                U.mesh().time().timeName(),
+                U.mesh()
+            ),
+            U.mesh(),
+            dimensionedScalar("one", dimless, 1.0)
+        )
+    );
+
+    surfaceScalarField& ddtCouplingCoeff = tddtCouplingCoeff.ref();
+
+    if (ddtPhiCoeff_ < 0)
+    {
         // See note below re: commented code
         ddtCouplingCoeff -= min
         (
             //  mag(phiCorr)
-            //*mesh().time().deltaT()*mag(mesh().deltaCoeffs())/mesh().magSf(),
+            // *mesh().time().deltaT()*mag(mesh().deltaCoeffs())/mesh().magSf(),
             //  scalar(1)
             mag(phiCorr)
             *mesh().time().deltaT()*mesh().deltaCoeffs()/mesh().magSf(),
@@ -231,7 +303,20 @@ tmp<surfaceScalarField> ddtScheme<Type>::fvcDdtPhiCoeff
     const volScalarField& rho
 )
 {
-    return fvcDdtPhiCoeff(U, phi, phiCorr/fvc::interpolate(rho));
+    if (experimentalDdtCorr)
+    {
+        return
+            fvcDdtPhiCoeffExperimental
+            (
+                U,
+                phi,
+                phiCorr/fvc::interpolate(rho)
+            );
+    }
+    else
+    {
+        return fvcDdtPhiCoeff(U, phi, phiCorr);
+    }
 }
 
 
@@ -242,7 +327,26 @@ tmp<surfaceScalarField> ddtScheme<Type>::fvcDdtPhiCoeff
     const fluxFieldType& phi
 )
 {
-    return fvcDdtPhiCoeff(U, phi, phi - fvc::dotInterpolate(mesh().Sf(), U));
+    if (experimentalDdtCorr)
+    {
+        return
+            fvcDdtPhiCoeffExperimental
+            (
+                U,
+                phi,
+                phi - fvc::dotInterpolate(mesh().Sf(), U)
+            );
+    }
+    else
+    {
+        return
+            fvcDdtPhiCoeff
+            (
+                U,
+                phi,
+                phi - fvc::dotInterpolate(mesh().Sf(), U)
+            );
+    }
 }
 
 
@@ -254,12 +358,25 @@ tmp<surfaceScalarField> ddtScheme<Type>::fvcDdtPhiCoeff
     const volScalarField& rho
 )
 {
+    if (experimentalDdtCorr)
+    {
+        return fvcDdtPhiCoeffExperimental
+        (
+            U,
+            phi,
+            (phi - fvc::dotInterpolate(mesh().Sf(), rho*U))
+           /fvc::interpolate(rho)
+        );
+    }
+    else
+    {
         return fvcDdtPhiCoeff
         (
             U,
             phi,
-        (phi - fvc::dotInterpolate(mesh().Sf(), U))/fvc::interpolate(rho)
+            (phi - fvc::dotInterpolate(mesh().Sf(), rho*U))
         );
+    }
 }
 
 

src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.H

@@ -64,10 +64,24 @@ namespace fv
                            Class ddtScheme Declaration
 \*---------------------------------------------------------------------------*/
 
+class ddtSchemeBase
+{
+public:
+
+        //- Flag to use experimental ddtCorr from org version
+        //- Default is off for backwards compatibility
+        static bool experimentalDdtCorr;
+
+        ddtSchemeBase()
+        {}
+};
+
+
 template<class Type>
 class ddtScheme
 :
-    public refCount
+    public refCount,
+    public ddtSchemeBase
 {
 
 protected:
@@ -76,7 +90,9 @@ protected:
 
         const fvMesh& mesh_;
 
-        //- Input for fvcDdtPhiCoeff (-1 default)
+        //- Input for fvcDdtPhiCoeff
+        //  If set to -1 (default) the code will calculate the coefficient
+        //  If > 0 the coupling coeff is set to this value
         scalar ddtPhiCoeff_;
 
 
@@ -218,6 +234,13 @@ public:
             const fluxFieldType& phiCorr
         );
 
+        tmp<surfaceScalarField> fvcDdtPhiCoeffExperimental
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& U,
+            const fluxFieldType& phi,
+            const fluxFieldType& phiCorr
+        );
+
         tmp<surfaceScalarField> fvcDdtPhiCoeff
         (
             const GeometricField<Type, fvPatchField, volMesh>& U,