Run-time bugfixes regarding time consistency

src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C

@@ -175,16 +175,17 @@ const Foam::dimensionedScalar Foam::pimpleControl::relaxFactor
     const Foam::volVectorField& U
 ) const
 {
-    return
+    scalar urf = 1;
-        dimensionedScalar
+
-        (
+    if (mesh_.solutionDict().relaxEquation(U.select(finalIter())))
-            "alphaU",
+    {
-            dimless,
+        urf = mesh_.solutionDict().equationRelaxationFactor
-            mesh_.solutionDict().equationRelaxationFactor
             (
                 U.select(finalIter())
-            )
+            );
-        );
+    }
+
+    return dimensionedScalar("alphaU", dimless, urf);
 }
 
 

src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.H

@@ -91,6 +91,15 @@ protected:
         virtual bool criteriaSatisfied();
 
 
+        // Time and under-relaxation consistency
+
+            //- Get relaxation factor for velocity field
+            virtual const dimensionedScalar relaxFactor
+            (
+                const volVectorField& U
+            ) const;
+
+
 public:
 
     // Static Data Members
@@ -123,15 +132,6 @@ public:
             inline label corrPISO() const;
 
 
-        // Time and under-relaxation consistency
-
-            //- Get relaxation factor for velocity field
-            virtual const dimensionedScalar relaxFactor
-            (
-                const volVectorField& U
-            ) const;
-
-
         // Solution control
 
             //- PIMPLE loop

src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C

@@ -222,6 +222,22 @@ Foam::scalar Foam::solutionControl::maxResidual
 }
 
 
+const Foam::dimensionedScalar Foam::solutionControl::relaxFactor
+(
+    const Foam::volVectorField& U
+) const
+{
+    scalar urf = 1;
+
+    if (mesh_.solutionDict().relaxEquation(U.name()))
+    {
+        urf = mesh_.solutionDict().equationRelaxationFactor(U.name());
+    }
+
+    return dimensionedScalar("alphaU", dimless, urf);
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::solutionControl::solutionControl(fvMesh& mesh, const word& algorithmName)
@@ -331,6 +347,9 @@ void Foam::solutionControl::calcTimeConsistentFlux
     // Interpolate original rAU on the faces
     const surfaceScalarField rAUf = fvc::interpolate(rAU);
 
+    // Store previous iteration for the correct handling of under-relaxation
+    phi.storePrevIter();
+
     // Calculate the ordinary part of the flux (H/A)
     phi = (faceU & mesh_.Sf());
 
@@ -471,19 +490,4 @@ const Foam::volScalarField& Foam::solutionControl::aCoeff() const
 }
 
 
-const Foam::dimensionedScalar Foam::solutionControl::relaxFactor
-(
-    const Foam::volVectorField& U
-) const
-{
-    return
-        dimensionedScalar
-        (
-            "alphaU",
-            dimless,
-            mesh_.solutionDict().equationRelaxationFactor(U.name())
-        );
-}
-
-
 // ************************************************************************* //

src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.H

@@ -145,6 +145,17 @@ protected:
         ) const;
 
 
+        // Time and under-relaxation consistency
+
+            //- Get relaxation factor for velocity field. Overriden in
+            //  pimpleControl since different relaxation factor may be used for
+            //  final iteration.
+            virtual const dimensionedScalar relaxFactor
+            (
+                const volVectorField& U
+            ) const;
+
+
 private:
 
     // Private data
@@ -251,14 +262,6 @@ public:
             //- Const access to aCoeff (needed for pressure equation)
             const volScalarField& aCoeff() const;
 
-            //- Get relaxation factor for velocity field. Overriden in
-            //  pimpleControl since different relaxation factor may be used for
-            //  final iteration.
-            virtual const dimensionedScalar relaxFactor
-            (
-                const volVectorField& U
-            ) const;
-
 
         // Evolution
 

src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C

@@ -660,7 +660,7 @@ CoEulerDdtScheme<Type>::fvcDdtConsistentPhiCorr
     const surfaceScalarField& rAUf
 ) const
 {
-    return (mesh().Sf() & faceU.oldTime())/rAUf*CofrDeltaT();
+    return (mesh().Sf() & faceU.oldTime())*rAUf*CofrDeltaT();
 }
 
 

src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C

@@ -522,7 +522,7 @@ EulerDdtScheme<Type>::fvcDdtConsistentPhiCorr
     const surfaceScalarField& rAUf
 ) const
 {
-    return (mesh().Sf() & faceU.oldTime())/rAUf/mesh().time().deltaT();
+    return (mesh().Sf() & faceU.oldTime())*rAUf/mesh().time().deltaT();
 }
 
 

src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C

@@ -665,7 +665,7 @@ SLTSDdtScheme<Type>::fvcDdtConsistentPhiCorr
     const surfaceScalarField& rAUf
 ) const
 {
-    return (mesh().Sf() & faceU.oldTime())/rAUf*fvc::interpolate(SLrDeltaT());
+    return (mesh().Sf() & faceU.oldTime())*rAUf*fvc::interpolate(SLrDeltaT());
 }
 
 

src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C

@@ -732,8 +732,8 @@ backwardDdtScheme<Type>::fvcDdtConsistentPhiCorr
 
     const scalar rDeltaT = 1.0/deltaT;
 
-    const dimensionedScalar beta("beta", dimTime, coefft0*rDeltaT);
+    const dimensionedScalar beta("beta", dimless/dimTime, coefft0*rDeltaT);
-    const dimensionedScalar gamma("gamma", dimTime, -coefft00*rDeltaT);
+    const dimensionedScalar gamma("gamma", dimless/dimTime, -coefft00*rDeltaT);
 
     return
         rAUf*

src/finiteVolume/finiteVolume/ddtSchemes/steadyInertialDdtScheme/steadyInertialDdtScheme.C

@@ -666,7 +666,7 @@ steadyInertialDdtScheme<Type>::fvcDdtConsistentPhiCorr
     const surfaceScalarField& rAUf
 ) const
 {
-    return (mesh().Sf() & faceU.oldTime())/rAUf*CofrDeltaT();
+    return (mesh().Sf() & faceU.oldTime())*rAUf*CofrDeltaT();
 }