diff --git a/applications/solvers/multiphase/VoF/alphaEqn.H b/applications/solvers/multiphase/VoF/alphaEqn.H
index ad9ddcb5e5..a9a959fe45 100644
--- a/applications/solvers/multiphase/VoF/alphaEqn.H
+++ b/applications/solvers/multiphase/VoF/alphaEqn.H
@@ -133,7 +133,15 @@
         if (alphaApplyPrevCorr && talphaPhi1Corr0.valid())
         {
             Info<< "Applying the previous iteration compression flux" << endl;
-            MULES::correct(alpha1, alphaPhi10, talphaPhi1Corr0.ref(), 1, 0);
+            MULES::correct
+            (
+                geometricOneField(),
+                alpha1,
+                alphaPhi10,
+                talphaPhi1Corr0.ref(),
+                oneField(),
+                zeroField()
+            );
 
             alphaPhi10 += talphaPhi1Corr0();
         }
@@ -182,8 +190,8 @@
                 talphaPhi1Corr.ref(),
                 Sp,
                 (-Sp*alpha1)(),
-                1,
-                0
+                oneField(),
+                zeroField()
             );
 
             // Under-relax the correction for all but the 1st corrector
@@ -209,8 +217,8 @@
                 alphaPhi10,
                 Sp,
                 (Su + divU*min(alpha1(), scalar(1)))(),
-                1,
-                0
+                oneField(),
+                zeroField()
             );
         }
 
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C
index 2263096cff..1eddd3e398 100644
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C
@@ -1051,8 +1051,8 @@ void Foam::multiphaseMixtureThermo::solveAlphas
             alphaPhiCorr,
             zeroField(),
             zeroField(),
-            1,
-            0,
+            oneField(),
+            zeroField(),
             true
         );
 
diff --git a/applications/solvers/multiphase/driftFluxFoam/alphaEqn.H b/applications/solvers/multiphase/driftFluxFoam/alphaEqn.H
index 3938fde30f..0a64efbcd8 100644
--- a/applications/solvers/multiphase/driftFluxFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/driftFluxFoam/alphaEqn.H
@@ -32,11 +32,12 @@
 
             MULES::correct
             (
+                geometricOneField(),
                 alpha1,
                 alphaPhi,
                 talphaPhiCorr0.ref(),
-                mixture.alphaMax(),
-                0
+                UniformField<scalar>(mixture.alphaMax()),
+                zeroField()
             );
 
             alphaPhi += talphaPhiCorr0();
@@ -71,11 +72,12 @@
 
             MULES::correct
             (
+                geometricOneField(),
                 alpha1,
                 talphaPhiUn(),
                 talphaPhiCorr.ref(),
-                mixture.alphaMax(),
-                0
+                UniformField<scalar>(mixture.alphaMax()),
+                zeroField()
             );
 
             // Under-relax the correction for all but the 1st corrector
@@ -95,11 +97,12 @@
 
             MULES::explicitSolve
             (
+                geometricOneField(),
                 alpha1,
                 phi,
                 alphaPhi,
-                mixture.alphaMax(),
-                0
+                UniformField<scalar>(mixture.alphaMax()),
+                zeroField()
             );
         }
     }
diff --git a/applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqn.H b/applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqn.H
index bbcd6affd8..e4e0e32f20 100644
--- a/applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqn.H
@@ -87,8 +87,8 @@
                 alphaPhi1,
                 zeroField(),
                 zeroField(),
-                1,
-                0
+                oneField(),
+                zeroField()
             );
         }
         else
@@ -103,8 +103,8 @@
                 alphaPhi1,
                 zeroField(),
                 zeroField(),
-                1,
-                0
+                oneField(),
+                zeroField()
             );
         }
 
@@ -150,8 +150,8 @@
                 alphaPhi2,
                 zeroField(),
                 zeroField(),
-                1,
-                0
+                oneField(),
+                zeroField()
             );
         }
         else
@@ -166,8 +166,8 @@
                 alphaPhi2,
                 zeroField(),
                 zeroField(),
-                1,
-                0
+                oneField(),
+                zeroField()
             );
         }
 
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H
index 150ccd1f78..e3bbba19d8 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H
@@ -78,8 +78,8 @@
                     divU*(alpha10 - alpha100)
                   - vDotvmcAlphal*alpha10
                 )(),
-                1,
-                0
+                oneField(),
+                zeroField()
             );
 
             // Under-relax the correction for all but the 1st corrector
@@ -103,8 +103,8 @@
                 talphaPhiCorr.ref(),
                 vDotvmcAlphal,
                 (divU*alpha1 + vDotcAlphal)(),
-                1,
-                0
+                oneField(),
+                zeroField()
             );
 
             talphaPhi = talphaPhiCorr;
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
index 7b82d3d1f5..5cc9422af5 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
@@ -132,8 +132,8 @@ void Foam::multiphaseSystem::solveAlphas()
             alphaPhiCorr,
             zeroField(),
             zeroField(),
-            1,
-            0,
+            oneField(),
+            zeroField(),
             true
         );
 
@@ -168,9 +168,7 @@ void Foam::multiphaseSystem::solveAlphas()
         (
             geometricOneField(),
             phase,
-            alphaPhi,
-            zeroField(),
-            zeroField()
+            alphaPhi
         );
 
         phase.alphaPhi() = alphaPhi;
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
index 20c8c6e8cf..119e9764c6 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
+++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
@@ -621,8 +621,8 @@ void Foam::multiphaseMixture::solveAlphas
             alphaPhiCorr,
             zeroField(),
             zeroField(),
-            1,
-            0,
+            oneField(),
+            zeroField(),
             true
         );
 
@@ -658,9 +658,7 @@ void Foam::multiphaseMixture::solveAlphas
         (
             geometricOneField(),
             alpha,
-            alphaPhi,
-            zeroField(),
-            zeroField()
+            alphaPhi
         );
 
         rhoPhi_ += alphaPhi*alpha.rho();
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
index 93d7ec3fb5..79aa52e03c 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
@@ -69,8 +69,6 @@ void Foam::multiphaseSystem::calcAlphas()
 
 void Foam::multiphaseSystem::solveAlphas()
 {
-    bool LTS = fv::localEulerDdt::enabled(mesh_);
-
     forAll(phases(), phasei)
     {
         phases()[phasei].correctBoundaryConditions();
@@ -138,40 +136,18 @@ void Foam::multiphaseSystem::solveAlphas()
 
         phase.correctInflowOutflow(alphaPhiCorr);
 
-        if (LTS)
-        {
-            MULES::limit
-            (
-                fv::localEulerDdt::localRDeltaT(mesh_),
-                geometricOneField(),
-                phase,
-                phi_,
-                alphaPhiCorr,
-                zeroField(),
-                zeroField(),
-                phase.alphaMax(),
-                0,
-                true
-            );
-        }
-        else
-        {
-            const scalar rDeltaT = 1.0/mesh_.time().deltaTValue();
-
-            MULES::limit
-            (
-                rDeltaT,
-                geometricOneField(),
-                phase,
-                phi_,
-                alphaPhiCorr,
-                zeroField(),
-                zeroField(),
-                phase.alphaMax(),
-                0,
-                true
-            );
-        }
+        MULES::limit
+        (
+            geometricOneField(),
+            phase,
+            phi_,
+            alphaPhiCorr,
+            zeroField(),
+            zeroField(),
+            UniformField<scalar>(phase.alphaMax()),
+            zeroField(),
+            true
+        );
     }
 
     MULES::limitSum(alphaPhiCorrs);
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
index 05e767253f..796df17ccf 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
@@ -274,8 +274,8 @@ void Foam::twoPhaseSystem::solve()
                     alphaPhi10,
                     (alphaSubCycle.index()*Sp)(),
                     (Su - (alphaSubCycle.index() - 1)*Sp*alpha1)(),
-                    phase1_.alphaMax(),
-                    0
+                    UniformField<scalar>(phase1_.alphaMax()),
+                    zeroField()
                 );
 
                 if (alphaSubCycle.index() == 1)
@@ -300,8 +300,8 @@ void Foam::twoPhaseSystem::solve()
                 alphaPhi1,
                 Sp,
                 Su,
-                phase1_.alphaMax(),
-                0
+                UniformField<scalar>(phase1_.alphaMax()),
+                zeroField()
             );
 
             phase1_.alphaPhi() = alphaPhi1;
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H
index 5344fb75f2..6e9888c1c5 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H
@@ -12,7 +12,15 @@
         )
     );
 
-    MULES::explicitSolve(alpha1, phi, alphaPhi, 1, 0);
+    MULES::explicitSolve
+    (
+        geometricOneField(),
+        alpha1,
+        phi,
+        alphaPhi,
+        oneField(),
+        zeroField()
+    );
 
     rhoPhi = alphaPhi*(rho1 - rho2) + phi*rho2;
 
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
index b5ad975720..5e9761b4fe 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
@@ -45,6 +45,7 @@ License
 #include "fixedValueFvsPatchFields.H"
 #include "blendingMethod.H"
 #include "HashPtrTable.H"
+#include "UniformField.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -464,8 +465,8 @@ void Foam::twoPhaseSystem::solve()
                     alphaPhic10,
                     (alphaSubCycle.index()*Sp)(),
                     (Su - (alphaSubCycle.index() - 1)*Sp*alpha1)(),
-                    phase1_.alphaMax(),
-                    0
+                    UniformField<scalar>(phase1_.alphaMax()),
+                    zeroField()
                 );
 
                 if (alphaSubCycle.index() == 1)
@@ -490,8 +491,8 @@ void Foam::twoPhaseSystem::solve()
                 alphaPhic1,
                 Sp,
                 Su,
-                phase1_.alphaMax(),
-                0
+                UniformField<scalar>(phase1_.alphaMax()),
+                zeroField()
             );
 
             phase1_.alphaPhi() = alphaPhic1;
diff --git a/src/OpenFOAM/fields/Fields/uniformField/UniformField.H b/src/OpenFOAM/fields/Fields/uniformField/UniformField.H
index 140b4b4fd6..3b9f958a2f 100644
--- a/src/OpenFOAM/fields/Fields/uniformField/UniformField.H
+++ b/src/OpenFOAM/fields/Fields/uniformField/UniformField.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -61,6 +61,8 @@ public:
 
     // Member Operators
 
+        inline operator Type() const;
+
         inline Type operator[](const label) const;
 
         inline UniformField field() const;
diff --git a/src/OpenFOAM/fields/Fields/uniformField/UniformFieldI.H b/src/OpenFOAM/fields/Fields/uniformField/UniformFieldI.H
index 9c7012a7ee..9bcf35579a 100644
--- a/src/OpenFOAM/fields/Fields/uniformField/UniformFieldI.H
+++ b/src/OpenFOAM/fields/Fields/uniformField/UniformFieldI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -34,6 +34,13 @@ inline Foam::UniformField<Type>::UniformField(const Type& value)
 {}
 
 
+template<class Type>
+inline Foam::UniformField<Type>::operator Type() const
+{
+    return value_;
+}
+
+
 template<class Type>
 inline Type Foam::UniformField<Type>::operator[](const label) const
 {
@@ -48,4 +55,65 @@ inline Foam::UniformField<Type> Foam::UniformField<Type>::field() const
 }
 
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+inline UniformField<Type> min
+(
+    const UniformField<Type>& u1,
+    const UniformField<Type>& u2
+)
+{
+    return UniformField<Type>(min(u1.operator Type(), u2.operator Type()));
+}
+
+
+template<class Type, class OtherType>
+inline OtherType min(const UniformField<Type>& u, const OtherType& o)
+{
+    return min(u.operator Type(), o);
+}
+
+
+template<class Type, class OtherType>
+inline OtherType min(const OtherType& o, const UniformField<Type>& u)
+{
+    return min(o, u.operator Type());
+}
+
+
+template<class Type>
+inline UniformField<Type> max
+(
+    const UniformField<Type>& u1,
+    const UniformField<Type>& u2
+)
+{
+    return UniformField<Type>(max(u1.operator Type(), u2.operator Type()));
+}
+
+
+template<class Type, class OtherType>
+inline OtherType max(const UniformField<Type>& u, const OtherType& o)
+{
+    return max(u.operator Type(), o);
+}
+
+
+template<class Type, class OtherType>
+inline OtherType max(const OtherType& o, const UniformField<Type>& u)
+{
+    return max(o, u.operator Type());
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
 // ************************************************************************* //
diff --git a/src/OpenFOAM/primitives/one/oneI.H b/src/OpenFOAM/primitives/one/oneI.H
index 05be7ad120..b85f7caba3 100644
--- a/src/OpenFOAM/primitives/one/oneI.H
+++ b/src/OpenFOAM/primitives/one/oneI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -81,6 +81,39 @@ inline const Type& operator/(const Type& t, const one&)
     return t;
 }
 
+inline const one& min(const one& o, const one&)
+{
+    return o;
+}
+
+template<class Type>
+inline Type min(const one&, const Type& t)
+{
+    return min(scalar(1), t);
+}
+
+template<class Type>
+inline Type min(const Type& t, const one&)
+{
+    return min(t, scalar(1));
+}
+
+inline const one& max(const one& o, const one&)
+{
+    return o;
+}
+
+template<class Type>
+inline Type max(const one&, const Type& t)
+{
+    return max(scalar(1), t);
+}
+
+template<class Type>
+inline Type max(const Type& t, const one&)
+{
+    return max(t, scalar(1));
+}
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/OpenFOAM/primitives/zero/zeroI.H b/src/OpenFOAM/primitives/zero/zeroI.H
index d69255f159..b4bcb0c2be 100644
--- a/src/OpenFOAM/primitives/zero/zeroI.H
+++ b/src/OpenFOAM/primitives/zero/zeroI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "zero.H"
+#include "scalar.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -89,6 +90,39 @@ inline zero operator/(const zero&, const Type& t)
     return Zero;
 }
 
+inline zero min(const zero&, const zero&)
+{
+    return Zero;
+}
+
+template<class Type>
+inline Type min(const zero&, const Type& t)
+{
+    return min(scalar(0), t);
+}
+
+template<class Type>
+inline Type min(const Type& t, const zero&)
+{
+    return min(t, scalar(0));
+}
+
+inline zero max(const zero&, const zero&)
+{
+    return Zero;
+}
+
+template<class Type>
+inline Type max(const zero&, const Type& t)
+{
+    return max(scalar(0), t);
+}
+
+template<class Type>
+inline Type max(const Type& t, const zero&)
+{
+    return max(t, scalar(0));
+}
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files
index da45466621..5f56d4e7f4 100644
--- a/src/finiteVolume/Make/files
+++ b/src/finiteVolume/Make/files
@@ -240,7 +240,6 @@ fields/surfaceFields/surfaceFields.C
 fvMatrices/fvMatrices.C
 fvMatrices/fvScalarMatrix/fvScalarMatrix.C
 fvMatrices/solvers/MULES/MULES.C
-fvMatrices/solvers/MULES/CMULES.C
 fvMatrices/solvers/GAMGSymSolver/GAMGAgglomerations/faceAreaPairGAMGAgglomeration/faceAreaPairGAMGAgglomeration.C
 
 interpolation = interpolation/interpolation
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/CMULES.C b/src/finiteVolume/fvMatrices/solvers/MULES/CMULES.C
deleted file mode 100644
index d5cd3fba5d..0000000000
--- a/src/finiteVolume/fvMatrices/solvers/MULES/CMULES.C
+++ /dev/null
@@ -1,51 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "CMULES.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-void Foam::MULES::correct
-(
-    volScalarField& psi,
-    const surfaceScalarField& phi,
-    surfaceScalarField& phiPsiCorr,
-    const scalar psiMax,
-    const scalar psiMin
-)
-{
-    correct
-    (
-        geometricOneField(),
-        psi,
-        phi,
-        phiPsiCorr,
-        zeroField(), zeroField(),
-        psiMax, psiMin
-    );
-}
-
-
-// ************************************************************************* //
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/CMULES.H b/src/finiteVolume/fvMatrices/solvers/MULES/CMULES.H
index fd76c9f2fd..cffdca571f 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/CMULES.H
+++ b/src/finiteVolume/fvMatrices/solvers/MULES/CMULES.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -66,14 +66,49 @@ void correct
     const RdeltaTType& rDeltaT,
     const RhoType& rho,
     volScalarField& psi,
-    const surfaceScalarField& phi,
     const surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su
 );
 
+template<class RhoType>
+void correct
+(
+    const RhoType& rho,
+    volScalarField& psi,
+    const surfaceScalarField& phiCorr
+);
+
 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 PsiMaxType, class PsiMinType>
+void correct
+(
+    const RhoType& rho,
+    volScalarField& psi,
+    const surfaceScalarField& phi,
+    surfaceScalarField& phiCorr,
+    const PsiMaxType& psiMax,
+    const PsiMinType& psiMin
+);
+
+template
+<
+    class RhoType,
+    class SpType,
+    class SuType,
+    class PsiMaxType,
+    class PsiMinType
+>
+void correct
 (
     const RhoType& rho,
     volScalarField& psi,
@@ -81,20 +116,19 @@ void correct
     surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su,
-    const scalar psiMax,
-    const scalar psiMin
+    const PsiMaxType& psiMax,
+    const PsiMinType& psiMin
 );
 
-void correct
-(
-    volScalarField& psi,
-    const surfaceScalarField& phi,
-    surfaceScalarField& phiCorr,
-    const scalar psiMax,
-    const scalar psiMin
-);
-
-template<class RdeltaTType, class RhoType, class SpType, class SuType>
+template
+<
+    class RdeltaTType,
+    class RhoType,
+    class SpType,
+    class SuType,
+    class PsiMaxType,
+    class PsiMinType
+>
 void limiterCorr
 (
     scalarField& allLambda,
@@ -105,11 +139,19 @@ void limiterCorr
     const surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su,
-    const scalar psiMax,
-    const scalar psiMin
+    const PsiMaxType& psiMax,
+    const PsiMinType& psiMin
 );
 
-template<class RdeltaTType, class RhoType, class SpType, class SuType>
+template
+<
+    class RdeltaTType,
+    class RhoType,
+    class SpType,
+    class SuType,
+    class PsiMaxType,
+    class PsiMinType
+>
 void limitCorr
 (
     const RdeltaTType& rDeltaT,
@@ -119,8 +161,8 @@ void limitCorr
     surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su,
-    const scalar psiMax,
-    const scalar psiMin
+    const PsiMaxType& psiMax,
+    const PsiMinType& psiMin
 );
 
 
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C b/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C
index beb0cf66ba..be51f1a9c5 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C
+++ b/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C
@@ -38,7 +38,6 @@ void Foam::MULES::correct
     const RdeltaTType& rDeltaT,
     const RhoType& rho,
     volScalarField& psi,
-    const surfaceScalarField& phi,
     const surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su
@@ -74,8 +73,67 @@ void Foam::MULES::correct
 }
 
 
+template<class RhoType>
+void Foam::MULES::correct
+(
+    const RhoType& rho,
+    volScalarField& psi,
+    const surfaceScalarField& phiCorr
+)
+{
+    correct(rho, psi, phiCorr, zeroField(), zeroField());
+}
+
+
 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
+)
+{
+    const fvMesh& mesh = psi.mesh();
+
+    if (fv::localEulerDdt::enabled(mesh))
+    {
+        const volScalarField& rDeltaT = fv::localEulerDdt::localRDeltaT(mesh);
+        correct(rDeltaT, rho, psi, phiCorr, Sp, Su);
+    }
+    else
+    {
+        const scalar rDeltaT = 1.0/mesh.time().deltaTValue();
+        correct(rDeltaT, rho, psi, phiCorr, Sp, Su);
+    }
+}
+
+
+template<class RhoType, class PsiMaxType, class PsiMinType>
+void Foam::MULES::correct
+(
+    const RhoType& rho,
+    volScalarField& psi,
+    const surfaceScalarField& phi,
+    surfaceScalarField& phiCorr,
+    const PsiMaxType& psiMax,
+    const PsiMinType& psiMin
+)
+{
+    correct(rho, psi, phi, phiCorr, zeroField(), zeroField(), psiMax, psiMin);
+}
+
+
+template
+<
+    class RhoType,
+    class SpType,
+    class SuType,
+    class PsiMaxType,
+    class PsiMinType
+>
+void Foam::MULES::correct
 (
     const RhoType& rho,
     volScalarField& psi,
@@ -83,8 +141,8 @@ void Foam::MULES::correct
     surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su,
-    const scalar psiMax,
-    const scalar psiMin
+    const PsiMaxType& psiMax,
+    const PsiMinType& psiMin
 )
 {
     const fvMesh& mesh = psi.mesh();
@@ -105,7 +163,8 @@ void Foam::MULES::correct
             psiMax,
             psiMin
         );
-        correct(rDeltaT, rho, psi, phi, phiCorr, Sp, Su);
+
+        correct(rDeltaT, rho, psi, phiCorr, Sp, Su);
     }
     else
     {
@@ -124,12 +183,20 @@ void Foam::MULES::correct
             psiMin
         );
 
-        correct(rDeltaT, rho, psi, phi, phiCorr, Sp, Su);
+        correct(rDeltaT, rho, psi, phiCorr, Sp, Su);
     }
 }
 
 
-template<class RdeltaTType, class RhoType, class SpType, class SuType>
+template
+<
+    class RdeltaTType,
+    class RhoType,
+    class SpType,
+    class SuType,
+    class PsiMaxType,
+    class PsiMinType
+>
 void Foam::MULES::limiterCorr
 (
     scalarField& allLambda,
@@ -140,8 +207,8 @@ void Foam::MULES::limiterCorr
     const surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su,
-    const scalar psiMax,
-    const scalar psiMin
+    const PsiMaxType& psiMax,
+    const PsiMinType& psiMin
 )
 {
     const scalarField& psiIf = psi;
@@ -199,8 +266,11 @@ void Foam::MULES::limiterCorr
     surfaceScalarField::Boundary& lambdaBf =
         lambda.boundaryFieldRef();
 
-    scalarField psiMaxn(psiIf.size(), psiMin);
-    scalarField psiMinn(psiIf.size(), psiMax);
+    scalarField psiMaxn(psiIf.size());
+    scalarField psiMinn(psiIf.size());
+
+    psiMaxn = psiMin;
+    psiMinn = psiMax;
 
     scalarField sumPhip(psiIf.size(), 0.0);
     scalarField mSumPhim(psiIf.size(), 0.0);
@@ -265,8 +335,8 @@ void Foam::MULES::limiterCorr
             {
                 const label pfCelli = pFaceCells[pFacei];
 
-                psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiMax);
-                psiMinn[pfCelli] = min(psiMinn[pfCelli], psiMin);
+                psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiMax[pfCelli]);
+                psiMinn[pfCelli] = min(psiMinn[pfCelli], psiMin[pfCelli]);
             }
         }
 
@@ -472,7 +542,15 @@ void Foam::MULES::limiterCorr
 }
 
 
-template<class RdeltaTType, class RhoType, class SpType, class SuType>
+template
+<
+    class RdeltaTType,
+    class RhoType,
+    class SpType,
+    class SuType,
+    class PsiMaxType,
+    class PsiMinType
+>
 void Foam::MULES::limitCorr
 (
     const RdeltaTType& rDeltaT,
@@ -482,8 +560,8 @@ void Foam::MULES::limitCorr
     surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su,
-    const scalar psiMax,
-    const scalar psiMin
+    const PsiMaxType& psiMax,
+    const PsiMinType& psiMin
 )
 {
     const fvMesh& mesh = psi.mesh();
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULES.C b/src/finiteVolume/fvMatrices/solvers/MULES/MULES.C
index 76ba38f969..77ab25c3b1 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/MULES.C
+++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULES.C
@@ -27,27 +27,6 @@ License
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-void Foam::MULES::explicitSolve
-(
-    volScalarField& psi,
-    const surfaceScalarField& phi,
-    surfaceScalarField& phiPsi,
-    const scalar psiMax,
-    const scalar psiMin
-)
-{
-    explicitSolve
-    (
-        geometricOneField(),
-        psi,
-        phi,
-        phiPsi,
-        zeroField(), zeroField(),
-        psiMax, psiMin
-    );
-}
-
-
 void Foam::MULES::limitSum(UPtrList<scalarField>& phiPsiCorrs)
 {
     forAll(phiPsiCorrs[0], facei)
@@ -97,4 +76,41 @@ void Foam::MULES::limitSum(UPtrList<scalarField>& phiPsiCorrs)
 }
 
 
+void Foam::MULES::limitSum
+(
+    const UPtrList<const scalarField>& alphas,
+    UPtrList<scalarField>& phiPsiCorrs,
+    const labelHashSet& fixed
+)
+{
+    labelHashSet notFixed(identity(phiPsiCorrs.size()));
+    notFixed -= fixed;
+
+    forAll(phiPsiCorrs[0], facei)
+    {
+        scalar alphaNotFixed = 0, corrNotFixed = 0;
+        forAllConstIter(labelHashSet, notFixed, iter)
+        {
+            alphaNotFixed += alphas[iter.key()][facei];
+            corrNotFixed += phiPsiCorrs[iter.key()][facei];
+        }
+
+        scalar corrFixed = 0;
+        forAllConstIter(labelHashSet, fixed, iter)
+        {
+            corrFixed += phiPsiCorrs[iter.key()][facei];
+        }
+
+        const scalar sumCorr = corrNotFixed + corrFixed;
+
+        const scalar lambda = - sumCorr/alphaNotFixed;
+
+        forAllConstIter(labelHashSet, notFixed, iter)
+        {
+            phiPsiCorrs[iter.key()][facei] += lambda*alphas[iter.key()][facei];
+        }
+    }
+}
+
+
 // ************************************************************************* //
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULES.H b/src/finiteVolume/fvMatrices/solvers/MULES/MULES.H
index 5ece4fefcb..749fbd2c3c 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/MULES.H
+++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULES.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -50,6 +50,8 @@ SourceFiles
 #include "zero.H"
 #include "zeroField.H"
 #include "UPtrList.H"
+#include "HashSet.H"
+#include "UniformField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -71,6 +73,14 @@ void explicitSolve
     const SuType& Su
 );
 
+template<class RhoType>
+void explicitSolve
+(
+    const RhoType& rho,
+    volScalarField& psi,
+    const surfaceScalarField& phiPsi
+);
+
 template<class RhoType, class SpType, class SuType>
 void explicitSolve
 (
@@ -81,7 +91,25 @@ void explicitSolve
     const SuType& Su
 );
 
-template<class RhoType, class SpType, class SuType>
+template<class RhoType, class PsiMaxType, class PsiMinType>
+void explicitSolve
+(
+    const RhoType& rho,
+    volScalarField& psi,
+    const surfaceScalarField& phiBD,
+    surfaceScalarField& phiPsi,
+    const PsiMaxType& psiMax,
+    const PsiMinType& psiMin
+);
+
+template
+<
+    class RhoType,
+    class SpType,
+    class SuType,
+    class PsiMaxType,
+    class PsiMinType
+>
 void explicitSolve
 (
     const RhoType& rho,
@@ -90,20 +118,19 @@ void explicitSolve
     surfaceScalarField& phiPsi,
     const SpType& Sp,
     const SuType& Su,
-    const scalar psiMax,
-    const scalar psiMin
+    const PsiMaxType& psiMax,
+    const PsiMinType& psiMin
 );
 
-void explicitSolve
-(
-    volScalarField& psi,
-    const surfaceScalarField& phiBD,
-    surfaceScalarField& phiPsi,
-    const scalar psiMax,
-    const scalar psiMin
-);
-
-template<class RdeltaTType, class RhoType, class SpType, class SuType>
+template
+<
+    class RdeltaTType,
+    class RhoType,
+    class SpType,
+    class SuType,
+    class PsiMaxType,
+    class PsiMinType
+>
 void limiter
 (
     scalarField& allLambda,
@@ -114,11 +141,19 @@ void limiter
     const surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su,
-    const scalar psiMax,
-    const scalar psiMin
+    const PsiMaxType& psiMax,
+    const PsiMinType& psiMin
 );
 
-template<class RdeltaTType, class RhoType, class SpType, class SuType>
+template
+<
+    class RdeltaTType,
+    class RhoType,
+    class SpType,
+    class SuType,
+    class PsiMaxType,
+    class PsiMinType
+>
 void limit
 (
     const RdeltaTType& rDeltaT,
@@ -128,17 +163,51 @@ void limit
     surfaceScalarField& phiPsi,
     const SpType& Sp,
     const SuType& Su,
-    const scalar psiMax,
-    const scalar psiMin,
+    const PsiMaxType& psiMax,
+    const PsiMinType& psiMin,
     const bool returnCorr
 );
 
+template
+<
+    class RhoType,
+    class SpType,
+    class SuType,
+    class PsiMaxType,
+    class PsiMinType
+>
+void limit
+(
+    const RhoType& rho,
+    const volScalarField& psi,
+    const surfaceScalarField& phi,
+    surfaceScalarField& phiPsi,
+    const SpType& Sp,
+    const SuType& Su,
+    const PsiMaxType& psiMax,
+    const PsiMinType& psiMin,
+    const bool returnCorr
+);
 
 void limitSum(UPtrList<scalarField>& phiPsiCorrs);
 
 template<class SurfaceScalarFieldList>
 void limitSum(SurfaceScalarFieldList& phiPsiCorrs);
 
+void limitSum
+(
+    const UPtrList<const scalarField>& alphas,
+    UPtrList<scalarField>& phiPsiCorrs,
+    const labelHashSet& fixed
+);
+
+template<class SurfaceScalarFieldList>
+void limitSum
+(
+    const SurfaceScalarFieldList& alphas,
+    SurfaceScalarFieldList& phiPsiCorrs,
+    const labelHashSet& fixed
+);
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
index a7f0c17055..ac4aa598c3 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
+++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
@@ -78,6 +78,18 @@ void Foam::MULES::explicitSolve
 }
 
 
+template<class RhoType>
+void Foam::MULES::explicitSolve
+(
+    const RhoType& rho,
+    volScalarField& psi,
+    const surfaceScalarField& phiPsi
+)
+{
+    explicitSolve(rho, psi, phiPsi, zeroField(), zeroField());
+}
+
+
 template<class RhoType, class SpType, class SuType>
 void Foam::MULES::explicitSolve
 (
@@ -103,7 +115,39 @@ void Foam::MULES::explicitSolve
 }
 
 
-template<class RhoType, class SpType, class SuType>
+template<class RhoType, class PsiMaxType, class PsiMinType>
+void Foam::MULES::explicitSolve
+(
+    const RhoType& rho,
+    volScalarField& psi,
+    const surfaceScalarField& phiBD,
+    surfaceScalarField& phiPsi,
+    const PsiMaxType& psiMax,
+    const PsiMinType& psiMin
+)
+{
+    explicitSolve
+    (
+        rho,
+        psi,
+        phiBD,
+        phiPsi,
+        zeroField(),
+        zeroField(),
+        psiMax,
+        psiMin
+    );
+}
+
+
+template
+<
+    class RhoType,
+    class SpType,
+    class SuType,
+    class PsiMaxType,
+    class PsiMinType
+>
 void Foam::MULES::explicitSolve
 (
     const RhoType& rho,
@@ -112,8 +156,8 @@ void Foam::MULES::explicitSolve
     surfaceScalarField& phiPsi,
     const SpType& Sp,
     const SuType& Su,
-    const scalar psiMax,
-    const scalar psiMin
+    const PsiMaxType& psiMax,
+    const PsiMinType& psiMin
 )
 {
     const fvMesh& mesh = psi.mesh();
@@ -123,47 +167,27 @@ void Foam::MULES::explicitSolve
     if (fv::localEulerDdt::enabled(mesh))
     {
         const volScalarField& rDeltaT = fv::localEulerDdt::localRDeltaT(mesh);
-
-        limit
-        (
-            rDeltaT,
-            rho,
-            psi,
-            phi,
-            phiPsi,
-            Sp,
-            Su,
-            psiMax,
-            psiMin,
-            false
-        );
-
+        limit(rDeltaT, rho, psi, phi, phiPsi, Sp, Su, psiMax, psiMin, false);
         explicitSolve(rDeltaT, rho, psi, phiPsi, Sp, Su);
     }
     else
     {
         const scalar rDeltaT = 1.0/mesh.time().deltaTValue();
-
-        limit
-        (
-            rDeltaT,
-            rho,
-            psi,
-            phi,
-            phiPsi,
-            Sp,
-            Su,
-            psiMax,
-            psiMin,
-            false
-        );
-
+        limit(rDeltaT, rho, psi, phi, phiPsi, Sp, Su, psiMax, psiMin, false);
         explicitSolve(rDeltaT, rho, psi, phiPsi, Sp, Su);
     }
 }
 
 
-template<class RdeltaTType, class RhoType, class SpType, class SuType>
+template
+<
+    class RdeltaTType,
+    class RhoType,
+    class SpType,
+    class SuType,
+    class PsiMaxType,
+    class PsiMinType
+>
 void Foam::MULES::limiter
 (
     scalarField& allLambda,
@@ -174,8 +198,8 @@ void Foam::MULES::limiter
     const surfaceScalarField& phiCorr,
     const SpType& Sp,
     const SuType& Su,
-    const scalar psiMax,
-    const scalar psiMin
+    const PsiMaxType& psiMax,
+    const PsiMinType& psiMin
 )
 {
     const scalarField& psiIf = psi;
@@ -235,8 +259,11 @@ void Foam::MULES::limiter
     scalarField& lambdaIf = lambda;
     surfaceScalarField::Boundary& lambdaBf = lambda.boundaryFieldRef();
 
-    scalarField psiMaxn(psiIf.size(), psiMin);
-    scalarField psiMinn(psiIf.size(), psiMax);
+    scalarField psiMaxn(psiIf.size());
+    scalarField psiMinn(psiIf.size());
+
+    psiMaxn = psiMin;
+    psiMinn = psiMax;
 
     scalarField sumPhiBD(psiIf.size(), 0.0);
 
@@ -307,8 +334,8 @@ void Foam::MULES::limiter
             {
                 const label pfCelli = pFaceCells[pFacei];
 
-                psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiMax);
-                psiMinn[pfCelli] = min(psiMinn[pfCelli], psiMin);
+                psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiMax[pfCelli]);
+                psiMinn[pfCelli] = min(psiMinn[pfCelli], psiMin[pfCelli]);
             }
         }
 
@@ -516,7 +543,15 @@ void Foam::MULES::limiter
 }
 
 
-template<class RdeltaTType, class RhoType, class SpType, class SuType>
+template
+<
+    class RdeltaTType,
+    class RhoType,
+    class SpType,
+    class SuType,
+    class PsiMaxType,
+    class PsiMinType
+>
 void Foam::MULES::limit
 (
     const RdeltaTType& rDeltaT,
@@ -526,8 +561,8 @@ void Foam::MULES::limit
     surfaceScalarField& phiPsi,
     const SpType& Sp,
     const SuType& Su,
-    const scalar psiMax,
-    const scalar psiMin,
+    const PsiMaxType& psiMax,
+    const PsiMinType& psiMin,
     const bool returnCorr
 )
 {
@@ -595,6 +630,42 @@ void Foam::MULES::limit
 }
 
 
+template
+<
+    class RhoType,
+    class SpType,
+    class SuType,
+    class PsiMaxType,
+    class PsiMinType
+>
+void Foam::MULES::limit
+(
+    const RhoType& rho,
+    const volScalarField& psi,
+    const surfaceScalarField& phi,
+    surfaceScalarField& phiPsi,
+    const SpType& Sp,
+    const SuType& Su,
+    const PsiMaxType& psiMax,
+    const PsiMinType& psiMin,
+    const bool rtnCorr
+)
+{
+    const fvMesh& mesh = psi.mesh();
+
+    if (fv::localEulerDdt::enabled(mesh))
+    {
+        const volScalarField& rDeltaT = fv::localEulerDdt::localRDeltaT(mesh);
+        limit(rDeltaT, rho, psi, phi, phiPsi, Sp, Su, psiMax, psiMin, rtnCorr);
+    }
+    else
+    {
+        const scalar rDeltaT = 1.0/mesh.time().deltaTValue();
+        limit(rDeltaT, rho, psi, phi, phiPsi, Sp, Su, psiMax, psiMin, rtnCorr);
+    }
+}
+
+
 template<class SurfaceScalarFieldList>
 void Foam::MULES::limitSum(SurfaceScalarFieldList& phiPsiCorrs)
 {
@@ -631,4 +702,59 @@ void Foam::MULES::limitSum(SurfaceScalarFieldList& phiPsiCorrs)
 }
 
 
+template<class SurfaceScalarFieldList>
+void Foam::MULES::limitSum
+(
+    const SurfaceScalarFieldList& alphas,
+    SurfaceScalarFieldList& phiPsiCorrs,
+    const labelHashSet& fixed
+)
+{
+    {
+        UPtrList<const scalarField> alphasInternal(alphas.size());
+        forAll(alphas, phasei)
+        {
+            alphasInternal.set(phasei, &alphas[phasei]);
+        }
+        UPtrList<scalarField> phiPsiCorrsInternal(phiPsiCorrs.size());
+        forAll(phiPsiCorrs, phasei)
+        {
+            phiPsiCorrsInternal.set(phasei, &phiPsiCorrs[phasei]);
+        }
+
+        limitSum(alphasInternal, phiPsiCorrsInternal, fixed);
+    }
+
+    const surfaceScalarField::Boundary& bfld =
+        phiPsiCorrs[0].boundaryField();
+
+    forAll(bfld, patchi)
+    {
+        if (bfld[patchi].coupled())
+        {
+            UPtrList<const scalarField> alphasPatch(alphas.size());
+            forAll(alphas, phasei)
+            {
+                alphasPatch.set
+                (
+                    phasei,
+                    &alphas[phasei].boundaryField()[patchi]
+                );
+            }
+            UPtrList<scalarField> phiPsiCorrsPatch(phiPsiCorrs.size());
+            forAll(phiPsiCorrs, phasei)
+            {
+                phiPsiCorrsPatch.set
+                (
+                    phasei,
+                    &phiPsiCorrs[phasei].boundaryFieldRef()[patchi]
+                );
+            }
+
+            limitSum(alphasPatch, phiPsiCorrsPatch, fixed);
+        }
+    }
+}
+
+
 // ************************************************************************* //