VoF Solvers: Relocate the correction, sub-cycling and compressions controls from PIMPLE to the alpha1 sub-dict

MULES: Add support for explicit correction
interPhaseChangeFoam: Add option for explicit MULES or as correction to an upwind solution
Deprecate implicit form of MULES

applications/solvers/multiphase/compressibleInterFoam/alphaEqnsSubCycle.H

@@ -1,7 +1,5 @@
 {
-    label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr")));
+    #include "alphaControls.H"
-
-    label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles")));
 
     surfaceScalarField phic(mag(phi/mesh.magSf()));
     phic = min(interface.cAlpha()*phic, max(phic));

applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/readControls.H

@@ -1,17 +1,5 @@
     #include "readTimeControls.H"
 
-    label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr")));
-
-    label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles")));
-
-    if (nAlphaSubCycles > 1 && pimple.nCorrPIMPLE() != 1)
-    {
-        FatalErrorIn(args.executable())
-            << "Sub-cycling alpha is only allowed for PISO, "
-               "i.e. when the number of outer-correctors = 1"
-            << exit(FatalError);
-    }
-
     bool correctPhi =
         pimple.dict().lookupOrDefault<Switch>("correctPhi", true);
 

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C

@@ -56,7 +56,7 @@ int main(int argc, char *argv[])
 
     pimpleControl pimple(mesh);
 
-    #include "readControls.H"
+    #include "readTimeControls.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
     #include "CourantNo.H"
@@ -68,7 +68,7 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readControls.H"
+        #include "readTimeControls.H"
         #include "CourantNo.H"
         #include "setDeltaT.H"
 

applications/solvers/multiphase/compressibleInterFoam/pEqn.H

@@ -85,8 +85,8 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            //p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
+            p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
-            //p_rgh = p - (alpha1*rho1 + alpha2*rho2)*gh;
+            p_rgh = p - (alpha1*rho1 + alpha2*rho2)*gh;
 
             dgdt =
             (
@@ -102,7 +102,7 @@
         }
     }
 
-    p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
+    // p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
 
     // Update densities from change in p_rgh
     rho1 += psi1*(p_rgh - p_rgh_0);

applications/solvers/multiphase/compressibleInterFoam/readControls.H

@@ -1,13 +0,0 @@
-    #include "readTimeControls.H"
-
-    label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr")));
-
-    label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles")));
-
-    if (nAlphaSubCycles > 1 && pimple.nCorrPIMPLE() != 1)
-    {
-        FatalErrorIn(args.executable())
-            << "Sub-cycling alpha is only allowed for PISO operation, "
-               "i.e. when the number of outer-correctors = 1"
-            << exit(FatalError);
-    }

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H

@@ -1,4 +1,2 @@
     #include "readTimeControls.H"
-
+    #include "alphaControls.H"
-    int nAlphaCorr(readInt(pimple.dict().lookup("nAlphaCorr")));
-    int nAlphaSubCycles(readInt(pimple.dict().lookup("nAlphaSubCycles")));

applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H

@@ -1,5 +1,4 @@
-label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr")));
+#include "alphaControls.H"
-label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles")));
 
 if (nAlphaSubCycles > 1)
 {

applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H

@@ -130,10 +130,7 @@
             << ", " << gMax(1/rDeltaT.internalField()) << endl;
     }
 
-    label nAlphaSubCycles
+    #include "alphaControls.H"
-    (
-        readLabel(pimpleDict.lookup("nAlphaSubCycles"))
-    );
 
     rSubDeltaT = rDeltaT*nAlphaSubCycles;
 }

applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H

@@ -1,5 +1,4 @@
-label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr")));
+#include "alphaControls.H"
-label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles")));
 
 if (nAlphaSubCycles > 1)
 {

applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqnsSubCycle.H

@@ -1,6 +1,4 @@
-const dictionary& pimpleDict = pimple.dict();
+#include "alphaControls.H"
-label nAlphaCorr(readLabel(pimpleDict.lookup("nAlphaCorr")));
-label nAlphaSubCycles(readLabel(pimpleDict.lookup("nAlphaSubCycles")));
 
 if (nAlphaSubCycles > 1)
 {

applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H

@@ -4,9 +4,41 @@
 
     surfaceScalarField phir("phir", phic*interface.nHatf());
 
-    for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
+    Pair<tmp<volScalarField> > vDotAlphal =
+        twoPhaseProperties->vDotAlphal();
+    const volScalarField& vDotcAlphal = vDotAlphal[0]();
+    const volScalarField& vDotvAlphal = vDotAlphal[1]();
+    const volScalarField vDotvmcAlphal(vDotvAlphal - vDotcAlphal);
+
+    tmp<surfaceScalarField> tphiAlpha;
+
+    if (MULESCorr)
     {
-        surfaceScalarField phiAlpha
+        fvScalarMatrix alpha1Eqn
+        (
+            fvm::ddt(alpha1)
+          + fvm::div(phi, alpha1, "UD") - fvm::Sp(divU, alpha1)
+         ==
+            fvm::Sp(vDotvmcAlphal, alpha1)
+          + vDotcAlphal
+        );
+
+        alpha1Eqn.solve();
+
+        Info<< "Phase-1 volume fraction = "
+            << alpha1.weightedAverage(mesh.Vsc()).value()
+            << "  Min(alpha1) = " << min(alpha1).value()
+            << "  Max(alpha1) = " << max(alpha1).value()
+            << endl;
+
+        tphiAlpha = alpha1Eqn.flux();
+    }
+
+    volScalarField alpha10("alpha10", alpha1);
+
+    for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
+    {
+        tmp<surfaceScalarField> tphiAlphaCorr
         (
             fvc::flux
             (
@@ -22,66 +54,52 @@
             )
         );
 
-        Pair<tmp<volScalarField> > vDotAlphal =
+        if (MULESCorr)
-            twoPhaseProperties->vDotAlphal();
+        {
-        const volScalarField& vDotcAlphal = vDotAlphal[0]();
+            tphiAlphaCorr() -= tphiAlpha();
-        const volScalarField& vDotvAlphal = vDotAlphal[1]();
 
-        volScalarField Sp
+
-        (
+            volScalarField alpha100("alpha100", alpha10);
-            IOobject
+            alpha10 = alpha1;
+
+            MULES::correct
             (
-                "Sp",
+                geometricOneField(),
-                runTime.timeName(),
+                alpha1,
-                mesh
+                tphiAlphaCorr(),
-            ),
+                vDotvmcAlphal,
-            vDotvAlphal - vDotcAlphal
+                (
-        );
+                    divU*(alpha10 - alpha100)
+                  - vDotvmcAlphal*alpha10
+                )(),
+                1,
+                0
+            );
 
-        volScalarField Su
+            tphiAlpha() += tphiAlphaCorr();
-        (
+        }
-            IOobject
+        else
+        {
+            MULES::explicitSolve
             (
-                "Su",
+                geometricOneField(),
-                runTime.timeName(),
+                alpha1,
-                mesh
+                phi,
-            ),
+                tphiAlphaCorr(),
-            // Divergence term is handled explicitly to be
+                vDotvmcAlphal,
-            // consistent with the explicit transport solution
+                (divU*alpha1 + vDotcAlphal)(),
-            divU*alpha1
+                1,
-          + vDotcAlphal
+                0
-        );
+            );
 
-        //MULES::explicitSolve
+            tphiAlpha = tphiAlphaCorr;
-        //(
+        }
-        //    geometricOneField(),
-        //    alpha1,
-        //    phi,
-        //    phiAlpha,
-        //    Sp,
-        //    Su,
-        //    1,
-        //    0
-        //);
-
-        MULES::implicitSolve
-        (
-            geometricOneField(),
-            alpha1,
-            phi,
-            phiAlpha,
-            Sp,
-            Su,
-            1,
-            0
-        );
 
         alpha2 = 1.0 - alpha1;
-        rhoPhi +=
-            (runTime.deltaT()/totalDeltaT)
-           *(phiAlpha*(rho1 - rho2) + phi*rho2);
     }
 
+    rhoPhi = tphiAlpha()*(rho1 - rho2) + phi*rho2;
+
     Info<< "Liquid phase volume fraction = "
         << alpha1.weightedAverage(mesh.V()).value()
         << "  Min(alpha1) = " << min(alpha1).value()

applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H

@@ -11,21 +11,18 @@ surfaceScalarField rhoPhi
 );
 
 {
-    const dictionary& pimpleDict = pimple.dict();
+    #include "alphaControls.H"
-
-    label nAlphaCorr(readLabel(pimpleDict.lookup("nAlphaCorr")));
-
-    label nAlphaSubCycles(readLabel(pimpleDict.lookup("nAlphaSubCycles")));
 
     surfaceScalarField phic(mag(phi/mesh.magSf()));
     phic = min(interface.cAlpha()*phic, max(phic));
 
     volScalarField divU(fvc::div(phi));
 
-    dimensionedScalar totalDeltaT = runTime.deltaT();
-
     if (nAlphaSubCycles > 1)
     {
+        dimensionedScalar totalDeltaT = runTime.deltaT();
+        surfaceScalarField rhoPhiSum("rhoPhiSum", rhoPhi);
+
         for
         (
             subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
@@ -33,7 +30,10 @@ surfaceScalarField rhoPhi
         )
         {
             #include "alphaEqn.H"
+            rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
         }
+
+        rhoPhi = rhoPhiSum;
     }
     else
     {

applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C

@@ -41,7 +41,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "IMULES.H"
+#include "MULES.H"
 #include "subCycle.H"
 #include "interfaceProperties.H"
 #include "phaseChangeTwoPhaseMixture.H"

applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C

@@ -812,9 +812,8 @@ void Foam::multiphaseSystem::solve()
 
     const Time& runTime = mesh_.time();
 
-    const dictionary& pimpleDict = mesh_.solutionDict().subDict("PIMPLE");
+    const dictionary& alphaControls = mesh_.solverDict(phases_.first().name());
-
+    label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
-    label nAlphaSubCycles(readLabel(pimpleDict.lookup("nAlphaSubCycles")));
 
     if (nAlphaSubCycles > 1)
     {

applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C

@@ -249,15 +249,12 @@ void Foam::multiphaseMixture::solve()
 
     const Time& runTime = mesh_.time();
 
-    const dictionary& pimpleDict = mesh_.solutionDict().subDict("PIMPLE");
-
-    label nAlphaSubCycles(readLabel(pimpleDict.lookup("nAlphaSubCycles")));
-
-    scalar cAlpha(readScalar(pimpleDict.lookup("cAlpha")));
-
-
     volScalarField& alpha = phases_.first();
 
+    const dictionary& alphaControls = mesh_.solverDict(alpha.name());
+    label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
+    scalar cAlpha(readScalar(alphaControls.lookup("cAlpha")));
+
     if (nAlphaSubCycles > 1)
     {
         surfaceScalarField rhoPhiSum(0.0*rhoPhi_);

applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqnSubCycle.H

@@ -1,5 +1,4 @@
-label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr")));
+#include "alphaControls.H"
-label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles")));
 
 if (nAlphaSubCycles > 1)
 {

applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H

@@ -1,5 +1,4 @@
     #include "readTimeControls.H"
+    #include "alphaControls.H"
 
-    int nAlphaCorr(readInt(pimple.dict().lookup("nAlphaCorr")));
-    int nAlphaSubCycles(readInt(pimple.dict().lookup("nAlphaSubCycles")));
     Switch correctAlpha(pimple.dict().lookup("correctAlpha"));

src/finiteVolume/fvMatrices/solvers/MULES/MULES.C

@@ -48,6 +48,25 @@ void Foam::MULES::explicitSolve
 }
 
 
+void Foam::MULES::correct
+(
+    volScalarField& psi,
+    surfaceScalarField& phiPsiCorr,
+    const scalar psiMax,
+    const scalar psiMin
+)
+{
+    correct
+    (
+        geometricOneField(),
+        psi,
+        phiPsiCorr,
+        zeroField(), zeroField(),
+        psiMax, psiMin
+    );
+}
+
+
 void Foam::MULES::limitSum(UPtrList<scalarField>& phiPsiCorrs)
 {
     forAll(phiPsiCorrs[0], facei)

src/finiteVolume/fvMatrices/solvers/MULES/MULES.H

@@ -92,6 +92,38 @@ void explicitSolve
     const scalar psiMin
 );
 
+
+template<class RhoType, class SpType, class SuType>
+void correct
+(
+    const RhoType& rho,
+    volScalarField& psi,
+    const surfaceScalarField& phiCorr,
+    const SpType& Sp,
+    const SuType& Su
+);
+
+template<class RhoType, class SpType, class SuType>
+void correct
+(
+    const RhoType& rho,
+    volScalarField& psi,
+    surfaceScalarField& phiCorr,
+    const SpType& Sp,
+    const SuType& Su,
+    const scalar psiMax,
+    const scalar psiMin
+);
+
+void correct
+(
+    volScalarField& psi,
+    surfaceScalarField& phiCorr,
+    const scalar psiMax,
+    const scalar psiMin
+);
+
+
 template<class RhoType, class SpType, class SuType>
 void limiter
 (
@@ -122,6 +154,35 @@ void limit
     const bool returnCorr
 );
 
+
+template<class RhoType, class SpType, class SuType>
+void limiterCorr
+(
+    scalarField& allLambda,
+    const RhoType& rho,
+    const volScalarField& psi,
+    const surfaceScalarField& phiCorr,
+    const SpType& Sp,
+    const SuType& Su,
+    const scalar psiMax,
+    const scalar psiMin,
+    const label nLimiterIter
+);
+
+template<class RhoType, class SpType, class SuType>
+void limitCorr
+(
+    const RhoType& rho,
+    const volScalarField& psi,
+    surfaceScalarField& phiCorr,
+    const SpType& Sp,
+    const SuType& Su,
+    const scalar psiMax,
+    const scalar psiMin,
+    const label nLimiterIter
+);
+
+
 void limitSum(UPtrList<scalarField>& phiPsiCorrs);
 
 template<class SurfaceScalarFieldList>

src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C

@@ -96,6 +96,74 @@ void Foam::MULES::explicitSolve
 }
 
 
+template<class RhoType, class SpType, class SuType>
+void Foam::MULES::correct
+(
+    const RhoType& rho,
+    volScalarField& psi,
+    const surfaceScalarField& phiCorr,
+    const SpType& Sp,
+    const SuType& Su
+)
+{
+    Info<< "MULES: Correcting " << psi.name() << endl;
+
+    const fvMesh& mesh = psi.mesh();
+
+    const scalar deltaT = mesh.time().deltaTValue();
+
+    scalarField psiIf(psi.size(), 0);
+    fvc::surfaceIntegrate(psiIf, phiCorr);
+
+    if (mesh.moving())
+    {
+        psi.internalField() =
+        (
+            rho.field()*psi.internalField()/deltaT
+          + Su.field()
+          - psiIf
+        )/(rho.field()/deltaT - Sp.field());
+    }
+    else
+    {
+        psi.internalField() =
+        (
+            rho.field()*psi.internalField()/deltaT
+          + Su.field()
+          - psiIf
+        )/(rho.field()/deltaT - Sp.field());
+    }
+
+    psi.correctBoundaryConditions();
+}
+
+
+template<class RhoType, class SpType, class SuType>
+void Foam::MULES::correct
+(
+    const RhoType& rho,
+    volScalarField& psi,
+    surfaceScalarField& phiCorr,
+    const SpType& Sp,
+    const SuType& Su,
+    const scalar psiMax,
+    const scalar psiMin
+)
+{
+    const fvMesh& mesh = psi.mesh();
+
+    const dictionary& MULEScontrols = mesh.solverDict(psi.name());
+
+    label nLimiterIter
+    (
+        readLabel(MULEScontrols.lookup("nLimiterIter"))
+    );
+
+    limitCorr(rho, psi, phiCorr, Sp, Su, psiMax, psiMin, nLimiterIter);
+    correct(rho, psi, phiCorr, Sp, Su);
+}
+
+
 template<class RhoType, class SpType, class SuType>
 void Foam::MULES::limiter
 (
@@ -511,6 +579,364 @@ void Foam::MULES::limit
 }
 
 
+template<class RhoType, class SpType, class SuType>
+void Foam::MULES::limiterCorr
+(
+    scalarField& allLambda,
+    const RhoType& rho,
+    const volScalarField& psi,
+    const surfaceScalarField& phiCorr,
+    const SpType& Sp,
+    const SuType& Su,
+    const scalar psiMax,
+    const scalar psiMin,
+    const label nLimiterIter
+)
+{
+    const scalarField& psiIf = psi;
+    const volScalarField::GeometricBoundaryField& psiBf = psi.boundaryField();
+
+    const fvMesh& mesh = psi.mesh();
+
+    const labelUList& owner = mesh.owner();
+    const labelUList& neighb = mesh.neighbour();
+    tmp<volScalarField::DimensionedInternalField> tVsc = mesh.Vsc();
+    const scalarField& V = tVsc();
+    const scalar deltaT = mesh.time().deltaTValue();
+
+    const scalarField& phiCorrIf = phiCorr;
+    const surfaceScalarField::GeometricBoundaryField& phiCorrBf =
+        phiCorr.boundaryField();
+
+    slicedSurfaceScalarField lambda
+    (
+        IOobject
+        (
+            "lambda",
+            mesh.time().timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        mesh,
+        dimless,
+        allLambda,
+        false   // Use slices for the couples
+    );
+
+    scalarField& lambdaIf = lambda;
+    surfaceScalarField::GeometricBoundaryField& lambdaBf =
+        lambda.boundaryField();
+
+    scalarField psiMaxn(psiIf.size(), psiMin);
+    scalarField psiMinn(psiIf.size(), psiMax);
+
+    scalarField sumPhip(psiIf.size(), VSMALL);
+    scalarField mSumPhim(psiIf.size(), VSMALL);
+
+    forAll(phiCorrIf, facei)
+    {
+        label own = owner[facei];
+        label nei = neighb[facei];
+
+        psiMaxn[own] = max(psiMaxn[own], psiIf[nei]);
+        psiMinn[own] = min(psiMinn[own], psiIf[nei]);
+
+        psiMaxn[nei] = max(psiMaxn[nei], psiIf[own]);
+        psiMinn[nei] = min(psiMinn[nei], psiIf[own]);
+
+        scalar phiCorrf = phiCorrIf[facei];
+
+        if (phiCorrf > 0.0)
+        {
+            sumPhip[own] += phiCorrf;
+            mSumPhim[nei] += phiCorrf;
+        }
+        else
+        {
+            mSumPhim[own] -= phiCorrf;
+            sumPhip[nei] -= phiCorrf;
+        }
+    }
+
+    forAll(phiCorrBf, patchi)
+    {
+        const fvPatchScalarField& psiPf = psiBf[patchi];
+        const scalarField& phiCorrPf = phiCorrBf[patchi];
+
+        const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
+
+        if (psiPf.coupled())
+        {
+            const scalarField psiPNf(psiPf.patchNeighbourField());
+
+            forAll(phiCorrPf, pFacei)
+            {
+                label pfCelli = pFaceCells[pFacei];
+
+                psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPNf[pFacei]);
+                psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPNf[pFacei]);
+            }
+        }
+        else
+        {
+            forAll(phiCorrPf, pFacei)
+            {
+                label pfCelli = pFaceCells[pFacei];
+
+                psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPf[pFacei]);
+                psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPf[pFacei]);
+            }
+        }
+
+        forAll(phiCorrPf, pFacei)
+        {
+            label pfCelli = pFaceCells[pFacei];
+
+            scalar phiCorrf = phiCorrPf[pFacei];
+
+            if (phiCorrf > 0.0)
+            {
+                sumPhip[pfCelli] += phiCorrf;
+            }
+            else
+            {
+                mSumPhim[pfCelli] -= phiCorrf;
+            }
+        }
+    }
+
+    psiMaxn = min(psiMaxn, psiMax);
+    psiMinn = max(psiMinn, psiMin);
+
+    // scalar smooth = 0.5;
+    // psiMaxn = min((1.0 - smooth)*psiIf + smooth*psiMaxn, psiMax);
+    // psiMinn = max((1.0 - smooth)*psiIf + smooth*psiMinn, psiMin);
+
+    psiMaxn =
+        V
+       *(
+           (rho.field()/deltaT - Sp.field())*psiMaxn
+         - Su.field()
+         - rho.field()*psi.internalField()/deltaT
+        );
+
+    psiMinn =
+        V
+       *(
+           Su.field()
+         - (rho.field()/deltaT - Sp.field())*psiMinn
+         + rho.field()*psi.internalField()/deltaT
+        );
+
+    scalarField sumlPhip(psiIf.size());
+    scalarField mSumlPhim(psiIf.size());
+
+    for (int j=0; j<nLimiterIter; j++)
+    {
+        sumlPhip = 0.0;
+        mSumlPhim = 0.0;
+
+        forAll(lambdaIf, facei)
+        {
+            label own = owner[facei];
+            label nei = neighb[facei];
+
+            scalar lambdaPhiCorrf = lambdaIf[facei]*phiCorrIf[facei];
+
+            if (lambdaPhiCorrf > 0.0)
+            {
+                sumlPhip[own] += lambdaPhiCorrf;
+                mSumlPhim[nei] += lambdaPhiCorrf;
+            }
+            else
+            {
+                mSumlPhim[own] -= lambdaPhiCorrf;
+                sumlPhip[nei] -= lambdaPhiCorrf;
+            }
+        }
+
+        forAll(lambdaBf, patchi)
+        {
+            scalarField& lambdaPf = lambdaBf[patchi];
+            const scalarField& phiCorrfPf = phiCorrBf[patchi];
+
+            const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
+
+            forAll(lambdaPf, pFacei)
+            {
+                label pfCelli = pFaceCells[pFacei];
+
+                scalar lambdaPhiCorrf = lambdaPf[pFacei]*phiCorrfPf[pFacei];
+
+                if (lambdaPhiCorrf > 0.0)
+                {
+                    sumlPhip[pfCelli] += lambdaPhiCorrf;
+                }
+                else
+                {
+                    mSumlPhim[pfCelli] -= lambdaPhiCorrf;
+                }
+            }
+        }
+
+        forAll(sumlPhip, celli)
+        {
+            sumlPhip[celli] =
+                max(min
+                (
+                    (sumlPhip[celli] + psiMaxn[celli])/mSumPhim[celli],
+                    1.0), 0.0
+                );
+
+            mSumlPhim[celli] =
+                max(min
+                (
+                    (mSumlPhim[celli] + psiMinn[celli])/sumPhip[celli],
+                    1.0), 0.0
+                );
+        }
+
+        const scalarField& lambdam = sumlPhip;
+        const scalarField& lambdap = mSumlPhim;
+
+        forAll(lambdaIf, facei)
+        {
+            if (phiCorrIf[facei] > 0.0)
+            {
+                lambdaIf[facei] = min
+                (
+                    lambdaIf[facei],
+                    min(lambdap[owner[facei]], lambdam[neighb[facei]])
+                );
+            }
+            else
+            {
+                lambdaIf[facei] = min
+                (
+                    lambdaIf[facei],
+                    min(lambdam[owner[facei]], lambdap[neighb[facei]])
+                );
+            }
+        }
+
+
+        forAll(lambdaBf, patchi)
+        {
+            fvsPatchScalarField& lambdaPf = lambdaBf[patchi];
+            const scalarField& phiCorrfPf = phiCorrBf[patchi];
+            const fvPatchScalarField& psiPf = psiBf[patchi];
+
+            if (isA<wedgeFvPatch>(mesh.boundary()[patchi]))
+            {
+                lambdaPf = 0;
+            }
+            else if (psiPf.coupled())
+            {
+                const labelList& pFaceCells =
+                    mesh.boundary()[patchi].faceCells();
+
+                forAll(lambdaPf, pFacei)
+                {
+                    label pfCelli = pFaceCells[pFacei];
+
+                    if (phiCorrfPf[pFacei] > 0.0)
+                    {
+                        lambdaPf[pFacei] =
+                            min(lambdaPf[pFacei], lambdap[pfCelli]);
+                    }
+                    else
+                    {
+                        lambdaPf[pFacei] =
+                            min(lambdaPf[pFacei], lambdam[pfCelli]);
+                    }
+                }
+            }
+            else
+            {
+                const labelList& pFaceCells =
+                    mesh.boundary()[patchi].faceCells();
+                const scalarField& phiCorrPf = phiCorrBf[patchi];
+
+                forAll(lambdaPf, pFacei)
+                {
+                    // Limit outlet faces only
+                    if (phiCorrPf[pFacei] > 0)
+                    {
+                        label pfCelli = pFaceCells[pFacei];
+
+                        if (phiCorrfPf[pFacei] > 0.0)
+                        {
+                            lambdaPf[pFacei] =
+                                min(lambdaPf[pFacei], lambdap[pfCelli]);
+                        }
+                        else
+                        {
+                            lambdaPf[pFacei] =
+                                min(lambdaPf[pFacei], lambdam[pfCelli]);
+                        }
+                    }
+                }
+            }
+        }
+
+        syncTools::syncFaceList(mesh, allLambda, minEqOp<scalar>());
+    }
+}
+
+
+template<class RhoType, class SpType, class SuType>
+void Foam::MULES::limitCorr
+(
+    const RhoType& rho,
+    const volScalarField& psi,
+    surfaceScalarField& phiCorr,
+    const SpType& Sp,
+    const SuType& Su,
+    const scalar psiMax,
+    const scalar psiMin,
+    const label nLimiterIter
+)
+{
+    const fvMesh& mesh = psi.mesh();
+
+    scalarField allLambda(mesh.nFaces(), 1.0);
+
+    slicedSurfaceScalarField lambda
+    (
+        IOobject
+        (
+            "lambda",
+            mesh.time().timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        mesh,
+        dimless,
+        allLambda,
+        false   // Use slices for the couples
+    );
+
+    limiterCorr
+    (
+        allLambda,
+        rho,
+        psi,
+        phiCorr,
+        Sp,
+        Su,
+        psiMax,
+        psiMin,
+        nLimiterIter
+    );
+
+    phiCorr *= lambda;
+}
+
+
 template<class SurfaceScalarFieldList>
 void Foam::MULES::limitSum(SurfaceScalarFieldList& phiPsiCorrs)
 {

src/transportModels/interfaceProperties/interfaceProperties.C

@@ -159,7 +159,7 @@ Foam::interfaceProperties::interfaceProperties
     (
         readScalar
         (
-            alpha1.mesh().solutionDict().subDict("PIMPLE").lookup("cAlpha")
+            alpha1.mesh().solverDict(alpha1.name()).lookup("cAlpha")
         )
     ),
     sigma_(dict.lookup("sigma")),

tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/controlDict

@@ -14,37 +14,37 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-application                interPhaseChangeFoam;
+application             interPhaseChangeFoam;
 
-startFrom                  latestTime;
+startFrom               latestTime;
 
-startTime                  0;
+startTime               0;
 
-stopAt                     endTime;
+stopAt                  endTime;
 
-endTime                    0.05;
+endTime                 0.05;
 
-deltaT                     1e-8;
+deltaT                  1e-8;
 
-writeControl               adjustableRunTime;
+writeControl            adjustableRunTime;
 
-writeInterval              0.001;
+writeInterval           0.001;
 
-purgeWrite                 0;
+purgeWrite              0;
 
-writeFormat                ascii;
+writeFormat             ascii;
 
-writePrecision             6;
+writePrecision          6;
 
-writeCompression           uncompressed;
+writeCompression        uncompressed;
 
-timeFormat                 general;
+timeFormat              general;
 
-runTimeModifiable          yes;
+runTimeModifiable       yes;
 
-adjustTimeStep             on;
+adjustTimeStep          on;
 
-maxCo                      1.0;
+maxCo                   5;
 
 
 // ************************************************************************* //

tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSchemes

@@ -32,6 +32,7 @@ divSchemes
     div(phi,k)           Gauss linearUpwind grad(k);
     div(phi,alpha)       Gauss vanLeer;
     div(phirb,alpha)     Gauss interfaceCompression;
+    UD                   Gauss upwind;
     div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
@@ -47,10 +48,7 @@ laplacianSchemes
 
 snGradSchemes
 {
-    default              none;
+    default              limited corrected 0.5;
-    snGrad(pd)           limited corrected 0.5;
-    snGrad(rho)          limited corrected 0.5;
-    snGrad(alpha1)       limited corrected 0.5;
 }
 
 fluxRequired

tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSolution

@@ -18,22 +18,24 @@ solvers
 {
     alpha1
     {
-        maxUnboundedness 1e-5;
+        cAlpha          0;
-        CoCoeff          2;
+        nAlphaCorr      2;
-        maxIter          5;
+        nAlphaSubCycles 1;
-        nLimiterIter     2;
 
-        solver           PBiCG;
+        MULESCorr       yes;
-        preconditioner   DILU;
+        nLimiterIter    5;
-        tolerance        1e-12;
+
-        relTol           0.1;
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-8;
+        relTol          0;
     };
 
-    U
+    "U.*"
     {
         solver           PBiCG;
         preconditioner   DILU;
-        tolerance        1e-06;
+        tolerance        1e-6;
         relTol           0;
     };
 
@@ -96,10 +98,6 @@ PIMPLE
     nOuterCorrectors           1;
     nCorrectors                3;
     nNonOrthogonalCorrectors   0;
-
-    cAlpha                     0;
-    nAlphaCorr                 1;
-    nAlphaSubCycles            1;
 }