Bugfix: interIsoFoam did not work properly when both nAlphaSubCycles > 1 and...

Bugfix: interIsoFoam did not work properly when both nAlphaSubCycles > 1 and nOuterCorrectors > 1. This has now ben fixed. Credits to Henning Scheufler from DLR Bremen.

applications/solvers/multiphase/interIsoFoam/alphaEqn.H

 {
-    if (pimple.nCorrPIMPLE() > 1)
-    {
-        // If nOuterCorrectors > 1 then for all but the first loop the advection
-        // of alpha is done using an average, 0.5*phi+0.5*phiNew where phi is
-        // the flux at the beginning of the time step and phiNew is the flux
-        // estimate at the end of the time step from the previous outer
-        // iteration. Similarly we use 0.5*U + 0.5*UNew in later iterations.
-        if (pimple.firstIter())
-        {
-            // To recalculate the alpha1 update in subsequent iterations, we
-            // must store its current value before overwriting with the new
-            // value
-            alpha1.storePrevIter();
-            // Storing initial phi and U for use in later outer iterations.
-            phi.storePrevIter();
-            U.storePrevIter();
-        }
-        else
-        {
-            // Resetting alpha1 to value before advection in first PIMPLE
-            // iteration.
-            alpha1 = alpha1.prevIter();
-
-            // Setting U and phi with which to advect interface.
-            U = 0.5*U.prevIter() + 0.5*U;
-            phi = 0.5*phi.prevIter() + 0.5*phi;
-        }
-    }
-
     // Updating alpha1
     advector.advect();
     #include "rhofs.H"
     rhoPhi = advector.getRhoPhi(rho1f, rho2f);
 
-    if (!pimple.firstIter())
-    {
-        // Resetting U and phi to value at latest iteration.
-        U = 2.0*U - U.prevIter();
-        phi = 2.0*phi - phi.prevIter();
-    }
-
     alpha2 = 1.0 - alpha1;
     mixture.correct();
-
 }
 
 Info<< "Phase-1 volume fraction = "

applications/solvers/multiphase/interIsoFoam/alphaEqnSubCycle.H

+if (pimple.nCorrPIMPLE() > 1)
+{
+    // If nOuterCorrectors > 1 then for all but the first loop the advection
+    // of alpha is done using an average, 0.5*phi+0.5*phiNew where phi is
+    // the flux at the beginning of the time step and phiNew is the flux
+    // estimate at the end of the time step from the previous outer
+    // iteration. Similarly we use 0.5*U + 0.5*UNew in later iterations.
+    if (pimple.firstIter())
+    {
+        // To recalculate the alpha1 update in subsequent iterations, we
+        // must store its current value before overwriting with the new
+        // value
+        alpha1.storePrevIter();
+        // Storing initial phi and U for use in later outer iterations.
+        phi.storePrevIter();
+        U.storePrevIter();
+    }
+    else
+    {
+        // Resetting alpha1 to value before advection in first PIMPLE
+        // iteration.
+        alpha1 = alpha1.prevIter();
+
+        // Temporarily setting U and phi with which to advect interface.
+        U = 0.5*U.prevIter() + 0.5*U;
+        phi = 0.5*phi.prevIter() + 0.5*phi;
+    }
+}
+
 if (nAlphaSubCycles > 1)
 {
     dimensionedScalar totalDeltaT = runTime.deltaT();
@@ -30,4 +59,11 @@ else
     #include "alphaEqn.H"
 }
 
+if (!pimple.firstIter())
+{
+    // Resetting U and phi to value at latest iteration.
+    U = 2.0*U - U.prevIter();
+    phi = 2.0*phi - phi.prevIter();
+}
+
 rho == alpha1*rho1 + alpha2*rho2;