MULES: Added optional boundaryExtremaCoeff

This provides more flexibility in specifying the allowed internal and boundary
extrema.

For driftFluxFoam and other settling problems it is beneficial to set the
boundaryExtremaCoeff to 1 to allow rapid accumulation of the partials on the
bottom wall (which was the previous default behaviour) but this is not suitable
for many Euler-Euler cases for which a uniform etrema coefficient is preferable,
either 0 or a small value.

Now by default boundaryExtremaCoeff is set to extremaCoeff which defaults to 0
which provides the behaviour before

OpenFOAM-dev commit cb2bc60fa5

and the driftFluxFoam tutorials have been updated adding

        boundaryExtremaCoeff 1;

to the MULES controls in fvSolution so reproduce the previous behaviour.

src/OpenFOAM/primitives/one/oneI.H

@@ -41,6 +41,30 @@ public:
     typedef arg2 type;
 };
 
+template<class Type>
+inline const Type& operator+(const Type& t, const one&)
+{
+    return t + 1;
+}
+
+template<class Type>
+inline const Type& operator+(const one&, const Type& t)
+{
+    return 1 + t;
+}
+
+template<class Type>
+inline const Type& operator-(const Type& t, const one&)
+{
+    return t - 1;
+}
+
+template<class Type>
+inline const Type& operator-(const one&, const Type& t)
+{
+    return 1 - t;
+}
+
 inline const one& operator*(const one& o, const one&)
 {
     return o;

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

@@ -233,6 +233,20 @@ void Foam::MULES::limiterCorr
         MULEScontrols.lookupOrDefault<scalar>("extremaCoeff", 0)
     );
 
+    const scalar boundaryExtremaCoeff
+    (
+        MULEScontrols.lookupOrDefault<scalar>
+        (
+            "boundaryExtremaCoeff",
+            extremaCoeff
+        )
+    );
+
+    const scalar boundaryDeltaExtremaCoeff
+    (
+        max(boundaryExtremaCoeff - extremaCoeff, 0)
+    );
+
     const labelUList& owner = mesh.owner();
     const labelUList& neighb = mesh.neighbour();
     tmp<volScalarField::Internal> tVsc = mesh.Vsc();
@@ -331,12 +345,20 @@ void Foam::MULES::limiterCorr
         }
         else
         {
-            forAll(phiCorrPf, pFacei)
+            // Add the optional additional allowed boundary extrema
+            if (boundaryDeltaExtremaCoeff > 0)
             {
-                const label pfCelli = pFaceCells[pFacei];
+                forAll(phiCorrPf, pFacei)
+                {
+                    const label pfCelli = pFaceCells[pFacei];
 
-                psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiMax[pfCelli]);
-                psiMinn[pfCelli] = min(psiMinn[pfCelli], psiMin[pfCelli]);
+                    const scalar extrema =
+                        boundaryDeltaExtremaCoeff
+                       *(psiMax[pfCelli] - psiMin[pfCelli]);
+
+                    psiMaxn[pfCelli] += extrema;
+                    psiMinn[pfCelli] -= extrema;
+                }
             }
         }
 

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

@@ -224,6 +224,20 @@ void Foam::MULES::limiter
         MULEScontrols.lookupOrDefault<scalar>("extremaCoeff", 0)
     );
 
+    const scalar boundaryExtremaCoeff
+    (
+        MULEScontrols.lookupOrDefault<scalar>
+        (
+            "boundaryExtremaCoeff",
+            extremaCoeff
+        )
+    );
+
+    const scalar boundaryDeltaExtremaCoeff
+    (
+        max(boundaryExtremaCoeff - extremaCoeff, 0)
+    );
+
     const scalarField& psi0 = psi.oldTime();
 
     const labelUList& owner = mesh.owner();
@@ -330,12 +344,20 @@ void Foam::MULES::limiter
         }
         else
         {
-            forAll(phiCorrPf, pFacei)
+            // Add the optional additional allowed boundary extrema
+            if (boundaryDeltaExtremaCoeff > 0)
             {
-                const label pfCelli = pFaceCells[pFacei];
+                forAll(phiCorrPf, pFacei)
+                {
+                    const label pfCelli = pFaceCells[pFacei];
 
-                psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiMax[pfCelli]);
-                psiMinn[pfCelli] = min(psiMinn[pfCelli], psiMin[pfCelli]);
+                    const scalar extrema =
+                        boundaryDeltaExtremaCoeff
+                       *(psiMax[pfCelli] - psiMin[pfCelli]);
+
+                    psiMaxn[pfCelli] += extrema;
+                    psiMinn[pfCelli] -= extrema;
+                }
             }
         }
 

tutorials/multiphase/driftFluxFoam/RAS/dahl/system/fvSolution

@@ -26,6 +26,8 @@ solvers
         nLimiterIter    3;
         alphaApplyPrevCorr  yes;
 
+        boundaryExtremaCoeff 1;
+
         solver          smoothSolver;
         smoother        symGaussSeidel;
         tolerance       1e-6;

tutorials/multiphase/driftFluxFoam/RAS/mixerVessel2D/system/fvSolution

@@ -26,6 +26,8 @@ solvers
         nLimiterIter    3;
         alphaApplyPrevCorr  yes;
 
+        boundaryExtremaCoeff 1;
+
         solver          smoothSolver;
         smoother        symGaussSeidel;
         tolerance       1e-6;

tutorials/multiphase/driftFluxFoam/RAS/tank3D/system/fvSolution

@@ -26,6 +26,8 @@ solvers
         nLimiterIter    3;
         alphaApplyPrevCorr  yes;
 
+        boundaryExtremaCoeff 1;
+
         solver          smoothSolver;
         smoother        symGaussSeidel;
         tolerance       1e-6;