Merge branch 'bugfix/missingMovingMeshTreatmentInConsistency' into CumulativeDevelopment-VukoVukcevic

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

@@ -660,7 +660,35 @@ CoEulerDdtScheme<Type>::fvcDdtConsistentPhiCorr
     const surfaceScalarField& rAUf
 )
 {
-    return (mesh().Sf() & faceU.oldTime())*rAUf*CofrDeltaT();
+    tmp<fluxFieldType> toldTimeFlux =
+        (mesh().Sf() & faceU.oldTime())*rAUf*CofrDeltaT();
+
+    if (mesh().moving())
+    {
+        // Mesh is moving, need to take into account the ratio between old and
+        // current cell volumes
+        volScalarField V0ByV
+        (
+            IOobject
+            (
+                "V0ByV",
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("one", dimless, 1.0),
+            zeroGradientFvPatchScalarField::typeName
+        );
+        V0ByV.internalField() = mesh().V0()/mesh().V();
+        V0ByV.correctBoundaryConditions();
+
+        // Correct the flux with interpolated volume ratio
+        toldTimeFlux() *= fvc::interpolate(V0ByV);
+    }
+
+    return toldTimeFlux;
 }
 
 

src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C

@@ -1201,16 +1201,39 @@ CrankNicolsonDdtScheme<Type>::fvcDdtConsistentPhiCorr
           - offCentre_(faceUDdt0())
         );
     }
+    
+    // Calculate old time flux
+    fluxFieldType oldTimeFlux =
+        rAUf*rDtCoef_(faceUDdt0)*(mesh().Sf() & faceU.oldTime());
+
+    if (mesh().moving())
+    {
+        // Mesh is moving, need to take into account the ratio between old and
+        // current cell volumes
+        volScalarField V0ByV
+        (
+            IOobject
+            (
+                "V0ByV",
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("one", dimless, 1.0),
+            zeroGradientFvPatchScalarField::typeName
+        );
+        V0ByV.internalField() = mesh().V0()/mesh().V();
+        V0ByV.correctBoundaryConditions();
+
+        // Correct the flux with interpolated volume ratio
+        oldTimeFlux *= fvc::interpolate(V0ByV);
+    }
 
     return
-        rAUf*
-        (
-            mesh().Sf()
-          & (
-                rDtCoef_(faceUDdt0)*faceU.oldTime()
-              + offCentre_(faceUDdt0())
-            )
-        );
+        oldTimeFlux
+      + rAUf*rDtCoef_(faceUDdt0)*(mesh().Sf() & offCentre_(faceUDdt0()));
 }
 
 

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

@@ -522,7 +522,35 @@ EulerDdtScheme<Type>::fvcDdtConsistentPhiCorr
     const surfaceScalarField& rAUf
 )
 {
-    return (mesh().Sf() & faceU.oldTime())*rAUf/mesh().time().deltaT();
+    tmp<fluxFieldType> toldTimeFlux =
+        (mesh().Sf() & faceU.oldTime())*rAUf/mesh().time().deltaT();
+
+    if (mesh().moving())
+    {
+        // Mesh is moving, need to take into account the ratio between old and
+        // current cell volumes
+        volScalarField V0ByV
+        (
+            IOobject
+            (
+                "V0ByV",
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("one", dimless, 1.0),
+            zeroGradientFvPatchScalarField::typeName
+        );
+        V0ByV.internalField() = mesh().V0()/mesh().V();
+        V0ByV.correctBoundaryConditions();
+
+        // Correct the flux with interpolated volume ratio
+        toldTimeFlux() *= fvc::interpolate(V0ByV);
+    }
+
+    return toldTimeFlux;
 }
 
 

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

@@ -665,7 +665,35 @@ SLTSDdtScheme<Type>::fvcDdtConsistentPhiCorr
     const surfaceScalarField& rAUf
 )
 {
-    return (mesh().Sf() & faceU.oldTime())*rAUf*fvc::interpolate(SLrDeltaT());
+    tmp<fluxFieldType> toldTimeFlux =
+        (mesh().Sf() & faceU.oldTime())*rAUf*fvc::interpolate(SLrDeltaT());
+
+    if (mesh().moving())
+    {
+        // Mesh is moving, need to take into account the ratio between old and
+        // current cell volumes
+        volScalarField V0ByV
+        (
+            IOobject
+            (
+                "V0ByV",
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("one", dimless, 1.0),
+            zeroGradientFvPatchScalarField::typeName
+        );
+        V0ByV.internalField() = mesh().V0()/mesh().V();
+        V0ByV.correctBoundaryConditions();
+
+        // Correct the flux with interpolated volume ratio
+        toldTimeFlux() *= fvc::interpolate(V0ByV);
+    }
+
+    return toldTimeFlux;
 }
 
 

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

@@ -732,18 +732,66 @@ backwardDdtScheme<Type>::fvcDdtConsistentPhiCorr
 
     const scalar rDeltaT = 1.0/deltaT;
 
+    // Note: minus sign in gamma coefficient so we can simply add the fluxes
+    // together at the end
     const dimensionedScalar beta("beta", dimless/dimTime, coefft0*rDeltaT);
     const dimensionedScalar gamma("gamma", dimless/dimTime, -coefft00*rDeltaT);
 
-    return
-        rAUf*
+    // Calculate old and old-old flux contributions
+    fluxFieldType oldTimeFlux =
+        beta*rAUf*(mesh().Sf() & faceU.oldTime());
+    fluxFieldType oldOldTimeFlux =
+        gamma*rAUf*(mesh().Sf() & faceU.oldTime().oldTime());
+
+    if (mesh().moving())
+    {
+        // Mesh is moving, need to take into account the ratio between old and
+        // current cell volumes for old flux contribution
+        volScalarField V0ByV
         (
-            mesh().Sf()
-          & (
-                beta*faceU.oldTime()
-              + gamma*faceU.oldTime().oldTime()
-            )
+            IOobject
+            (
+                "V0ByV",
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("one", dimless, 1.0),
+            zeroGradientFvPatchScalarField::typeName
         );
+        V0ByV.internalField() = mesh().V0()/mesh().V();
+        V0ByV.correctBoundaryConditions();
+
+        // Correct old time flux contribution
+        oldTimeFlux *= fvc::interpolate(V0ByV);
+
+
+        // Also need to take into account the ratio between old-old and current
+        // cell volumes for old-old time flux contribution
+        volScalarField V00ByV
+        (
+            IOobject
+            (
+                "V00ByV",
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("one", dimless, 1.0),
+            zeroGradientFvPatchScalarField::typeName
+        );
+        V00ByV.internalField() = mesh().V00()/mesh().V();
+        V00ByV.correctBoundaryConditions();
+
+        // Correct old-old time flux contribution
+        oldOldTimeFlux *= fvc::interpolate(V00ByV);
+    }
+
+    return oldTimeFlux + oldOldTimeFlux;
 }
 
 

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

@@ -667,7 +667,35 @@ steadyInertialDdtScheme<Type>::fvcDdtConsistentPhiCorr
     const surfaceScalarField& rAUf
 )
 {
-    return (mesh().Sf() & faceU.oldTime())*rAUf*CofrDeltaT();
+    tmp<fluxFieldType> toldTimeFlux =
+        (mesh().Sf() & faceU.oldTime())*rAUf*CofrDeltaT();
+
+    if (mesh().moving())
+    {
+        // Mesh is moving, need to take into account the ratio between old and
+        // current cell volumes
+        volScalarField V0ByV
+        (
+            IOobject
+            (
+                "V0ByV",
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("one", dimless, 1.0),
+            zeroGradientFvPatchScalarField::typeName
+        );
+        V0ByV.internalField() = mesh().V0()/mesh().V();
+        V0ByV.correctBoundaryConditions();
+
+        // Correct the flux with interpolated volume ratio
+        toldTimeFlux() *= fvc::interpolate(V0ByV);
+    }
+
+    return toldTimeFlux;
 }
 
 

validationAndVerificationSuite/verification/incompressible/icoFoam/consistencyTest/cavity/Allrun

@@ -18,4 +18,4 @@ done
 
 # Print out the converged pressure for all time steps for visual check whether
 # the solution does not depend on the time step
-tail -n 1 */probes/0/p
+tail -n 1 */postProcessing/probes/0/p

validationAndVerificationSuite/verification/incompressible/icoFoam/consistencyTest/cavity/explanationAndReference

@@ -5,6 +5,6 @@ simulations are performed with four different time steps spanning four orders of
 magnitude: 0.01, 0.001, 0.0001 and 0.00001 s. Tolerances for all equations are
 very small, yielding extremely small (O(1e-11)) differences in converged
 pressure field. If the differences are larger - it means that the converged
-solution is not independent to time step size.
+solution depends on the time step size.
 
 Author: Vuko Vukcevic, vuko.vukcevic@fsb.hr

validationAndVerificationSuite/verification/incompressible/simpleFoam/consistencyTest/cavity/Allrun

@@ -18,4 +18,4 @@ done
 
 # Print out the converged pressure for all relaxation factors for visual check
 # whether the solution does not depend on the under-relaxation factors
-tail -n 1 */probes/0/p
+tail -n 1 */postProcessing/probes/0/p

validationAndVerificationSuite/verification/incompressible/simpleFoam/consistencyTest/cavity/explanationAndReference

@@ -9,6 +9,6 @@ simulations are performed with five different under-relaxation pairs:
     5. alphaU = 0.4, alphap = 0.6
 Tolerances for all equations are very small, yielding extremely small (O(1e-11))
 differences in converged pressure field. If the differences are larger - it
-means that the converged solution is not independent to relaxation factors.
+means that the converged solution is dependend on relaxation factors.
 
 Author: Vuko Vukcevic, vuko.vukcevic@fsb.hr