ENH: use GeometricField clamp_min, clamp_max, clamp_range

- newer naming allows for less confusing code.
  Eg,
      max(lower)  ->  clamp_min(lower)
      min(upper)  ->  clamp_max(upper)

- prefer combined method, for few operations.
  Eg,
      max(lower) + min(upper) -> clamp_range(lower, upper)

  The updated naming also helps avoid some obvious coding errors.
  Eg,

     Re.min(1200.0);
     Re.max(18800.0);

  instead of
     Re.clamp_range(1200.0, 18800.0);

- can also use implicit conversion of zero_one to MinMax<Type> for
  this type of code:

      lambda_.clamp_range(zero_one{});

applications/solvers/combustion/XiFoam/bEqn.H

@@ -171,10 +171,7 @@ if (ign.ignited())
 
         fvOptions.correct(Su);
 
-        // Limit the maximum Su
+        Su.clamp_range(SuMin, SuMax);
-        // ~~~~~~~~~~~~~~~~~~~~
-        Su.min(SuMax);
-        Su.max(SuMin);
     }
     else
     {

applications/solvers/combustion/fireFoam/YEEqn.H

@@ -39,13 +39,13 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 
             fvOptions.correct(Yi);
 
-            Yi.max(0.0);
+            Yi.clamp_min(0);
             Yt += Yi;
         }
     }
 
     Y[inertIndex] = scalar(1) - Yt;
-    Y[inertIndex].max(0.0);
+    Y[inertIndex].clamp_min(0);
 
     radiation->correct();
 

applications/solvers/combustion/reactingFoam/YEqn.H

@@ -38,11 +38,11 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 
             fvOptions.correct(Yi);
 
-            Yi.max(0.0);
+            Yi.clamp_min(0);
             Yt += Yi;
         }
     }
 
     Y[inertIndex] = scalar(1) - Yt;
-    Y[inertIndex].max(0.0);
+    Y[inertIndex].clamp_min(0);
 }

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/YEqn.H

@@ -48,7 +48,7 @@ if (Y.size())
 
             fvOptions.correct(Yi);
 
-            Yi.max(0.0);
+            Yi.clamp_min(0);
             Yt += Yi;
         }
     }
@@ -56,6 +56,6 @@ if (Y.size())
     if (Y.size())
     {
         Y[inertIndex] = scalar(1) - Yt;
-        Y[inertIndex].max(0.0);
+        Y[inertIndex].clamp_min(0);
     }
 }

applications/solvers/lagrangian/coalChemistryFoam/YEqn.H

@@ -40,11 +40,11 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 
             fvOptions.correct(Yi);
 
-            Yi.max(0.0);
+            Yi.clamp_min(0);
             Yt += Yi;
         }
     }
 
     Y[inertIndex] = scalar(1) - Yt;
-    Y[inertIndex].max(0.0);
+    Y[inertIndex].clamp_min(0);
 }

applications/solvers/lagrangian/reactingParcelFoam/YEqn.H

@@ -41,11 +41,11 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 
             fvOptions.correct(Yi);
 
-            Yi.max(0.0);
+            Yi.clamp_min(0);
             Yt += Yi;
         }
     }
 
     Y[inertIndex] = scalar(1) - Yt;
-    Y[inertIndex].max(0.0);
+    Y[inertIndex].clamp_min(0);
 }

applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/YEqn.H

@@ -38,11 +38,11 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 
             fvOptions.correct(Yi);
 
-            Yi.max(0.0);
+            Yi.clamp_min(0);
             Yt += Yi;
         }
     }
 
     Y[inertIndex] = scalar(1) - Yt;
-    Y[inertIndex].max(0.0);
+    Y[inertIndex].clamp_min(0);
 }

applications/solvers/lagrangian/simpleCoalParcelFoam/YEqn.H

@@ -38,11 +38,11 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 
             fvOptions.correct(Yi);
 
-            Yi.max(0.0);
+            Yi.clamp_min(0);
             Yt += Yi;
         }
     }
 
     Y[inertIndex] = scalar(1) - Yt;
-    Y[inertIndex].max(0.0);
+    Y[inertIndex].clamp_min(0);
 }

applications/solvers/lagrangian/sprayFoam/YEqn.H

@@ -39,11 +39,11 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 
             fvOptions.correct(Yi);
 
-            Yi.max(0.0);
+            Yi.clamp_min(0);
             Yt += Yi;
         }
     }
 
     Y[inertIndex] = scalar(1) - Yt;
-    Y[inertIndex].max(0.0);
+    Y[inertIndex].clamp_min(0);
 }

applications/utilities/preProcessing/applyBoundaryLayer/applyBoundaryLayer.C

@@ -92,7 +92,7 @@ void blendField
         volScalarField fld(fieldHeader, mesh);
         scalarField& pf = fld.primitiveFieldRef();
         pf = (1 - mask)*pf + mask*boundaryLayerField;
-        fld.max(SMALL);
+        fld.clamp_min(SMALL);
 
         // Correct the processor patches only.
         // Do not correct BC
@@ -131,7 +131,7 @@ void calcOmegaField
         scalarField& pf = omega.primitiveFieldRef();
 
         pf = (1 - mask)*pf + mask*epsilonBL/(Cmu*kBL + SMALL);
-        omega.max(SMALL);
+        omega.clamp_min(SMALL);
 
         // Correct the processor patches only.
         // Do not correct BC

applications/utilities/preProcessing/setTurbulenceFields/setTurbulenceFields.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2022 OpenCFD Ltd.
+    Copyright (C) 2022-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -241,13 +241,7 @@ void calcF
     solve(fEqn);
 
     // (M:p. 2)
-    const dimensioned<scalarMinMax> fMinMax
+    f.clamp_range(0, scalar(1) - Foam::exp(-scalar(400)/scalar(50)));
-    (
-        dimless,
-        scalarMinMax(Zero, scalar(1) - Foam::exp(-scalar(400)/scalar(50)))
-    );
-
-    f.clip(fMinMax);
 }
 
 

src/TurbulenceModels/turbulenceModels/LES/dynamicKEqn/dynamicKEqn.C

@@ -102,7 +102,7 @@ volScalarField dynamicKEqn<BasicTurbulenceModel>::Ce() const
     (
         0.5*(filter_(magSqr(this->U_)) - magSqr(filter_(this->U_)))
     );
-    KK.max(dimensionedScalar("small", KK.dimensions(), SMALL));
+    KK.clamp_min(SMALL);
 
     return Ce(D, KK);
 }
@@ -245,7 +245,7 @@ void dynamicKEqn<BasicTurbulenceModel>::correct()
     tgradU.clear();
 
     volScalarField KK(0.5*(filter_(magSqr(U)) - magSqr(filter_(U))));
-    KK.max(dimensionedScalar("small", KK.dimensions(), SMALL));
+    KK.clamp_min(SMALL);
 
     tmp<fvScalarMatrix> kEqn
     (

src/TurbulenceModels/turbulenceModels/ReynoldsStress/ReynoldsStress.C

@@ -39,20 +39,15 @@ void Foam::ReynoldsStress<BasicTurbulenceModel>::boundNormalStress
     volSymmTensorField& R
 ) const
 {
-    scalar kMin = this->kMin_.value();
+    const scalar kMin = this->kMin_.value();
 
-    R.max
+    R.clamp_min
     (
-        dimensionedSymmTensor
+        symmTensor
         (
-            "zero",
+            kMin, -GREAT, -GREAT,
-            R.dimensions(),
+            kMin, -GREAT,
-            symmTensor
+            kMin
-            (
-                kMin, -GREAT, -GREAT,
-                kMin, -GREAT,
-                kMin
-            )
         )
     );
 }

src/combustionModels/FSD/reactionRateFlameAreaModels/relaxation/relaxation.C

@@ -136,8 +136,7 @@ void Foam::reactionRateFlameAreaModels::relaxation::correct
        - fvm::SuSp(rho*(tau + Rc), omega_)
     );
 
-    omega_.min(omega0);
+    omega_.clamp_range(0, omega0.value());
-    omega_.max(0.0);
 }
 
 

src/finiteVolume/fvMesh/wallDist/patchDistMethods/Poisson/PoissonPatchDistMethod.C

@@ -120,7 +120,7 @@ bool Foam::patchDistMethods::Poisson::correct
     // Need to stabilise the y for overset meshes since the holed cells
     // keep the initial value (0.0) so the gradient of that will be
     // zero as well. Turbulence models do not like zero wall distance.
-    y.max(SMALL);
+    y.clamp_min(SMALL);
 
     // For overset: enforce smooth y field (yPsi is smooth, magGradyPsi is
     // not)

src/finiteVolume/fvMesh/wallDist/patchDistMethods/advectionDiffusion/advectionDiffusionPatchDistMethod.C

@@ -153,7 +153,7 @@ bool Foam::patchDistMethods::advectionDiffusion::correct
     // Need to stabilise the y for overset meshes since the holed cells
     // keep the initial value (0.0) so the gradient of that will be
     // zero as well. Turbulence models do not like zero wall distance.
-    y.max(SMALL);
+    y.clamp_min(SMALL);
 
     // Only calculate n if the field is defined
     if (notNull(n))

src/functionObjects/field/comfort/comfort.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2019 OpenFOAM Foundation
-    Copyright (C) 2021 OpenCFD Ltd.
+    Copyright (C) 2021-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -41,6 +41,11 @@ namespace functionObjects
 }
 }
 
+
+// Temperature bounds based on EN ISO 7730 (10 - 40 degC)
+static const Foam::scalarMinMax Tbounds(283.15, 313.15);
+
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 Foam::tmp<Foam::volScalarField> Foam::functionObjects::comfort::magU() const
@@ -93,7 +98,7 @@ Foam::dimensionedScalar Foam::functionObjects::comfort::Trad() const
     }
 
     // Bounds based on EN ISO 7730
-    if ((Trad.value() < 283.15) || (Trad.value() > 313.15))
+    if (!Tbounds.contains(Trad.value()))
     {
         WarningInFunction
             << "The calculated mean wall radiation temperature is out of the\n"
@@ -177,9 +182,6 @@ Foam::tmp<Foam::volScalarField> Foam::functionObjects::comfort::Tcloth
 
     Tcl.storePrevIter();
 
-    // Same temperatures as for the radiation
-    const dimensionedScalar Tmin(dimTemperature, 283.15);
-    const dimensionedScalar Tmax(dimTemperature, 313.15);
 
     // Iterative solving of equation (2)
     do
@@ -208,7 +210,7 @@ Foam::tmp<Foam::volScalarField> Foam::functionObjects::comfort::Tcloth
 
         // Make sure that Tcl is in some physical limit (same range as we used
         // for the radiative estimation - based on ISO EN 7730:2005)
-        Tcl.clip(Tmin, Tmax);
+        Tcl.clamp_range(Tbounds);
 
     } while (!converged(Tcl) && i++ < maxClothIter_);
 
@@ -422,7 +424,7 @@ bool Foam::functionObjects::comfort::execute()
 
     // Limit the velocity field to the values given in EN ISO 7733
     volScalarField Umag(mag(lookupObject<volVectorField>("U")));
-    Umag.clip(Umin, Umax);
+    Umag.clamp_range(Umin, Umax);
 
     // Calculate the turbulent intensity if turbulent kinetic energy field k
     // exists

src/functionObjects/field/limitFields/limitFields.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019 OpenCFD Ltd
+    Copyright (C) 2019-2023 OpenCFD Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -46,9 +46,10 @@ const Foam::Enum
 >
 Foam::functionObjects::limitFields::limitTypeNames_
 ({
-    { limitType::MIN, "min" },
+    { limitType::CLAMP_MIN, "min" },
-    { limitType::MAX, "max" },
+    { limitType::CLAMP_MAX, "max" },
-    { limitType::BOTH, "both" },
+    { limitType::CLAMP_RANGE, "range" },
+    { limitType::CLAMP_RANGE, "both" },
 });
 
 
@@ -67,16 +68,31 @@ bool Foam::functionObjects::limitFields::limitScalarField
 
     auto& field = *fieldPtr;
 
-    if (limit_ & MIN)
+    if (limitType::CLAMP_NONE != withBounds_)
     {
-        Log << ": min(" << gMin(field) << ")";
+        MinMax<scalar> currentRange = gMinMax(field);
-        field.max(dimensionedScalar("", field.dimensions(), min_));
+
+        if (withBounds_ & limitType::CLAMP_MIN)
+        {
+            Log << ": min(" << currentRange.min() << ')';
+        }
+        if (withBounds_ & limitType::CLAMP_MAX)
+        {
+            Log << ": max(" << currentRange.max() << ')';
+        }
     }
 
-    if (limit_ & MAX)
+    if (limitType::CLAMP_MIN == withBounds_)
     {
-        Log << ": max(" << gMax(field) << ")";
+        field.clamp_min(min_);
-        field.min(dimensionedScalar("", field.dimensions(), max_));
+    }
+    else if (limitType::CLAMP_MAX == withBounds_)
+    {
+        field.clamp_max(max_);
+    }
+    else if (limitType::CLAMP_RANGE == withBounds_)
+    {
+        field.clamp_range(min_, max_);
     }
 
     return true;
@@ -93,8 +109,8 @@ Foam::functionObjects::limitFields::limitFields
 )
 :
     fvMeshFunctionObject(name, runTime, dict),
-    limit_(MIN),
     fieldSet_(mesh_),
+    withBounds_(limitType::CLAMP_NONE),
     min_(-VGREAT),
     max_(VGREAT)
 {
@@ -106,19 +122,21 @@ Foam::functionObjects::limitFields::limitFields
 
 bool Foam::functionObjects::limitFields::read(const dictionary& dict)
 {
+    withBounds_ = limitType::CLAMP_NONE;
+
     if (fvMeshFunctionObject::read(dict))
     {
         Info<< type() << " " << name() << ":" << nl;
 
-        limit_ = limitTypeNames_.get("limit", dict);
+        withBounds_ = limitTypeNames_.get("limit", dict);
 
-        if (limit_ & MIN)
+        if (withBounds_ & limitType::CLAMP_MIN)
         {
             min_ = dict.get<scalar>("min");
             Info<< "    Imposing lower limit " << min_ << nl;
         }
 
-        if (limit_ & MAX)
+        if (withBounds_ & limitType::CLAMP_MAX)
         {
             max_ = dict.get<scalar>("max");
             Info<< "    Imposing upper limit " << max_ << nl;

src/functionObjects/field/limitFields/limitFields.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -127,9 +127,10 @@ public:
 
         enum limitType : unsigned
         {
-            MIN = 0x1,          //!< limit by minimum value
+            CLAMP_NONE = 0,     //!< No limit
-            MAX = 0x2,          //!< limit by maximum value
+            CLAMP_MIN = 0x1,    //!< Clamp minimum value
-            BOTH = (MIN | MAX)  //!< limit by both minimum and maximum values
+            CLAMP_MAX = 0x2,    //!< Clamp maximum value
+            CLAMP_RANGE = (CLAMP_MIN | CLAMP_MAX)  //!< Clamp min/max
         };
 
 
@@ -140,12 +141,12 @@ protected:
         //- Limit type names
         static const Enum<limitType> limitTypeNames_;
 
-        //- Limiting type
-        limitType limit_;
-
         //- Fields to limit
         volFieldSelection fieldSet_;
 
+        //- Limiting type
+        limitType withBounds_;
+
         //- Minimum limit
         scalar min_;
 

src/functionObjects/field/limitFields/limitFieldsTemplates.C

@@ -47,23 +47,23 @@ bool Foam::functionObjects::limitFields::limitField(const word& fieldName)
 
     const dimensionedScalar eps("eps", field.dimensions(), ROOTVSMALL);
 
-    if (limit_ & MIN)
+    if (withBounds_ & limitType::CLAMP_MIN)
     {
         volScalarField mField(typeName + ":mag" + field.name(), mag(field));
         Log << " min(|" << gMin(mField) << "|)";
         //field.normalise();
         field /= mag(field) + eps;
-        mField.max(dimensionedScalar("min", field.dimensions(), min_));
+        mField.clamp_min(min_);
         field *= mField;
     }
 
-    if (limit_ & MAX)
+    if (withBounds_ & limitType::CLAMP_MAX)
     {
         volScalarField mField(typeName + ":mag" + field.name(), mag(field));
         Log << " max(|" << gMax(mField) << "|)";
         //field.normalise();
         field /= mag(field) + eps;
-        mField.min(dimensionedScalar("max", field.dimensions(), max_));
+        mField.clamp_max(max_);
         field *= mField;
     }
 

src/functionObjects/field/stabilityBlendingFactor/stabilityBlendingFactor.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2021 OpenCFD Ltd.
+    Copyright (C) 2018-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -448,8 +448,7 @@ bool Foam::functionObjects::stabilityBlendingFactor::init(bool first)
     }
 
     indicator.correctBoundaryConditions();
-    indicator.min(1.0);
+    indicator.clamp_range(zero_one{});
-    indicator.max(0.0);
 
     // Update the blended surface field
     auto& surBlended = mesh_.lookupObjectRef<surfaceScalarField>(resultName_);

src/phaseSystemModels/multiphaseInter/phasesSystem/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C

@@ -428,7 +428,7 @@ void Foam::MultiComponentPhaseModel<BasePhaseModel, phaseThermo>::solveYi
     }
 
     X_[inertIndex_] = scalar(1) - Yt;
-    X_[inertIndex_].max(0.0);
+    X_[inertIndex_].clamp_min(0);
 
     calculateMassFractions();
 }

src/phaseSystemModels/reactingEuler/multiphaseSystem/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C

@@ -425,8 +425,8 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctInterfaceThermo()
         volScalarField H2(heatTransferModelIter().second()->K());
 
         // Limit the H[12] to avoid /0
-        H1.max(SMALL);
+        H1.clamp_min(SMALL);
-        H2.max(SMALL);
+        H2.clamp_min(SMALL);
 
         Tf = (H1*T1 + H2*T2 + iDmdtNew*L)/(H1 + H2);
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/PhaseSystems/TwoResistanceHeatTransferPhaseSystem/TwoResistanceHeatTransferPhaseSystem.C

@@ -297,8 +297,8 @@ correctInterfaceThermo()
         );
 
         // Limit the H[12] to avoid /0
-        H1.max(SMALL);
+        H1.clamp_min(SMALL);
-        H2.max(SMALL);
+        H2.clamp_min(SMALL);
 
         Tf = (H1*T1 + H2*T2 + dmdt*L)/(H1 + H2);
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/interfacialModels/dragModels/IshiiZuber/IshiiZuber.C

@@ -86,7 +86,7 @@ Foam::dragModels::IshiiZuber::CdRe() const
     );
 
     volScalarField F((muc/muMix)*sqrt(1 - pair_.dispersed()));
-    F.max(1e-3);
+    F.clamp_min(1e-3);
 
     const volScalarField Ealpha((1 + 17.67*pow(F, 0.8571428))/(18.67*F));
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/interfacialModels/liftModels/Moraga/Moraga.C

@@ -87,11 +87,8 @@ Foam::tmp<Foam::volScalarField> Foam::liftModels::Moraga::Cl() const
             << endl;
     }
 
-    Re.min(1200.0);
+    Re.clamp_range(1200.0, 18800.0);
-    Re.max(18800.0);
+    sqrSr.clamp_range(0.0016, 0.04);
-
-    sqrSr.min(0.0016);
-    sqrSr.max(0.04);
 
     return 0.2*exp(- Re*sqrSr/3.6e5 - 0.12)*exp(Re*sqrSr/3.0e7);
 }

src/phaseSystemModels/reactingEuler/multiphaseSystem/multiphaseSystem/multiphaseSystem.C

@@ -708,7 +708,7 @@ void Foam::multiphaseSystem::solve()
         phase.alphaRhoPhiRef() =
             fvc::interpolate(phase.rho())*phase.alphaPhi();
 
-        phase.clip(SMALL, 1 - SMALL);
+        phase.clamp_range(SMALL, 1 - SMALL);
     }
 
     calcAlphas();

src/phaseSystemModels/reactingEuler/multiphaseSystem/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C

@@ -124,14 +124,14 @@ void Foam::MultiComponentPhaseModel<BasePhaseModel>::correctThermo()
     if (inertIndex_ != -1)
     {
         Yi[inertIndex_] = scalar(1) - Yt;
-        Yi[inertIndex_].max(0);
+        Yi[inertIndex_].clamp_min(0);
     }
     else
     {
         forAll(Yi, i)
         {
             Yi[i] /= Yt;
-            Yi[i].max(0);
+            Yi[i].clamp_min(0);
         }
     }
 

src/phaseSystemModels/reactingEuler/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C

@@ -526,8 +526,7 @@ void Foam::RASModels::kineticTheoryModel::correct()
         kappa_ = conductivityModel_->kappa(alpha, Theta_, gs0_, rho, da, e_);
     }
 
-    Theta_.max(0);
+    Theta_.clamp_range(0, 100);
-    Theta_.min(100);
 
     {
         // particle viscosity (Table 3.2, p.47)
@@ -559,7 +558,8 @@ void Foam::RASModels::kineticTheoryModel::correct()
         );
 
         // Limit viscosity and add frictional viscosity
-        nut_.min(maxNut_);
+        nut_.clamp_max(maxNut_);
+
         nuFric_ = min(nuFric_, maxNut_ - nut_);
         nut_ += nuFric_;
     }

src/phaseSystemModels/reactingEuler/twoPhaseSystem/twoPhaseSystem.C

@@ -340,7 +340,7 @@ void Foam::twoPhaseSystem::solve()
             << endl;
 
         // Ensure the phase-fractions are bounded
-        alpha1.clip(SMALL, 1 - SMALL);
+        alpha1.clamp_range(SMALL, 1 - SMALL);
 
         // Update the phase-fraction of the other phase
         alpha2 = scalar(1) - alpha1;

src/phaseSystemModels/twoPhaseEuler/phaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C

@@ -530,8 +530,7 @@ void Foam::RASModels::kineticTheoryModel::correct()
         kappa_ = conductivityModel_->kappa(alpha, Theta_, gs0_, rho, da, e_);
     }
 
-    Theta_.max(0);
+    Theta_.clamp_range(0, 100);
-    Theta_.min(100);
 
     {
         // particle viscosity (Table 3.2, p.47)
@@ -563,7 +562,8 @@ void Foam::RASModels::kineticTheoryModel::correct()
         );
 
         // Limit viscosity and add frictional viscosity
-        nut_.min(maxNut_);
+        nut_.clamp_max(maxNut_);
+
         nuFric_ = min(nuFric_, maxNut_ - nut_);
         nut_ += nuFric_;
     }

src/phaseSystemModels/twoPhaseEuler/twoPhaseSystem/interfacialModels/dragModels/IshiiZuber/IshiiZuber.C

@@ -87,7 +87,7 @@ Foam::dragModels::IshiiZuber::CdRe() const
     );
 
     volScalarField F((muc/muMix)*sqrt(1 - pair_.dispersed()));
-    F.max(1e-3);
+    F.clamp_min(1e-3);
 
     volScalarField Ealpha((1 + 17.67*pow(F, 0.8571428))/(18.67*F));
 

src/phaseSystemModels/twoPhaseEuler/twoPhaseSystem/interfacialModels/liftModels/Moraga/Moraga.C

@@ -87,11 +87,8 @@ Foam::tmp<Foam::volScalarField> Foam::liftModels::Moraga::Cl() const
             << endl;
     }
 
-    Re.min(1200.0);
+    Re.clamp_range(1200.0, 18800.0);
-    Re.max(18800.0);
+    sqrSr.clamp_range(0.0016, 0.04);
-
-    sqrSr.min(0.0016);
-    sqrSr.max(0.04);
 
     return 0.2*exp(- Re*sqrSr/3.6e5 - 0.12)*exp(Re*sqrSr/3.0e7);
 }

src/phaseSystemModels/twoPhaseEuler/twoPhaseSystem/twoPhaseSystem/twoPhaseSystem.C

@@ -530,8 +530,7 @@ void Foam::twoPhaseSystem::solve()
             << endl;
 
         // Ensure the phase-fractions are bounded
-        alpha1.max(0);
+        alpha1.clamp_range(zero_one{});
-        alpha1.min(1);
 
         alpha2 = scalar(1) - alpha1;
     }

src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.C

@@ -128,7 +128,7 @@ tmp<areaScalarField> filmTurbulenceModel::Cw() const
             const scalar h0 = film_.h0().value();
 
             Cw.primitiveFieldRef() = 3*mu/((h + h0)*rho);
-            Cw.min(5000.0);
+            Cw.clamp_max(5000.0);
 
             break;
         }

src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C

@@ -276,7 +276,7 @@ void reactingOneDim::solveSpeciesMass()
         }
 
         YiEqn.solve(regionMesh().solver("Yi"));
-        Yi.max(0.0);
+        Yi.clamp_min(0);
         Yt += Yi;
     }
 

src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C

@@ -412,7 +412,7 @@ void kinematicSingleLayer::solveThickness
     }
 
     // Bound film thickness by a minimum of zero
-    delta_.max(0.0);
+    delta_.clamp_min(0);
 
     // Update U field
     U_ -= fvc::reconstruct(deltarUAf*phiAdd);

src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayerI.H

@@ -250,7 +250,7 @@ inline tmp<volScalarField> kinematicSingleLayer::gNormClipped() const
     );
 
     volScalarField& gNormClipped = tgNormClipped.ref();
-    gNormClipped.min(0.0);
+    gNormClipped.clamp_max(0);
 
     return tgNormClipped;
 }

src/regionModels/surfaceFilmModels/submodels/kinematic/filmTurbulenceModel/laminar/laminar.C

@@ -136,7 +136,7 @@ tmp<fvVectorMatrix> laminar::Su(volVectorField& U) const
     // employ simple coeff-based model
     volScalarField Cs("Cs", Cf_*rhop*mag(Up - U));
     volScalarField Cw("Cw", mu/((1.0/3.0)*(delta + film.deltaSmall())));
-    Cw.min(5000.0);
+    Cw.clamp_max(5000.0);
 
     return
     (

src/regionModels/surfaceFilmModels/submodels/thermo/filmViscosityModel/thixotropicViscosity/thixotropicViscosity.C

@@ -85,8 +85,7 @@ thixotropicViscosity::thixotropicViscosity
         film.regionMesh()
     )
 {
-    lambda_.min(1);
+    lambda_.clamp_range(zero_one{});
-    lambda_.max(0);
 
     // Initialise viscosity to inf value because it cannot be evaluated yet
     mu_ = muInf_;
@@ -166,8 +165,7 @@ void thixotropicViscosity::correct
     lambdaEqn.relax();
     lambdaEqn.solve();
 
-    lambda_.min(1);
+    lambda_.clamp_range(zero_one{});
-    lambda_.max(0);
 
     mu_ = muInf_/(sqr(1 - K_*lambda_) + ROOTVSMALL);
     mu_.correctBoundaryConditions();

src/regionModels/surfaceFilmModels/submodels/thermo/phaseChangeModel/waxSolventEvaporation/waxSolventEvaporation.C

@@ -357,8 +357,7 @@ void waxSolventEvaporation::correctModel
         YsolventEqn.relax();
         YsolventEqn.solve();
 
-        Ysolvent_.min(1);
+        Ysolvent_.clamp_range(zero_one{});
-        Ysolvent_.max(0);
 
         scalarField dm
         (

src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2017-2020 OpenCFD Ltd.
+    Copyright (C) 2017-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -478,18 +478,23 @@ thermoSingleLayer::thermoSingleLayer
     ),
     phaseChange_(phaseChangeModel::New(*this, coeffs())),
     radiation_(filmRadiationModel::New(*this, coeffs())),
-    Tmin_(-VGREAT),
+    withTbounds_(limitType::CLAMP_NONE),
-    Tmax_(VGREAT)
+    Tbounds_(0, 5000)
 {
-    if (coeffs().readIfPresent("Tmin", Tmin_))
+    unsigned userLimits(limitType::CLAMP_NONE);
+
+    if (coeffs().readIfPresent("Tmin", Tbounds_.min()))
     {
-        Info<< "    limiting minimum temperature to " << Tmin_ << endl;
+        userLimits |= limitType::CLAMP_MIN;
+        Info<< "    limiting minimum temperature to " << Tbounds_.min() << nl;
     }
 
-    if (coeffs().readIfPresent("Tmax", Tmax_))
+    if (coeffs().readIfPresent("Tmax", Tbounds_.max()))
     {
-        Info<< "    limiting maximum temperature to " << Tmax_ << endl;
+        userLimits |= limitType::CLAMP_MAX;
+        Info<< "    limiting maximum temperature to " << Tbounds_.max() << nl;
     }
+    withTbounds_ = limitType(userLimits);
 
     if (thermo_.hasMultiComponentCarrier())
     {

src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -40,8 +41,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef thermoSingleLayer_H
+#ifndef Foam_thermoSingleLayer_H
-#define thermoSingleLayer_H
+#define Foam_thermoSingleLayer_H
 
 #include "kinematicSingleLayer.H"
 #include "SLGThermo.H"
@@ -69,7 +70,7 @@ class thermoSingleLayer
 :
     public kinematicSingleLayer
 {
-    // Private member functions
+    // Private Member Functions
 
         //- No copy construct
         thermoSingleLayer(const thermoSingleLayer&) = delete;
@@ -83,7 +84,17 @@ class thermoSingleLayer
 
 protected:
 
-    // Protected data
+    // Protected Data
+
+        //- Enumerated limiter type
+        enum limitType : unsigned
+        {
+            CLAMP_NONE = 0,     //!< No limit
+            CLAMP_MIN = 0x1,    //!< Clamp minimum value
+            CLAMP_MAX = 0x2,    //!< Clamp maximum value
+            CLAMP_RANGE = (CLAMP_MIN | CLAMP_MAX)  //!< Clamp min/max
+        };
+
 
         // Thermo properties
 
@@ -184,14 +195,14 @@ protected:
 
         // Limits
 
-            //- Minimum temperature limit (optional)
+            //- Limiting type
-            scalar Tmin_;
+            limitType withTbounds_;
 
-            //- Maximum temperature limit (optional)
+            //- Temperature limits (optional)
-            scalar Tmax_;
+            scalarMinMax Tbounds_;
 
 
-    // Protected member functions
+    // Protected Member Functions
 
         //- Read control parameters from dictionary
         virtual bool read();

src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayerI.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -102,8 +103,18 @@ inline tmp<volScalarField> thermoSingleLayer::T
         )
     );
 
-    tT.ref().min(Tmax_);
+    if (limitType::CLAMP_MIN == withTbounds_)
-    tT.ref().max(Tmin_);
+    {
+        tT.ref().clamp_min(Tbounds_.min());
+    }
+    else if (limitType::CLAMP_MAX == withTbounds_)
+    {
+        tT.ref().clamp_max(Tbounds_.max());
+    }
+    else if (limitType::CLAMP_RANGE == withTbounds_)
+    {
+        tT.ref().clamp_range(Tbounds_);
+    }
 
     return tT;
 }

src/thermophysicalModels/chemistryModel/functionObjects/BilgerMixtureFraction/BilgerMixtureFraction.C

@@ -89,11 +89,7 @@ void Foam::functionObjects::BilgerMixtureFraction::calcBilgerMixtureFraction()
            /thermo_.W(i);
     }
     f_Bilger /= o2RequiredFuelOx_;
-    f_Bilger.clip
+    f_Bilger.clamp_range(zero_one{});
-    (
-        dimensionedScalar(dimless, 0),
-        dimensionedScalar(dimless, 1)
-    );
 }