CMULES: Corrected the handling of outflow faces

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

@@ -30,6 +30,7 @@ License
 void Foam::MULES::correct
 (
     volScalarField& psi,
+    const surfaceScalarField& phi,
     surfaceScalarField& phiPsiCorr,
     const scalar psiMax,
     const scalar psiMin
@@ -39,6 +40,7 @@ void Foam::MULES::correct
     (
         geometricOneField(),
         psi,
+        phi,
         phiPsiCorr,
         zeroField(), zeroField(),
         psiMax, psiMin
@@ -49,6 +51,7 @@ void Foam::MULES::correct
 void Foam::MULES::LTScorrect
 (
     volScalarField& psi,
+    const surfaceScalarField& phi,
     surfaceScalarField& phiPsiCorr,
     const scalar psiMax,
     const scalar psiMin
@@ -58,6 +61,7 @@ void Foam::MULES::LTScorrect
     (
         geometricOneField(),
         psi,
+        phi,
         phiPsiCorr,
         zeroField(), zeroField(),
         psiMax, psiMin

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

@@ -66,6 +66,7 @@ void correct
     const RdeltaTType& rDeltaT,
     const RhoType& rho,
     volScalarField& psi,
+    const surfaceScalarField& phi,
     const surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su
@@ -76,6 +77,7 @@ void correct
 (
     const RhoType& rho,
     volScalarField& psi,
+    const surfaceScalarField& phi,
     surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su,
@@ -86,6 +88,7 @@ void correct
 void correct
 (
     volScalarField& psi,
+    const surfaceScalarField& phi,
     surfaceScalarField& phiCorr,
     const scalar psiMax,
     const scalar psiMin
@@ -96,6 +99,7 @@ void LTScorrect
 (
     const RhoType& rho,
     volScalarField& psi,
+    const surfaceScalarField& phi,
     surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su,
@@ -106,6 +110,7 @@ void LTScorrect
 void LTScorrect
 (
     volScalarField& psi,
+    const surfaceScalarField& phi,
     surfaceScalarField& phiCorr,
     const scalar psiMax,
     const scalar psiMin
@@ -119,6 +124,7 @@ void limiterCorr
     const RdeltaTType& rDeltaT,
     const RhoType& rho,
     const volScalarField& psi,
+    const surfaceScalarField& phi,
     const surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su,
@@ -133,6 +139,7 @@ void limitCorr
     const RdeltaTType& rDeltaT,
     const RhoType& rho,
     const volScalarField& psi,
+    const surfaceScalarField& phi,
     surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su,

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

@@ -37,6 +37,7 @@ void Foam::MULES::correct
     const RdeltaTType& rDeltaT,
     const RhoType& rho,
     volScalarField& psi,
+    const surfaceScalarField& phi,
     const surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su
@@ -77,6 +78,7 @@ void Foam::MULES::correct
 (
     const RhoType& rho,
     volScalarField& psi,
+    const surfaceScalarField& phi,
     surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su,
@@ -94,8 +96,18 @@ void Foam::MULES::correct
         readLabel(MULEScontrols.lookup("nLimiterIter"))
     );
 
-    limitCorr(rDeltaT, rho, psi, phiCorr, Sp, Su, psiMax, psiMin, nLimiterIter);
-    correct(rDeltaT, rho, psi, phiCorr, Sp, Su);
+    limitCorr
+    (
+        rDeltaT,
+        rho,
+        psi,
+        phi,
+        phiCorr,
+        Sp, Su,
+        psiMax, psiMin,
+        nLimiterIter
+    );
+    correct(rDeltaT, rho, psi, phi, phiCorr, Sp, Su);
 }
 
 
@@ -104,6 +116,7 @@ void Foam::MULES::LTScorrect
 (
     const RhoType& rho,
     volScalarField& psi,
+    const surfaceScalarField& phi,
     surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su,
@@ -123,8 +136,18 @@ void Foam::MULES::LTScorrect
         readLabel(MULEScontrols.lookup("nLimiterIter"))
     );
 
-    limitCorr(rDeltaT, rho, psi, phiCorr, Sp, Su, psiMax, psiMin, nLimiterIter);
-    correct(rDeltaT, rho, psi, phiCorr, Sp, Su);
+    limitCorr
+    (
+        rDeltaT,
+        rho,
+        psi,
+        phi,
+        phiCorr,
+        Sp, Su,
+        psiMax, psiMin,
+        nLimiterIter
+    );
+    correct(rDeltaT, rho, psi, phi, phiCorr, Sp, Su);
 }
 
 
@@ -135,6 +158,7 @@ void Foam::MULES::limiterCorr
     const RdeltaTType& rDeltaT,
     const RhoType& rho,
     const volScalarField& psi,
+    const surfaceScalarField& phi,
     const surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su,
@@ -153,6 +177,9 @@ void Foam::MULES::limiterCorr
     tmp<volScalarField::DimensionedInternalField> tVsc = mesh.Vsc();
     const scalarField& V = tVsc();
 
+    const surfaceScalarField::GeometricBoundaryField& phiBf =
+        phi.boundaryField();
+
     const scalarField& phiCorrIf = phiCorr;
     const surfaceScalarField::GeometricBoundaryField& phiCorrBf =
         phiCorr.boundaryField();
@@ -408,12 +435,12 @@ void Foam::MULES::limiterCorr
             {
                 const labelList& pFaceCells =
                     mesh.boundary()[patchi].faceCells();
-                const scalarField& phiCorrPf = phiCorrBf[patchi];
+                const scalarField& phiPf = phiBf[patchi];
 
                 forAll(lambdaPf, pFacei)
                 {
                     // Limit outlet faces only
-                    if (phiCorrPf[pFacei] > SMALL*SMALL)
+                    if (phiPf[pFacei] > SMALL*SMALL)
                     {
                         label pfCelli = pFaceCells[pFacei];
 
@@ -443,6 +470,7 @@ void Foam::MULES::limitCorr
     const RdeltaTType& rDeltaT,
     const RhoType& rho,
     const volScalarField& psi,
+    const surfaceScalarField& phi,
     surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su,
@@ -478,6 +506,7 @@ void Foam::MULES::limitCorr
         rDeltaT,
         rho,
         psi,
+        phi,
         phiCorr,
         Sp,
         Su,