diff --git a/applications/solvers/combustion/reactingFoam/setRDeltaT.H b/applications/solvers/combustion/reactingFoam/setRDeltaT.H
index 6bf4cf8375f4db1054d32aa70556e2656ec9bd8d..7b449fd90077b8bf292d6b18e1c1956e7816d5df 100644
--- a/applications/solvers/combustion/reactingFoam/setRDeltaT.H
+++ b/applications/solvers/combustion/reactingFoam/setRDeltaT.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2016 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020,2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -57,10 +57,22 @@ License
     // (relative to reference value)
     scalar alphaY(pimpleDict.getOrDefault<scalar>("alphaY", 1.0));
 
-    Info<< "Time scales min/max:" << endl;
 
-    // Cache old reciprocal time scale field
-    volScalarField rDeltaT0("rDeltaT0", rDeltaT);
+    // The old reciprocal time scale field, with any damping factor
+    tmp<volScalarField> rDeltaT0_damped;
+
+    // Calculate damped value before applying any other changes
+    if
+    (
+        rDeltaTDampingCoeff < 1
+     && runTime.timeIndex() > runTime.startTimeIndex() + 1
+    )
+    {
+        rDeltaT0_damped = (scalar(1) - rDeltaTDampingCoeff)*(rDeltaT);
+    }
+
+
+    Info<< "Time scales min/max:" << endl;
 
     // Flow time scale
     {
@@ -70,8 +82,8 @@ License
            /((2*maxCo)*mesh.V()*rho())
         );
 
-        // Limit the largest time scale
-        rDeltaT.max(1/maxDeltaT);
+        // Limit the largest time scale (=> smallest reciprocal time)
+        rDeltaT.clamp_min(1/maxDeltaT);
 
         Info<< "    Flow        = "
             << 1/gMax(rDeltaT.primitiveField()) << ", "
@@ -86,11 +98,11 @@ License
             mag(Qdot)/(alphaTemp*rho*thermo.Cp()*T)
         );
 
+        rDeltaT.primitiveFieldRef().clamp_min(rDeltaTT);
+
         Info<< "    Temperature = "
             << 1/(gMax(rDeltaTT.field()) + VSMALL) << ", "
             << 1/(gMin(rDeltaTT.field()) + VSMALL) << endl;
-
-        rDeltaT.ref() = max(rDeltaT(), rDeltaTT);
     }
 
     // Reaction rate time scale
@@ -138,11 +150,11 @@ License
 
         if (foundY)
         {
+            rDeltaT.primitiveFieldRef().clamp_min(rDeltaTY);
+
             Info<< "    Composition = "
                 << 1/(gMax(rDeltaTY.field()) + VSMALL) << ", "
                 << 1/(gMin(rDeltaTY.field()) + VSMALL) << endl;
-
-            rDeltaT.ref() = max(rDeltaT(), rDeltaTY);
         }
         else
         {
@@ -161,20 +173,12 @@ License
         fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
     }
 
-    // Limit rate of change of time scale
+    // Limit rate of change of time scale (=> smallest reciprocal time)
     // - reduce as much as required
     // - only increase at a fraction of old time scale
-    if
-    (
-        rDeltaTDampingCoeff < 1
-     && runTime.timeIndex() > runTime.startTimeIndex() + 1
-    )
+    if (rDeltaT0_damped)
     {
-        rDeltaT = max
-        (
-            rDeltaT,
-            (scalar(1) - rDeltaTDampingCoeff)*rDeltaT0
-        );
+        rDeltaT.clamp_min(rDeltaT0_damped());
     }
 
     // Update tho boundary values of the reciprocal time-step
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
index fe3d7428164790fd1f20427da40daab7730920f3..d4dd47300ec87197c5056d06c8bb322b55571b45 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
@@ -78,8 +78,8 @@
     }
 
     rho = thermo.rho();
-    rho = max(rho, rhoMin[i]);
-    rho = min(rho, rhoMax[i]);
+
+    rho.clamp_range(rhoMin[i], rhoMax[i]);
     rho.relax();
 
     Info<< "Min/max rho:" << min(rho).value() << ' '
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/setRDeltaT.H b/applications/solvers/lagrangian/coalChemistryFoam/setRDeltaT.H
index bda1cc2d8837ef033a07444d808a33cdf816167e..aabfd0e0156b66fd84f4188d3ff8cea260694980 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/setRDeltaT.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/setRDeltaT.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020,2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -54,10 +54,21 @@ License
     scalar alphaTemp(pimpleDict.getOrDefault("alphaTemp", 0.05));
 
 
-    Info<< "Time scales min/max:" << endl;
+    // The old reciprocal time scale field, with any damping factor
+    tmp<volScalarField> rDeltaT0_damped;
+
+    // Calculate damped value before applying any other changes
+    if
+    (
+        rDeltaTDampingCoeff < 1
+     && runTime.timeIndex() > runTime.startTimeIndex() + 1
+    )
+    {
+        rDeltaT0_damped = (scalar(1) - rDeltaTDampingCoeff)*(rDeltaT);
+    }
 
-    // Cache old reciprocal time scale field
-    volScalarField rDeltaT0("rDeltaT0", rDeltaT);
+
+    Info<< "Time scales min/max:" << endl;
 
     // Flow time scale
     {
@@ -67,8 +78,8 @@ License
            /((2*maxCo)*mesh.V()*rho())
         );
 
-        // Limit the largest time scale
-        rDeltaT.max(1/maxDeltaT);
+        // Limit the largest time scale (=> smallest reciprocal time)
+        rDeltaT.clamp_min(1/maxDeltaT);
 
         Info<< "    Flow        = "
             << gMin(1/rDeltaT.primitiveField()) << ", "
@@ -93,15 +104,11 @@ License
            )
         );
 
+        rDeltaT.primitiveFieldRef().clamp_min(rDeltaTT);
+
         Info<< "    Temperature = "
             << gMin(1/(rDeltaTT.field() + VSMALL)) << ", "
             << gMax(1/(rDeltaTT.field() + VSMALL)) << endl;
-
-        rDeltaT.ref() = max
-        (
-            rDeltaT(),
-            rDeltaTT
-        );
     }
 
     // Update tho boundary values of the reciprocal time-step
@@ -113,20 +120,12 @@ License
         fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
     }
 
-    // Limit rate of change of time scale
+    // Limit rate of change of time scale (=> smallest reciprocal time)
     // - reduce as much as required
     // - only increase at a fraction of old time scale
-    if
-    (
-        rDeltaTDampingCoeff < 1.0
-     && runTime.timeIndex() > runTime.startTimeIndex() + 1
-    )
+    if (rDeltaT0_damped)
     {
-        rDeltaT = max
-        (
-            rDeltaT,
-            (scalar(1) - rDeltaTDampingCoeff)*rDeltaT0
-        );
+        rDeltaT.clamp_min(rDeltaT0_damped());
     }
 
     Info<< "    Overall     = "
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/setRDeltaT.H b/applications/solvers/lagrangian/reactingParcelFoam/setRDeltaT.H
index 5977b18fe67dde0922adfaf598f3fa510ef72233..8ef4cd1eb2317085bd0d837df2a0fa8cc0f510af 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/setRDeltaT.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/setRDeltaT.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020,2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -54,10 +54,21 @@ License
     scalar alphaTemp(pimpleDict.getOrDefault("alphaTemp", 0.05));
 
 
-    Info<< "Time scales min/max:" << endl;
+    // The old reciprocal time scale field, with any damping factor
+    tmp<volScalarField> rDeltaT0_damped;
+
+    // Calculate damped value before applying any other changes
+    if
+    (
+        rDeltaTDampingCoeff < 1
+     && runTime.timeIndex() > runTime.startTimeIndex() + 1
+    )
+    {
+        rDeltaT0_damped = (scalar(1) - rDeltaTDampingCoeff)*(rDeltaT);
+    }
+
 
-    // Cache old reciprocal time scale field
-    volScalarField rDeltaT0("rDeltaT0", rDeltaT);
+    Info<< "Time scales min/max:" << endl;
 
     // Flow time scale
     {
@@ -67,8 +78,8 @@ License
            /((2*maxCo)*mesh.V()*rho())
         );
 
-        // Limit the largest time scale
-        rDeltaT.max(1/maxDeltaT);
+        // Limit the largest time scale (=> smallest reciprocal time)
+        rDeltaT.clamp_min(1/maxDeltaT);
 
         Info<< "    Flow        = "
             << gMin(1/rDeltaT.primitiveField()) << ", "
@@ -92,15 +103,11 @@ License
            )
         );
 
+        rDeltaT.primitiveFieldRef().clamp_min(rDeltaTT);
+
         Info<< "    Temperature = "
             << gMin(1/(rDeltaTT.field() + VSMALL)) << ", "
             << gMax(1/(rDeltaTT.field() + VSMALL)) << endl;
-
-        rDeltaT.ref() = max
-        (
-            rDeltaT(),
-            rDeltaTT
-        );
     }
 
     // Update the boundary values of the reciprocal time-step
@@ -112,22 +119,17 @@ License
         fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
     }
 
-    // Limit rate of change of time scale
+    // Limit rate of change of time scale (=> smallest reciprocal time)
     // - reduce as much as required
     // - only increase at a fraction of old time scale
-    if
-    (
-        rDeltaTDampingCoeff < 1.0
-     && runTime.timeIndex() > runTime.startTimeIndex() + 1
-    )
+    if (rDeltaT0_damped)
     {
-        rDeltaT = max
-        (
-            rDeltaT,
-            (scalar(1) - rDeltaTDampingCoeff)*rDeltaT0
-        );
+        rDeltaT.clamp_min(rDeltaT0_damped());
     }
 
+    // Update the boundary values of the reciprocal time-step
+    rDeltaT.correctBoundaryConditions();
+
     Info<< "    Overall     = "
         << gMin(1/rDeltaT.primitiveField())
         << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/pEqn.H
index 1ce5db0ec8f1fbd4758d126be362ebb0dd2d3354..2b3d264c62610cadd097c79d19d39408e876818a 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/pEqn.H
@@ -48,8 +48,7 @@ U.correctBoundaryConditions();
 fvOptions.correct(U);
 
 rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
+rho.clamp_range(rhoMin, rhoMax);
 rho.relax();
 
 Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;
diff --git a/applications/solvers/lagrangian/simpleCoalParcelFoam/pEqn.H b/applications/solvers/lagrangian/simpleCoalParcelFoam/pEqn.H
index 055eff6f0b34ec7e30fb3da5d8ac9f6f96386d25..1c85828e24a632853aa7272e79ed5cf744e4e0f3 100644
--- a/applications/solvers/lagrangian/simpleCoalParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/simpleCoalParcelFoam/pEqn.H
@@ -49,8 +49,7 @@
     fvOptions.correct(U);
 
     rho = thermo.rho();
-    rho = max(rho, rhoMin);
-    rho = min(rho, rhoMax);
+    rho.clamp_range(rhoMin, rhoMax);
     rho.relax();
 
     Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;
diff --git a/applications/solvers/lagrangian/sprayFoam/pEqn.H b/applications/solvers/lagrangian/sprayFoam/pEqn.H
index 198bf00d52caab6aae55d8e2075ddacae20bfa94..fbdcc577908d41ea50663918e6b777722503de7b 100644
--- a/applications/solvers/lagrangian/sprayFoam/pEqn.H
+++ b/applications/solvers/lagrangian/sprayFoam/pEqn.H
@@ -1,6 +1,5 @@
 rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
+rho.clamp_range(rhoMin, rhoMax);
 rho.relax();
 
 volScalarField rAU(1.0/UEqn.A());
@@ -94,8 +93,7 @@ p.relax();
 
 // Recalculate density from the relaxed pressure
 rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
+rho.clamp_range(rhoMin, rhoMax);
 rho.relax();
 Info<< "rho min/max : " << min(rho).value() << " " << max(rho).value() << endl;
 
diff --git a/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/pEqn.H b/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/pEqn.H
index a7850a5e903257f9a921d475fc91edb2148ca694..92bd8c2564b6cf0a77515eb7de0d7afcca7a6318 100644
--- a/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/pEqn.H
+++ b/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/pEqn.H
@@ -1,6 +1,5 @@
 rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
+rho.clamp_range(rhoMin, rhoMax);
 rho.relax();
 
 volScalarField rAU(1.0/UEqn.A());
@@ -94,8 +93,7 @@ p.relax();
 
 // Recalculate density from the relaxed pressure
 rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
+rho.clamp_range(rhoMin, rhoMax);
 rho.relax();
 Info<< "rho min/max : " << min(rho).value() << " " << max(rho).value() << endl;
 
diff --git a/applications/solvers/multiphase/VoF/setRDeltaT.H b/applications/solvers/multiphase/VoF/setRDeltaT.H
index 313fca2e33b1d9162125a7f32d770fd15df93ff2..b690f0accc69e416b7fc384784539517e4701f08 100644
--- a/applications/solvers/multiphase/VoF/setRDeltaT.H
+++ b/applications/solvers/multiphase/VoF/setRDeltaT.H
@@ -53,6 +53,21 @@
         pimpleDict.getOrDefault<scalar>("maxDeltaT", GREAT)
     );
 
+
+    // The old reciprocal time scale field, with any damping factor
+    tmp<volScalarField> rDeltaT0_damped;
+
+    // Calculate damped value before applying any other changes
+    if
+    (
+        rDeltaTDampingCoeff < 1
+     && runTime.timeIndex() > runTime.startTimeIndex() + 1
+    )
+    {
+        rDeltaT0_damped = (scalar(1) - rDeltaTDampingCoeff)*(rDeltaT);
+    }
+
+
     volScalarField rDeltaT0("rDeltaT0", rDeltaT);
 
     // Set the reciprocal time-step from the local Courant number
@@ -114,20 +129,12 @@
         << gMin(1/rDeltaT.primitiveField())
         << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
 
-    // Limit rate of change of time scale
+    // Limit rate of change of time scale (=> smallest reciprocal time)
     // - reduce as much as required
     // - only increase at a fraction of old time scale
-    if
-    (
-        rDeltaTDampingCoeff < 1.0
-     && runTime.timeIndex() > runTime.startTimeIndex() + 1
-    )
+    if (rDeltaT0_damped)
     {
-        rDeltaT = max
-        (
-            rDeltaT,
-            (scalar(1) - rDeltaTDampingCoeff)*rDeltaT0
-        );
+        rDeltaT.clamp_min(rDeltaT0_damped());
 
         Info<< "Damped flow time scale min/max = "
             << gMin(1/rDeltaT.primitiveField())
diff --git a/applications/test/minMax2/Test-minMax2.cxx b/applications/test/minMax2/Test-minMax2.cxx
index e228f0dce2cc0415876f398bce5e95bbe75bd1f7..b58363c74485f3514049703a8202fb918c57a3ec 100644
--- a/applications/test/minMax2/Test-minMax2.cxx
+++ b/applications/test/minMax2/Test-minMax2.cxx
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2023 OpenCFD Ltd.
+    Copyright (C) 2019-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -175,6 +175,36 @@ int main(int argc, char *argv[])
             << nl;
     }
 
+    {
+        scalarField field1(10);
+        scalarField field2(10);
+        scalarField work;
+
+        Random rnd(4567);
+        for (scalar& val : field1)
+        {
+            val = rnd.position(scalar(-0.2), scalar(1.2));
+        }
+        for (scalar& val : field2)
+        {
+            val = rnd.position(scalar(-0.1), scalar(1.1));
+        }
+
+        Info<< nl
+            << "field1: " << flatOutput(field1) << nl
+            << "field2: " << flatOutput(field2) << nl;
+
+        work = field1;
+        work.clamp_min(field2);
+
+        Info<< "clamp_min: " << flatOutput(work) << nl;
+
+        work = field1;
+        work.clamp_max(field2);
+
+        Info<< "clamp_max: " << flatOutput(work) << nl;
+    }
+
     Info<< nl << "\nDone\n" << endl;
     return 0;
 }
diff --git a/src/OpenFOAM/fields/FieldFields/FieldField/FieldField.C b/src/OpenFOAM/fields/FieldFields/FieldField/FieldField.C
index 17d07fb2f069ccb5bd79f6c2bd16ea2a7277a5cb..1cd6a54c6b2d67c444d87b9b2de88605036f1ff0 100644
--- a/src/OpenFOAM/fields/FieldFields/FieldField/FieldField.C
+++ b/src/OpenFOAM/fields/FieldFields/FieldField/FieldField.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2023 OpenCFD Ltd.
+    Copyright (C) 2019-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -308,6 +308,36 @@ void FieldField<Field, Type>::clamp_max
 }
 
 
+template<template<class> class Field, class Type>
+void FieldField<Field, Type>::clamp_min
+(
+    const FieldField<Field, Type>& lower
+)
+{
+    const label loopLen = this->size();
+
+    for (label i = 0; i < loopLen; ++i)
+    {
+        (*this)[i].clamp_min(lower[i]);
+    }
+}
+
+
+template<template<class> class Field, class Type>
+void FieldField<Field, Type>::clamp_max
+(
+    const FieldField<Field, Type>& upper
+)
+{
+    const label loopLen = this->size();
+
+    for (label i = 0; i < loopLen; ++i)
+    {
+        (*this)[i].clamp_max(upper[i]);
+    }
+}
+
+
 template<template<class> class Field, class Type>
 void FieldField<Field, Type>::clamp_range
 (
diff --git a/src/OpenFOAM/fields/FieldFields/FieldField/FieldField.H b/src/OpenFOAM/fields/FieldFields/FieldField/FieldField.H
index d8dbbcdf0e04f56a4146034a8e96d8389350455a..f1c9fda937f710a6d9561aed34bd4c97183dfb96 100644
--- a/src/OpenFOAM/fields/FieldFields/FieldField/FieldField.H
+++ b/src/OpenFOAM/fields/FieldFields/FieldField/FieldField.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2022-2023 OpenCFD Ltd.
+    Copyright (C) 2022-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -152,9 +152,15 @@ public:
         //- Impose lower (floor) clamp on the field values (in-place)
         void clamp_min(const Type& lower);
 
+        //- Impose lower (floor) clamp on the field values (in-place)
+        void clamp_min(const FieldField<Field, Type>& lower);
+
         //- Impose upper (ceiling) clamp on the field values (in-place)
         void clamp_max(const Type& upper);
 
+        //- Impose upper (ceiling) clamp on the field values (in-place)
+        void clamp_max(const FieldField<Field, Type>& upper);
+
         //- Clamp field values (in-place) to the specified range.
         //  Does not check if range is valid or not.
         void clamp_range(const Type& lower, const Type& upper);
diff --git a/src/OpenFOAM/fields/Fields/Field/Field.C b/src/OpenFOAM/fields/Fields/Field/Field.C
index 0372523102a2411f043ceaead0e67f6ee004341d..1b500ff6c97cbd1211444a285c10fe27ae88b05b 100644
--- a/src/OpenFOAM/fields/Fields/Field/Field.C
+++ b/src/OpenFOAM/fields/Fields/Field/Field.C
@@ -669,6 +669,7 @@ void Foam::Field<Type>::clamp_min(const Type& lower)
     }
 }
 
+
 template<class Type>
 void Foam::Field<Type>::clamp_max(const Type& upper)
 {
@@ -681,6 +682,40 @@ void Foam::Field<Type>::clamp_max(const Type& upper)
 }
 
 
+template<class Type>
+void Foam::Field<Type>::clamp_min(const UList<Type>& lower)
+{
+    // Use free function max() [sic] to impose component-wise clamp_min
+    std::transform
+    (
+        // Can use (std::execution::par_unseq | std::execution::unseq)
+        this->begin(),
+        // this->end() but with some extra range safety
+        this->begin(lower.size()),
+        lower.begin(),
+        this->begin(),
+        maxOp<Type>()
+    );
+}
+
+
+template<class Type>
+void Foam::Field<Type>::clamp_max(const UList<Type>& upper)
+{
+    // Use free function min() [sic] to impose component-wise clamp_max
+    std::transform
+    (
+        // Can use (std::execution::par_unseq | std::execution::unseq)
+        this->begin(),
+        // this->end() but with some extra range safety
+        this->begin(upper.size()),
+        upper.begin(),
+        this->begin(),
+        minOp<Type>()
+    );
+}
+
+
 template<class Type>
 void Foam::Field<Type>::clamp_range(const Type& lower, const Type& upper)
 {
diff --git a/src/OpenFOAM/fields/Fields/Field/Field.H b/src/OpenFOAM/fields/Fields/Field/Field.H
index fb7851edc6dfc78953a95c81c413dd79d35b44fa..3441a50e0b5949f53be62a8a2ffbf717804afd3b 100644
--- a/src/OpenFOAM/fields/Fields/Field/Field.H
+++ b/src/OpenFOAM/fields/Fields/Field/Field.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2023 OpenCFD Ltd.
+    Copyright (C) 2015-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -50,6 +50,7 @@ SourceFiles
 #include "direction.H"
 #include "labelList.H"
 #include "keyType.H"
+#include "ops.H"
 #include "scalarList.H"
 #include "VectorSpace.H"
 #include "IOobjectOption.H"
@@ -441,9 +442,15 @@ public:
         //- Impose lower (floor) clamp on the field values (in-place)
         void clamp_min(const Type& lower);
 
+        //- Impose lower (floor) clamp on the field values (in-place)
+        void clamp_min(const UList<Type>& lower);
+
         //- Impose upper (ceiling) clamp on the field values (in-place)
         void clamp_max(const Type& upper);
 
+        //- Impose upper (ceiling) clamp on the field values (in-place)
+        void clamp_max(const UList<Type>& upper);
+
         //- Clamp field values (in-place) to the specified range.
         //  Does not check if range is valid or not.
         void clamp_range(const Type& lower, const Type& upper);
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C
index 2ef9c028adee99b0137cb4e6692c72f0cee0ba30..7667b784355bb41f04d586792849a6810d9ac5cf 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C
@@ -1276,6 +1276,38 @@ void Foam::GeometricField<Type, PatchField, GeoMesh>::clamp_max
 }
 
 
+template<class Type, template<class> class PatchField, class GeoMesh>
+void Foam::GeometricField<Type, PatchField, GeoMesh>::clamp_min
+(
+    const GeometricField<Type, PatchField, GeoMesh>& lower
+)
+{
+    primitiveFieldRef().clamp_min(lower.primitiveField());
+    boundaryFieldRef().clamp_min(lower.boundaryField());
+    correctLocalBoundaryConditions();
+    if (GeometricBoundaryField<Type, PatchField, GeoMesh>::debug)
+    {
+        boundaryField().check();
+    }
+}
+
+
+template<class Type, template<class> class PatchField, class GeoMesh>
+void Foam::GeometricField<Type, PatchField, GeoMesh>::clamp_max
+(
+    const GeometricField<Type, PatchField, GeoMesh>& upper
+)
+{
+    primitiveFieldRef().clamp_max(upper.primitiveField());
+    boundaryFieldRef().clamp_max(upper.boundaryField());
+    correctLocalBoundaryConditions();
+    if (GeometricBoundaryField<Type, PatchField, GeoMesh>::debug)
+    {
+        boundaryField().check();
+    }
+}
+
+
 template<class Type, template<class> class PatchField, class GeoMesh>
 void Foam::GeometricField<Type, PatchField, GeoMesh>::clamp_range
 (
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H
index 5ff69df5cab49a8c63d6823a0d030bf04ecf81a2..2345159762d3334e190c7cb087fd6d824df6adae 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H
@@ -936,17 +936,25 @@ public:
         //- Impose lower (floor) clamp on the field values (in-place)
         void clamp_min(const Type& lower);
 
-        //- Impose upper (ceiling) clamp on the field values (in-place)
-        void clamp_max(const Type& upper);
-
         //- Impose lower (floor) clamp on the field values (in-place)
         //  No dimension checking
         void clamp_min(const dimensioned<Type>& lower);
 
+        //- Impose lower (floor) clamp on the field values (in-place)
+        //  No dimension checking
+        void clamp_min(const GeometricField<Type, PatchField, GeoMesh>& lower);
+
+        //- Impose upper (ceiling) clamp on the field values (in-place)
+        void clamp_max(const Type& upper);
+
         //- Impose upper (ceiling) clamp on the field values (in-place)
         //  No dimension checking
         void clamp_max(const dimensioned<Type>& upper);
 
+        //- Impose upper (ceiling) clamp on the field values (in-place)
+        //  No dimension checking
+        void clamp_max(const GeometricField<Type, PatchField, GeoMesh>& upper);
+
         //- Clamp field values (in-place) to the specified range.
         //  Does not check if range is valid or not. No dimension checking.
         void clamp_range(const dimensioned<MinMax<Type>>& range);
@@ -1055,11 +1063,13 @@ public:
         //- Use minimum of the field and specified value. ie, clamp_max().
         //  This sets the \em ceiling on the field values
         //  \deprecated(2023-01) prefer clamp_max()
+        FOAM_DEPRECATED_STRICT(2023-01, "clamp_max() method")
         void min(const dimensioned<Type>& upper) { this->clamp_max(upper); }
 
         //- Use maximum of the field and specified value. ie, clamp_min().
         //  This sets the \em floor on the field values
         //  \deprecated(2023-01) prefer clamp_min()
+        FOAM_DEPRECATED_STRICT(2023-01, "clamp_min() method")
         void max(const dimensioned<Type>& lower) { this->clamp_min(lower); }
 
         //- Deprecated(2019-01) identical to clamp_range()
diff --git a/src/finiteArea/faMatrices/faScalarMatrix/faScalarMatrix.C b/src/finiteArea/faMatrices/faScalarMatrix/faScalarMatrix.C
index f9d39a7d00d1f713efa73ffb39d7fd75c0db3cbe..d944b639285ff1a39a3379e559d293d56bcc5036 100644
--- a/src/finiteArea/faMatrices/faScalarMatrix/faScalarMatrix.C
+++ b/src/finiteArea/faMatrices/faScalarMatrix/faScalarMatrix.C
@@ -85,7 +85,7 @@ Foam::solverPerformance Foam::faMatrix<Foam::scalar>::solve
         internalCoeffs_,
         psi_.boundaryField().scalarInterfaces(),
         solverControls
-    )->solve(psi.ref(), totalSource);
+    )->solve(psi.primitiveFieldRef(), totalSource);
 
     if (logLevel)
     {
@@ -146,7 +146,7 @@ Foam::tmp<Foam::areaScalarField> Foam::faMatrix<Foam::scalar>::H() const
     Hphi.primitiveFieldRef() = (lduMatrix::H(psi_.primitiveField()) + source_);
     addBoundarySource(Hphi.primitiveFieldRef());
 
-    Hphi.ref() /= psi_.mesh().S();
+    Hphi.primitiveFieldRef() /= psi_.mesh().S();
     Hphi.correctBoundaryConditions();
 
     return tHphi;
diff --git a/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C b/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C
index 074cccfde8f83d418f8ad06a8e7b567fabc306b2..9ad65f87485bfd9df30af8e075c0737bfb3f6ba9 100644
--- a/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C
+++ b/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C
@@ -256,8 +256,7 @@ void Foam::rhoThermo::correctRho
 )
 {
     rho_ += deltaRho;
-    rho_ = max(rho_, rhoMin);
-    rho_ = min(rho_, rhoMax);
+    rho_.clamp_range(rhoMin, rhoMax);
 }
 
 void Foam::rhoThermo::correctRho(const Foam::volScalarField& deltaRho)