From 01727c84f13e567e35d813245cc7a50892ff1aa4 Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Wed, 2 Apr 2025 12:01:28 +0200
Subject: [PATCH] ENH: use gMinMax() instead of separate gMin(), gMax()

- for reciprocal values, gMinMax() first and then calculate the
  reciprocal, which avoids creating temporaries

STYLE: prefer MinMax to separate min/max accounting

COMP: namespace qualify min/max for deltaT, CourantNo, etc (#3348)
---
 .../solvers/combustion/PDRFoam/setDeltaT.H    |   6 +-
 .../solvers/combustion/chemFoam/setDeltaT.H   |   3 +-
 .../fireFoam/setMultiRegionDeltaT.H           |  13 +-
 .../combustion/reactingFoam/setRDeltaT.H      |  24 ++--
 .../compressible/rhoCentralFoam/setRDeltaT.H  |  10 +-
 .../compressible/rhoPimpleFoam/setRDeltaT.H   |  26 ++--
 .../fluid/compressibleMultiRegionCourantNo.H  |   2 +-
 ...woPhaseRadCoupledMixedFvPatchScalarField.C |  35 +++--
 .../fluid/compressibleMultiRegionCourantNo.H  |   6 +-
 .../fluid/compressibleMultiRegionCourantNo.H  |   4 +-
 .../include/setInitialMultiRegionDeltaT.H     |   6 +-
 .../include/setMultiRegionDeltaT.H            |  12 +-
 .../solid/solidRegionDiffusionNo.H            |   2 +-
 .../heatTransfer/solidFoam/solidDiffusionNo.H |   2 +-
 .../incompressible/pimpleFoam/setRDeltaT.H    |  26 ++--
 .../lagrangian/coalChemistryFoam/setRDeltaT.H |  18 ++-
 .../reactingParcelFoam/setMultiRegionDeltaT.H |  11 +-
 .../reactingParcelFoam/setRDeltaT.H           |  18 ++-
 .../solvers/multiphase/VoF/setDeltaT.H        |   7 +-
 .../solvers/multiphase/VoF/setRDeltaT.H       |  25 ++--
 .../multiphase/cavitatingFoam/setDeltaT.H     |   7 +-
 .../cavitatingFoam/setInitialDeltaT.H         |   8 +-
 .../setDeltaT.H                               |  12 +-
 .../reactingMultiphaseEulerFoam/setRDeltaT.H  |   6 +-
 .../reactingTwoPhaseEulerFoam/CourantNos.H    |   2 +-
 .../reactingTwoPhaseEulerFoam/setRDeltaT.H    |   6 +-
 .../multiphase/twoPhaseEulerFoam/CourantNos.H |   2 +-
 .../test/fieldMapping/Test-fieldMapping.C     |  27 ++--
 applications/test/hexRef8/Test-hexRef8.C      |  27 ++--
 .../advanced/refineHexMesh/refineHexMesh.C    |  27 ++--
 .../mesh/advanced/selectCells/edgeStats.C     | 121 +++++++++--------
 .../snappyRefineMesh/snappyRefineMesh.C       | 117 ++++++++--------
 .../mesh/manipulation/refineMesh/refineMesh.C | 122 +++++++----------
 .../createBoxTurb/createBoxTurb.C             |  14 +-
 .../setAlphaField/setAlphaField.C             |  13 +-
 src/combustionModels/EDC/EDC.C                |   5 +-
 .../interfaceTrackingFvMesh.C                 |   7 +-
 .../motionSmoother/motionSmootherAlgo.C       |  17 ++-
 .../componentDisplacementMotionSolver.C       |   4 +-
 .../displacementInterpolationMotionSolver.C   |  23 ++--
 .../displacementLayeredMotionMotionSolver.C   |  24 ++--
 .../edgeLimitedFaGrad/edgeLimitedFaGrads.C    |  18 ++-
 .../faceLimitedFaGrad/faceLimitedFaGrads.C    |  18 ++-
 .../cfdTools/general/include/setDeltaT.H      |   6 +-
 .../general/include/setInitialDeltaT.H        |   4 +-
 ...lectrostaticDepositionFvPatchScalarField.C |  13 +-
 .../mappedField/mappedFieldFvPatchField.C     |   9 +-
 .../mappedMixedFieldFvPatchField.C            |   9 +-
 .../mappedFixedValueFvPatchField.C            |   9 +-
 .../mappedMixed/mappedMixedFvPatchField.C     |   9 +-
 .../timeVaryingMappedFixedValueFvPatchField.C |   9 +-
 .../turbulentDFSEMInletFvPatchVectorField.C   |   4 +-
 .../waveSurfacePressureFvPatchScalarField.C   |   8 +-
 .../ddtSchemes/ddtScheme/ddtScheme.C          |  16 +--
 .../cellLimitedGrad/cellLimitedGrad.C         |   9 +-
 .../faceLimitedGrad/faceLimitedGrads.C        |  18 ++-
 .../averageNeighbourFvGeometryScheme.C        |   8 +-
 .../highAspectRatio/cellAspectRatio.C         |   7 +-
 .../highAspectRatioFvGeometryScheme.C         |   7 +-
 .../field/AMIWeights/AMIWeights.C             | 128 +++++++++---------
 src/functionObjects/field/yPlus/yPlus.C       |  17 ++-
 ...meVaryingMappedFixedValuePointPatchField.C |   9 +-
 .../effectivenessTable/effectivenessTable.C   |  16 +--
 .../Templates/KinematicCloud/KinematicCloud.C |  11 +-
 .../clouds/Templates/MPPICCloud/MPPICCloud.C  |  11 +-
 .../snappyHexMeshDriver/snappyLayerDriver.C   |   3 +-
 .../snappyHexMeshDriver/snappySnapDriver.C    |   6 +-
 .../AMIInterpolation/AMIInterpolation.C       |  42 +++---
 .../cyclicPeriodicAMIPolyPatch.C              |  35 +++--
 .../PatchFunction1/MappedFile/MappedFile.C    |  10 +-
 .../regularisationPDE/Helmoltz/Helmholtz.C    |   2 +-
 .../ThermalPhaseChangePhaseSystem.C           |  33 +++--
 ...allBoilingWallFunctionFvPatchScalarField.C |  33 ++---
 ...ixedMultiPhaseHeatFluxFvPatchScalarField.C |   7 +-
 .../reactingOneDim/reactingOneDim.C           |   7 +-
 .../thermoSingleLayer/thermoSingleLayer.C     |   7 +-
 ...pedMassWallTemperatureFvPatchScalarField.C |   9 +-
 .../thermalBaffle1DFvPatchScalarField.C       |  10 +-
 ...tureCoupledBaffleMixedFvPatchScalarField.C |   9 +-
 .../enthalpySorptionFvPatchScalarField.C      |   8 +-
 .../speciesSorptionFvPatchScalarField.C       |   8 +-
 .../movingConeTopoFvMesh.C                    |   4 +-
 .../isoAdvection/isoAdvection.C               |   8 +-
 .../isoAdvection/isoAdvectionTemplates.C      |  22 +--
 84 files changed, 837 insertions(+), 646 deletions(-)

diff --git a/applications/solvers/combustion/PDRFoam/setDeltaT.H b/applications/solvers/combustion/PDRFoam/setDeltaT.H
index 365126ce90..177fc1fa7f 100644
--- a/applications/solvers/combustion/PDRFoam/setDeltaT.H
+++ b/applications/solvers/combustion/PDRFoam/setDeltaT.H
@@ -36,11 +36,13 @@ Description
 if (adjustTimeStep)
 {
     scalar maxDeltaTFact = maxCo/(CoNum + StCoNum + SMALL);
-    scalar deltaTFact = min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
+
+    const scalar deltaTFact =
+        Foam::min(Foam::min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
 
     runTime.setDeltaT
     (
-        min
+        Foam::min
         (
             deltaTFact*runTime.deltaTValue(),
             maxDeltaT
diff --git a/applications/solvers/combustion/chemFoam/setDeltaT.H b/applications/solvers/combustion/chemFoam/setDeltaT.H
index ef50a97026..93f68ceb6c 100644
--- a/applications/solvers/combustion/chemFoam/setDeltaT.H
+++ b/applications/solvers/combustion/chemFoam/setDeltaT.H
@@ -1,5 +1,6 @@
 if (adjustTimeStep)
 {
-    runTime.setDeltaT(min(dtChem, maxDeltaT));
+    runTime.setDeltaT(Foam::min(dtChem, maxDeltaT));
+
     Info<< "deltaT = " <<  runTime.deltaTValue() << endl;
 }
diff --git a/applications/solvers/combustion/fireFoam/setMultiRegionDeltaT.H b/applications/solvers/combustion/fireFoam/setMultiRegionDeltaT.H
index 31c9e675fc..a5ac0fcc9a 100644
--- a/applications/solvers/combustion/fireFoam/setMultiRegionDeltaT.H
+++ b/applications/solvers/combustion/fireFoam/setMultiRegionDeltaT.H
@@ -54,9 +54,18 @@ if (adjustTimeStep)
 
     runTime.setDeltaT
     (
-        min
+        Foam::min
         (
-            dt0*min(min(TFactorFluid, min(TFactorFilm, TFactorSolid)), 1.2),
+            dt0
+          * Foam::min
+            (
+                Foam::min
+                (
+                    TFactorFluid,
+                    Foam::min(TFactorFilm, TFactorSolid)
+                ),
+                1.2
+            ),
             maxDeltaT
         )
     );
diff --git a/applications/solvers/combustion/reactingFoam/setRDeltaT.H b/applications/solvers/combustion/reactingFoam/setRDeltaT.H
index 7b449fd900..67a2b00f91 100644
--- a/applications/solvers/combustion/reactingFoam/setRDeltaT.H
+++ b/applications/solvers/combustion/reactingFoam/setRDeltaT.H
@@ -85,9 +85,11 @@ License
         // Limit the largest time scale (=> smallest reciprocal time)
         rDeltaT.clamp_min(1/maxDeltaT);
 
+        auto limits = gMinMax(rDeltaT.primitiveField());
+        limits.reset(1/(limits.max()+VSMALL), 1/(limits.min()+VSMALL));
+
         Info<< "    Flow        = "
-            << 1/gMax(rDeltaT.primitiveField()) << ", "
-            << 1/gMin(rDeltaT.primitiveField()) << endl;
+            << limits.min() << ", " << limits.max() << endl;
     }
 
     // Heat release rate time scale
@@ -100,9 +102,11 @@ License
 
         rDeltaT.primitiveFieldRef().clamp_min(rDeltaTT);
 
+        auto limits = gMinMax(rDeltaTT.field());
+        limits.reset(1/(limits.max()+VSMALL), 1/(limits.min()+VSMALL));
+
         Info<< "    Temperature = "
-            << 1/(gMax(rDeltaTT.field()) + VSMALL) << ", "
-            << 1/(gMin(rDeltaTT.field()) + VSMALL) << endl;
+            << limits.min() << ", " << limits.max() << endl;
     }
 
     // Reaction rate time scale
@@ -152,9 +156,11 @@ License
         {
             rDeltaT.primitiveFieldRef().clamp_min(rDeltaTY);
 
+            auto limits = gMinMax(rDeltaTY.field());
+            limits.reset(1/(limits.max()+VSMALL), 1/(limits.min()+VSMALL));
+
             Info<< "    Composition = "
-                << 1/(gMax(rDeltaTY.field()) + VSMALL) << ", "
-                << 1/(gMin(rDeltaTY.field()) + VSMALL) << endl;
+                << limits.min() << ", " << limits.max() << endl;
         }
         else
         {
@@ -184,9 +190,11 @@ License
     // Update tho boundary values of the reciprocal time-step
     rDeltaT.correctBoundaryConditions();
 
+    auto limits = gMinMax(rDeltaT.field());
+    limits.reset(1/(limits.max()+VSMALL), 1/(limits.min()+VSMALL));
+
     Info<< "    Overall     = "
-        << 1/gMax(rDeltaT.primitiveField())
-        << ", " << 1/gMin(rDeltaT.primitiveField()) << endl;
+        << limits.min() << ", " << limits.max() << endl;
 }
 
 
diff --git a/applications/solvers/compressible/rhoCentralFoam/setRDeltaT.H b/applications/solvers/compressible/rhoCentralFoam/setRDeltaT.H
index 5e2a2b302c..8f9560b742 100644
--- a/applications/solvers/compressible/rhoCentralFoam/setRDeltaT.H
+++ b/applications/solvers/compressible/rhoCentralFoam/setRDeltaT.H
@@ -23,7 +23,11 @@
 
     fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
 
-    Info<< "Flow time scale min/max = "
-        << gMin(1/rDeltaT.primitiveField())
-        << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
+    {
+        auto limits = gMinMax(rDeltaT.primitiveField());
+        limits.reset(1/(limits.max()+VSMALL), 1/(limits.min()+VSMALL));
+
+        Info<< "Flow time scale min/max = "
+            << limits.min() << ", " << limits.max() << endl;
+    }
 }
diff --git a/applications/solvers/compressible/rhoPimpleFoam/setRDeltaT.H b/applications/solvers/compressible/rhoPimpleFoam/setRDeltaT.H
index e818402eff..51ecac9f07 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/setRDeltaT.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/setRDeltaT.H
@@ -52,18 +52,26 @@
     // Update the boundary values of the reciprocal time-step
     rDeltaT.correctBoundaryConditions();
 
-    Info<< "Flow time scale min/max = "
-        << gMin(1/rDeltaT.primitiveField())
-        << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
+    {
+        auto limits = gMinMax(rDeltaT.primitiveField());
+        limits.reset(1/(limits.max()+VSMALL), 1/(limits.min()+VSMALL));
+
+        Info<< "Flow time scale min/max = "
+            << limits.min() << ", " << limits.max() << endl;
+    }
 
     if (rDeltaTSmoothingCoeff < 1.0)
     {
         fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
     }
 
-    Info<< "Smoothed flow time scale min/max = "
-        << gMin(1/rDeltaT.primitiveField())
-        << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
+    {
+        auto limits = gMinMax(rDeltaT.primitiveField());
+        limits.reset(1/(limits.max()+VSMALL), 1/(limits.min()+VSMALL));
+
+        Info<< "Smoothed flow time scale min/max = "
+            << limits.min() << ", " << limits.max() << endl;
+    }
 
     // Limit rate of change of time scale
     // - reduce as much as required
@@ -78,8 +86,10 @@
             rDeltaT0
            *max(rDeltaT/rDeltaT0, scalar(1) - rDeltaTDampingCoeff);
 
+        auto limits = gMinMax(rDeltaT.primitiveField());
+        limits.reset(1/(limits.max()+VSMALL), 1/(limits.min()+VSMALL));
+
         Info<< "Damped flow time scale min/max = "
-            << gMin(1/rDeltaT.primitiveField())
-            << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
+            << limits.min() << ", " << limits.max() << endl;
     }
 }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/compressibleMultiRegionCourantNo.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/compressibleMultiRegionCourantNo.H
index 3ca2f68581..5ccc00d3ad 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/compressibleMultiRegionCourantNo.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/compressibleMultiRegionCourantNo.H
@@ -1,7 +1,7 @@
     scalar CoNum = -GREAT;
     forAll(fluidRegions, regionI)
     {
-        CoNum = max
+        CoNum = Foam::max
         (
             compressibleCourantNo
             (
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/derivedFvPatchFields/turbulentTemperatureTwoPhaseRadCoupledMixed/turbulentTemperatureTwoPhaseRadCoupledMixedFvPatchScalarField.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/derivedFvPatchFields/turbulentTemperatureTwoPhaseRadCoupledMixed/turbulentTemperatureTwoPhaseRadCoupledMixedFvPatchScalarField.C
index 33d1c05548..a0232c0b1c 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/derivedFvPatchFields/turbulentTemperatureTwoPhaseRadCoupledMixed/turbulentTemperatureTwoPhaseRadCoupledMixedFvPatchScalarField.C
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/derivedFvPatchFields/turbulentTemperatureTwoPhaseRadCoupledMixed/turbulentTemperatureTwoPhaseRadCoupledMixedFvPatchScalarField.C
@@ -387,15 +387,18 @@ updateCoeffs()
         {
             scalar Q = gSum(kappa(Tp)*patch().magSf()*snGrad());
 
-                Info<< "T solid : " << nl << endl;
+            auto limits = gMinMax(Tp);
+            auto avg = gAverage(Tp);
 
-                Info
-                    << " heat transfer rate from solid:" << Q
-                    << " walltemperature "
-                    << " min:" << gMin(Tp)
-                    << " max:" << gMax(Tp)
-                    << " avg:" << gAverage(Tp) << nl
-                    << endl;
+            Info<< "T solid : " << nl << endl;
+
+            Info
+                << " heat transfer rate from solid:" << Q
+                << " walltemperature "
+                << " min:" << limits.min()
+                << " max:" << limits.max()
+                << " avg:" << avg << nl
+                << endl;
         }
     }
     else if (regionType_ == fluid)
@@ -445,10 +448,16 @@ updateCoeffs()
             scalarField qLiq((Tp - Tc)*KdeltaLiq);
             scalarField qVap((Tp - Tv.patchInternalField())*KdeltaVap);
 
+            auto infoT = gMinMax(Tp);
+            auto avgT = gAverage(Tp);
+
+            auto infoLiq = gMinMax(qLiq);
+            auto infoVap = gMinMax(qVap);
+
             Info<< "T flow : " << nl << endl;
 
-            Info<< "  qLiq: " << gMin(qLiq) << " - " << gMax(qLiq) << endl;
-            Info<< "  qVap: " << gMin(qVap) << " - " << gMax(qVap) << endl;
+            Info<< "  qLiq: " << infoLiq.min() << " - " << infoLiq.max() << nl
+                << "  qVap: " << infoVap.min() << " - " << infoVap.max() << nl;
 
             scalar QLiq = gSum(qLiq*patch().magSf());
             scalar QVap = gSum(qVap*patch().magSf());
@@ -457,9 +466,9 @@ updateCoeffs()
             Info<<  " Heat transfer to Vap: " << QVap << endl;
 
             Info<< " walltemperature "
-                << " min:" << gMin(Tp)
-                << " max:" << gMax(Tp)
-                << " avg:" << gAverage(Tp)
+                << " min:" << infoT.min()
+                << " max:" << infoT.max()
+                << " avg:" << avgT
                 << endl;
         }
     }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/compressibleMultiRegionCourantNo.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/compressibleMultiRegionCourantNo.H
index 290bb596cf..392d74f6a3 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/compressibleMultiRegionCourantNo.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/compressibleMultiRegionCourantNo.H
@@ -31,7 +31,7 @@
             );
 
 
-            CoNum =
+            scalar regionCoNum =
                 0.5*gMax
                 (
                     sumPhi/fluidRegions[regioni].V().field()
@@ -41,9 +41,9 @@
             (
                 fvc::surfaceSum(mag(phi1 - phi2))().primitiveField()
               / fluidRegions[regioni].V().field()
-            )*runTime.deltaTValue(),
+            )*runTime.deltaTValue();
 
-            CoNum = max(UrCoNum, CoNum);
+            CoNum = Foam::max(CoNum, Foam::max(regionCoNum, UrCoNum));
         }
     }
 
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H
index dd129047dc..d5b5c7ceaa 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H
@@ -2,7 +2,7 @@
 
     forAll(fluidRegions, regioni)
     {
-        CoNum = max
+        CoNum = Foam::max
         (
             compressibleCourantNo
             (
@@ -17,7 +17,7 @@
 /*
     forAll(porousFluidRegions, porousi)
     {
-        CoNum = max
+        CoNum = Foam::max
         (
             compressibleCourantNo
             (
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/include/setInitialMultiRegionDeltaT.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/include/setInitialMultiRegionDeltaT.H
index 09f4a717a9..8be6345765 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/include/setInitialMultiRegionDeltaT.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/include/setInitialMultiRegionDeltaT.H
@@ -47,10 +47,10 @@ if (adjustTimeStep)
 
         runTime.setDeltaT
         (
-            min
+            Foam::min
             (
-                min(maxCo/CoNum, maxDi/DiNum)*runTime.deltaTValue(),
-                min(runTime.deltaTValue(), maxDeltaT)
+                Foam::min(maxCo/CoNum, maxDi/DiNum)*runTime.deltaTValue(),
+                Foam::min(runTime.deltaTValue(), maxDeltaT)
             )
         );
         Info<< "deltaT = " <<  runTime.deltaTValue() << endl;
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/include/setMultiRegionDeltaT.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/include/setMultiRegionDeltaT.H
index 8eed82dd5c..d81d374f41 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/include/setMultiRegionDeltaT.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/include/setMultiRegionDeltaT.H
@@ -48,18 +48,14 @@ if (adjustTimeStep)
     scalar maxDeltaTFluid = maxCo/(CoNum + SMALL);
     scalar maxDeltaTSolid = maxDi/(DiNum + SMALL);
 
-    scalar deltaTFluid =
-        min
-        (
-            min(maxDeltaTFluid, 1.0 + 0.1*maxDeltaTFluid),
-            1.2
-        );
+    const scalar deltaTFluid =
+        Foam::min(Foam::min(maxDeltaTFluid, 1.0 + 0.1*maxDeltaTFluid), 1.2);
 
     runTime.setDeltaT
     (
-        min
+        Foam::min
         (
-            min(deltaTFluid, maxDeltaTSolid)*runTime.deltaTValue(),
+            Foam::min(deltaTFluid, maxDeltaTSolid)*runTime.deltaTValue(),
             maxDeltaT
         )
     );
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffusionNo.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffusionNo.H
index 6e87fb9203..8c6111d276 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffusionNo.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffusionNo.H
@@ -22,7 +22,7 @@ forAll(solidRegions, i)
     tmp<volScalarField> trho = thermo.rho();
     const volScalarField& rho = trho();
 
-    DiNum = max
+    DiNum = Foam::max
     (
         solidRegionDiffNo
         (
diff --git a/applications/solvers/heatTransfer/solidFoam/solidDiffusionNo.H b/applications/solvers/heatTransfer/solidFoam/solidDiffusionNo.H
index 90317f67c6..3484c35b9f 100644
--- a/applications/solvers/heatTransfer/solidFoam/solidDiffusionNo.H
+++ b/applications/solvers/heatTransfer/solidFoam/solidDiffusionNo.H
@@ -17,7 +17,7 @@ scalar DiNum = -GREAT;
     tmp<volScalarField> trho = thermo.rho();
     const volScalarField& rho = trho();
 
-    DiNum = max
+    DiNum = Foam::max
     (
         solidRegionDiffNo
         (
diff --git a/applications/solvers/incompressible/pimpleFoam/setRDeltaT.H b/applications/solvers/incompressible/pimpleFoam/setRDeltaT.H
index 962bfa8766..3acc01c5fa 100644
--- a/applications/solvers/incompressible/pimpleFoam/setRDeltaT.H
+++ b/applications/solvers/incompressible/pimpleFoam/setRDeltaT.H
@@ -36,18 +36,26 @@
     // Update the boundary values of the reciprocal time-step
     rDeltaT.correctBoundaryConditions();
 
-    Info<< "Flow time scale min/max = "
-        << gMin(1/rDeltaT.primitiveField())
-        << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
+    {
+        auto limits = gMinMax(rDeltaT.primitiveField());
+        limits.reset(1/(limits.max()+VSMALL), 1/(limits.min()+VSMALL));
+
+        Info<< "Flow time scale min/max = "
+            << limits.min() << ", " << limits.max() << endl;
+    }
 
     if (rDeltaTSmoothingCoeff < 1.0)
     {
         fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
     }
 
-    Info<< "Smoothed flow time scale min/max = "
-        << gMin(1/rDeltaT.primitiveField())
-        << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
+    {
+        auto limits = gMinMax(rDeltaT.primitiveField());
+        limits.reset(1/(limits.max()+VSMALL), 1/(limits.min()+VSMALL));
+
+        Info<< "Smoothed flow time scale min/max = "
+            << limits.min() << ", " << limits.max() << endl;
+    }
 
     // Limit rate of change of time scale
     // - reduce as much as required
@@ -62,8 +70,10 @@
             rDeltaT0
            *max(rDeltaT/rDeltaT0, scalar(1) - rDeltaTDampingCoeff);
 
+        auto limits = gMinMax(rDeltaT.primitiveField());
+        limits.reset(1/(limits.max()+VSMALL), 1/(limits.min()+VSMALL));
+
         Info<< "Damped flow time scale min/max = "
-            << gMin(1/rDeltaT.primitiveField())
-            << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
+            << limits.min() << ", " << limits.max() << endl;
     }
 }
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/setRDeltaT.H b/applications/solvers/lagrangian/coalChemistryFoam/setRDeltaT.H
index aabfd0e015..745249cfdf 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/setRDeltaT.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/setRDeltaT.H
@@ -81,9 +81,11 @@ License
         // Limit the largest time scale (=> smallest reciprocal time)
         rDeltaT.clamp_min(1/maxDeltaT);
 
+        auto limits = gMinMax(rDeltaT.primitiveField());
+        limits.reset(1/(limits.max()+VSMALL), 1/(limits.min()+VSMALL));
+
         Info<< "    Flow        = "
-            << gMin(1/rDeltaT.primitiveField()) << ", "
-            << gMax(1/rDeltaT.primitiveField()) << endl;
+            << limits.min() << ", " << limits.max() << endl;
     }
 
     // Reaction source time scale
@@ -106,9 +108,11 @@ License
 
         rDeltaT.primitiveFieldRef().clamp_min(rDeltaTT);
 
+        auto limits = gMinMax(rDeltaTT.field());
+        limits.reset(1/(limits.max()+VSMALL), 1/(limits.min()+VSMALL));
+
         Info<< "    Temperature = "
-            << gMin(1/(rDeltaTT.field() + VSMALL)) << ", "
-            << gMax(1/(rDeltaTT.field() + VSMALL)) << endl;
+            << limits.min() << ", " << limits.max() << endl;
     }
 
     // Update tho boundary values of the reciprocal time-step
@@ -128,9 +132,11 @@ License
         rDeltaT.clamp_min(rDeltaT0_damped());
     }
 
+    auto limits = gMinMax(rDeltaT.primitiveField());
+    limits.reset(1/(limits.max()+VSMALL), 1/(limits.min()+VSMALL));
+
     Info<< "    Overall     = "
-        << gMin(1/rDeltaT.primitiveField())
-        << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
+        << limits.min() << ", " << limits.max() << endl;
 }
 
 
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/setMultiRegionDeltaT.H b/applications/solvers/lagrangian/reactingParcelFoam/setMultiRegionDeltaT.H
index 70d55cb8c8..1257447f3e 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/setMultiRegionDeltaT.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/setMultiRegionDeltaT.H
@@ -36,13 +36,18 @@ Description
 if (adjustTimeStep)
 {
     const scalar maxDeltaTFact =
-        min(maxCo/(CoNum + SMALL), maxCo/(surfaceFilm.CourantNumber() + SMALL));
+        Foam::min
+        (
+            maxCo/(CoNum + SMALL),
+            maxCo/(surfaceFilm.CourantNumber() + SMALL)
+        );
+
     const scalar deltaTFact =
-        min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
+        Foam::min(Foam::min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
 
     runTime.setDeltaT
     (
-        min
+        Foam::min
         (
             deltaTFact*runTime.deltaTValue(),
             maxDeltaT
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/setRDeltaT.H b/applications/solvers/lagrangian/reactingParcelFoam/setRDeltaT.H
index 8ef4cd1eb2..42e187ed3c 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/setRDeltaT.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/setRDeltaT.H
@@ -81,9 +81,11 @@ License
         // Limit the largest time scale (=> smallest reciprocal time)
         rDeltaT.clamp_min(1/maxDeltaT);
 
+        auto limits = gMinMax(rDeltaT.primitiveField());
+        limits.reset(1/(limits.max()+VSMALL), 1/(limits.min()+VSMALL));
+
         Info<< "    Flow        = "
-            << gMin(1/rDeltaT.primitiveField()) << ", "
-            << gMax(1/rDeltaT.primitiveField()) << endl;
+            << limits.min() << ", " << limits.max() << endl;
     }
 
     // Reaction source time scale
@@ -105,9 +107,11 @@ License
 
         rDeltaT.primitiveFieldRef().clamp_min(rDeltaTT);
 
+        auto limits = gMinMax(rDeltaTT.field());
+        limits.reset(1/(limits.max()+VSMALL), 1/(limits.min()+VSMALL));
+
         Info<< "    Temperature = "
-            << gMin(1/(rDeltaTT.field() + VSMALL)) << ", "
-            << gMax(1/(rDeltaTT.field() + VSMALL)) << endl;
+            << limits.min() << ", " << limits.max() << endl;
     }
 
     // Update the boundary values of the reciprocal time-step
@@ -130,9 +134,11 @@ License
     // Update the boundary values of the reciprocal time-step
     rDeltaT.correctBoundaryConditions();
 
+    auto limits = gMinMax(rDeltaT.primitiveField());
+    limits.reset(1/(limits.max()+VSMALL), 1/(limits.min()+VSMALL));
+
     Info<< "    Overall     = "
-        << gMin(1/rDeltaT.primitiveField())
-        << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
+        << limits.min() << ", " << limits.max() << endl;
 }
 
 
diff --git a/applications/solvers/multiphase/VoF/setDeltaT.H b/applications/solvers/multiphase/VoF/setDeltaT.H
index ac59647e7e..79fbe90c4f 100644
--- a/applications/solvers/multiphase/VoF/setDeltaT.H
+++ b/applications/solvers/multiphase/VoF/setDeltaT.H
@@ -36,13 +36,14 @@ Description
 if (adjustTimeStep)
 {
     scalar maxDeltaTFact =
-        min(maxCo/(CoNum + SMALL), maxAlphaCo/(alphaCoNum + SMALL));
+        Foam::min(maxCo/(CoNum + SMALL), maxAlphaCo/(alphaCoNum + SMALL));
 
-    scalar deltaTFact = min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
+    const scalar deltaTFact =
+        Foam::min(Foam::min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
 
     runTime.setDeltaT
     (
-        min
+        Foam::min
         (
             deltaTFact*runTime.deltaTValue(),
             maxDeltaT
diff --git a/applications/solvers/multiphase/VoF/setRDeltaT.H b/applications/solvers/multiphase/VoF/setRDeltaT.H
index b690f0accc..bcc6e3b56d 100644
--- a/applications/solvers/multiphase/VoF/setRDeltaT.H
+++ b/applications/solvers/multiphase/VoF/setRDeltaT.H
@@ -98,10 +98,13 @@
     // Update tho boundary values of the reciprocal time-step
     rDeltaT.correctBoundaryConditions();
 
-    Info<< "Flow time scale min/max = "
-        << gMin(1/rDeltaT.primitiveField())
-        << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
+    {
+        auto limits = gMinMax(rDeltaT.primitiveField());
+        limits.reset(1/(limits.max()+VSMALL), 1/(limits.min()+VSMALL));
 
+        Info<< "Flow time scale min/max = "
+            << limits.min() << ", " << limits.max() << endl;
+    }
     if (rDeltaTSmoothingCoeff < 1.0)
     {
         fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
@@ -125,9 +128,13 @@
         fvc::sweep(rDeltaT, alpha1, nAlphaSweepIter, alphaSpreadDiff);
     }
 
-    Info<< "Smoothed flow time scale min/max = "
-        << gMin(1/rDeltaT.primitiveField())
-        << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
+    {
+        auto limits = gMinMax(rDeltaT.primitiveField());
+        limits.reset(1/(limits.max()+VSMALL), 1/(limits.min()+VSMALL));
+
+        Info<< "Smoothed flow time scale min/max = "
+            << limits.min() << ", " << limits.max() << endl;
+    }
 
     // Limit rate of change of time scale (=> smallest reciprocal time)
     // - reduce as much as required
@@ -136,8 +143,10 @@
     {
         rDeltaT.clamp_min(rDeltaT0_damped());
 
+        auto limits = gMinMax(rDeltaT.primitiveField());
+        limits.reset(1/(limits.max()+VSMALL), 1/(limits.min()+VSMALL));
+
         Info<< "Damped flow time scale min/max = "
-            << gMin(1/rDeltaT.primitiveField())
-            << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
+            << limits.min() << ", " << limits.max() << endl;
     }
 }
diff --git a/applications/solvers/multiphase/cavitatingFoam/setDeltaT.H b/applications/solvers/multiphase/cavitatingFoam/setDeltaT.H
index e1e7427fe0..1d2dd28f6b 100644
--- a/applications/solvers/multiphase/cavitatingFoam/setDeltaT.H
+++ b/applications/solvers/multiphase/cavitatingFoam/setDeltaT.H
@@ -36,13 +36,14 @@ Description
 if (adjustTimeStep)
 {
     scalar maxDeltaTFact =
-        min(maxCo/(CoNum + SMALL), maxAcousticCo/(acousticCoNum + SMALL));
+        Foam::min(maxCo/(CoNum + SMALL), maxAcousticCo/(acousticCoNum + SMALL));
 
-    scalar deltaTFact = min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
+    const scalar deltaTFact =
+        Foam::min(Foam::min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
 
     runTime.setDeltaT
     (
-        min
+        Foam::min
         (
             deltaTFact*runTime.deltaTValue(),
             maxDeltaT
diff --git a/applications/solvers/multiphase/cavitatingFoam/setInitialDeltaT.H b/applications/solvers/multiphase/cavitatingFoam/setInitialDeltaT.H
index 69fbe57333..435b2e4b0b 100644
--- a/applications/solvers/multiphase/cavitatingFoam/setInitialDeltaT.H
+++ b/applications/solvers/multiphase/cavitatingFoam/setInitialDeltaT.H
@@ -37,11 +37,15 @@ if (adjustTimeStep)
     if (CoNum > SMALL)
     {
         scalar maxDeltaTFact =
-            min(maxCo/(CoNum + SMALL), maxAcousticCo/(acousticCoNum + SMALL));
+            Foam::min
+            (
+                maxCo/(CoNum + SMALL),
+                maxAcousticCo/(acousticCoNum + SMALL)
+            );
 
         runTime.setDeltaT
         (
-            min
+            Foam::min
             (
                 maxDeltaTFact*runTime.deltaTValue(),
                 maxDeltaT
diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/setDeltaT.H b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/setDeltaT.H
index 026fb506bb..9a67a63d13 100644
--- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/setDeltaT.H
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/setDeltaT.H
@@ -36,13 +36,13 @@ Description
 if (adjustTimeStep)
 {
     scalar maxDeltaTFact =
-        min
+        Foam::min
         (
             maxCo/(CoNum + SMALL),
-            min
+            Foam::min
             (
                 maxAlphaCo/(alphaCoNum + SMALL),
-                min
+                Foam::min
                 (
                     maxAlphaDdt/(ddtAlphaNum + SMALL),
                     maxDi/(DiNum + SMALL)
@@ -50,16 +50,18 @@ if (adjustTimeStep)
             )
         );
 
-    scalar deltaTFact = min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
+    const scalar deltaTFact =
+        Foam::min(Foam::min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
 
     runTime.setDeltaT
     (
-        min
+        Foam::min
         (
             deltaTFact*runTime.deltaTValue(),
             maxDeltaT
         )
     );
+
     Info<< "deltaT = " <<  runTime.deltaTValue() << endl;
 }
 
diff --git a/applications/solvers/multiphase/reactingMultiphaseEulerFoam/setRDeltaT.H b/applications/solvers/multiphase/reactingMultiphaseEulerFoam/setRDeltaT.H
index 1d0c65b153..c15dfe4a9f 100644
--- a/applications/solvers/multiphase/reactingMultiphaseEulerFoam/setRDeltaT.H
+++ b/applications/solvers/multiphase/reactingMultiphaseEulerFoam/setRDeltaT.H
@@ -38,7 +38,9 @@
 
     fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
 
+    auto limits = gMinMax(rDeltaT.primitiveField());
+    limits.reset(1/(limits.max()+VSMALL), 1/(limits.min()+VSMALL));
+
     Info<< "Flow time scale min/max = "
-        << gMin(1/rDeltaT.primitiveField())
-        << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
+        << limits.min() << ", " << limits.max() << endl;
 }
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/CourantNos.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/CourantNos.H
index 13c45ad51c..1bd1d9d036 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/CourantNos.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/CourantNos.H
@@ -8,5 +8,5 @@
 
     Info<< "Max Ur Courant Number = " << UrCoNum << endl;
 
-    CoNum = max(CoNum, UrCoNum);
+    CoNum = Foam::max(CoNum, UrCoNum);
 }
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/setRDeltaT.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/setRDeltaT.H
index 10a7c9e0f0..45d9e6c74a 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/setRDeltaT.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/setRDeltaT.H
@@ -31,7 +31,9 @@
 
     fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
 
+    auto limits = gMinMax(rDeltaT.primitiveField());
+    limits.reset(1/(limits.max()+VSMALL), 1/(limits.min()+VSMALL));
+
     Info<< "Flow time scale min/max = "
-        << gMin(1/rDeltaT.primitiveField())
-        << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
+        << limits.min() << ", " << limits.max() << endl;
 }
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/CourantNos.H b/applications/solvers/multiphase/twoPhaseEulerFoam/CourantNos.H
index d828f32769..96e1ece1ac 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/CourantNos.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/CourantNos.H
@@ -8,5 +8,5 @@
 
     Info<< "Max Ur Courant Number = " << UrCoNum << endl;
 
-    CoNum = max(CoNum, UrCoNum);
+    CoNum = Foam::max(CoNum, UrCoNum);
 }
diff --git a/applications/test/fieldMapping/Test-fieldMapping.C b/applications/test/fieldMapping/Test-fieldMapping.C
index bb4ad0576b..f2bd3a87c9 100644
--- a/applications/test/fieldMapping/Test-fieldMapping.C
+++ b/applications/test/fieldMapping/Test-fieldMapping.C
@@ -262,13 +262,12 @@ int main(int argc, char *argv[])
 
         // Check constant profile
         {
-            const scalar max = gMax(one);
-            const scalar min = gMin(one);
+            auto limits =  gMinMax(one);
 
-            Info<< "Uniform one field min = " << min
-                << "  max = " << max << endl;
+            Info<< "Uniform one field min = "
+                << limits.min() << "  max = " << limits.max() << endl;
 
-            if (isNotEqual(max, 1) || isNotEqual(min, 1))
+            if (isNotEqual(limits.min(), 1) || isNotEqual(limits.max(), 1))
             {
                 FatalErrorInFunction
                     << "Uniform volVectorField not preserved."
@@ -286,13 +285,12 @@ int main(int argc, char *argv[])
         {
             const scalarField diff = ccX-mesh.C().component(0);
 
-            const scalar max = gMax(diff);
-            const scalar min = gMin(diff);
+            auto limits =  gMinMax(diff);
 
-            Info<< "Linear profile field min = " << min
-                << "  max = " << max << endl;
+            Info<< "Linear profile field min = "
+                << limits.min() << "  max = " << limits.max() << endl;
 
-            if (isNotEqual(max, 0) || isNotEqual(min, 0))
+            if (isNotEqual(limits.min(), 0) || isNotEqual(limits.max(), 0))
             {
                 FatalErrorInFunction
                     << "Linear profile not preserved."
@@ -309,13 +307,12 @@ int main(int argc, char *argv[])
         // Check face field mapping
         if (surfaceOne.size())
         {
-            const scalar max = gMax(surfaceOne.primitiveField());
-            const scalar min = gMin(surfaceOne.primitiveField());
+            auto limits =  gMinMax(surfaceOne.primitiveField());
 
-            Info<< "Uniform surface field min = " << min
-                << "  max = " << max << endl;
+            Info<< "Uniform surface field min = "
+                << limits.min() << "  max = " << limits.max() << endl;
 
-            if (isNotEqual(max, 1) || isNotEqual(min, 1))
+            if (isNotEqual(limits.min(), 1) || isNotEqual(limits.max(), 1))
             {
                 FatalErrorInFunction
                     << "Uniform surfaceScalarField not preserved."
diff --git a/applications/test/hexRef8/Test-hexRef8.C b/applications/test/hexRef8/Test-hexRef8.C
index 12c6d5e38e..f122a6a92e 100644
--- a/applications/test/hexRef8/Test-hexRef8.C
+++ b/applications/test/hexRef8/Test-hexRef8.C
@@ -344,13 +344,12 @@ int main(int argc, char *argv[])
 
         // Check constant profile
         {
-            const scalar max = gMax(one);
-            const scalar min = gMin(one);
+            auto limits =  gMinMax(one);
 
-            Info<< "Uniform one field min = " << min
-                << "  max = " << max << endl;
+            Info<< "Uniform one field min = "
+                << limits.min() << "  max = " << limits.max() << endl;
 
-            if (isNotEqual(min, 1) || isNotEqual(max, 1))
+            if (isNotEqual(limits.min(), 1) || isNotEqual(limits.max(), 1))
             {
                 FatalErrorInFunction
                     << "Uniform volVectorField not preserved."
@@ -368,13 +367,12 @@ int main(int argc, char *argv[])
         {
             const scalarField diff = ccX-mesh.C().component(0);
 
-            const scalar max = gMax(diff);
-            const scalar min = gMin(diff);
+            auto limits =  gMinMax(diff);
 
-            Info<< "Linear profile field min = " << min
-                << "  max = " << max << endl;
+            Info<< "Linear profile field min = "
+                << limits.min() << "  max = " << limits.max() << endl;
 
-            if (isNotEqual(min, 0) || isNotEqual(max, 0))
+            if (isNotEqual(limits.min(), 0) || isNotEqual(limits.max(), 0))
             {
                 Info<< "Linear profile not preserved."
                     << " Min and max should both be 0.0. min:" << min
@@ -389,13 +387,12 @@ int main(int argc, char *argv[])
         // Check face field mapping
         if (surfaceOne.size())
         {
-            const scalar max = gMax(surfaceOne.primitiveField());
-            const scalar min = gMin(surfaceOne.primitiveField());
+            auto limits =  gMinMax(surfaceOne.primitiveField());
 
-            Info<< "Uniform surface field min = " << min
-                << "  max = " << max << endl;
+            Info<< "Uniform surface field min = "
+                << limits.min() << "  max = " << limits.max() << endl;
 
-            if (isNotEqual(min, 1) || isNotEqual(max, 1))
+            if (isNotEqual(limits.min(), 1) || isNotEqual(limits.max(), 1))
             {
                 FatalErrorInFunction
                     << "Uniform surfaceScalarField not preserved."
diff --git a/applications/utilities/mesh/advanced/refineHexMesh/refineHexMesh.C b/applications/utilities/mesh/advanced/refineHexMesh/refineHexMesh.C
index 0da12172db..fa387d0116 100644
--- a/applications/utilities/mesh/advanced/refineHexMesh/refineHexMesh.C
+++ b/applications/utilities/mesh/advanced/refineHexMesh/refineHexMesh.C
@@ -145,17 +145,22 @@ int main(int argc, char *argv[])
     hexRef8 meshCutter(mesh);
 
     // Some stats
-    Info<< "Read mesh:" << nl
-        << "    cells:" << mesh.globalData().nTotalCells() << nl
-        << "    faces:" << mesh.globalData().nTotalFaces() << nl
-        << "    points:" << mesh.globalData().nTotalPoints() << nl
-        << "    cellLevel :"
-        << " min:" << gMin(meshCutter.cellLevel())
-        << " max:" << gMax(meshCutter.cellLevel()) << nl
-        << "    pointLevel :"
-        << " min:" << gMin(meshCutter.pointLevel())
-        << " max:" << gMax(meshCutter.pointLevel()) << nl
-        << endl;
+    {
+        auto cellLimits = gMinMax(meshCutter.cellLevel());
+        auto pointLimits = gMinMax(meshCutter.pointLevel());
+
+        Info<< "Read mesh:" << nl
+            << "    cells:" << mesh.globalData().nTotalCells() << nl
+            << "    faces:" << mesh.globalData().nTotalFaces() << nl
+            << "    points:" << mesh.globalData().nTotalPoints() << nl
+            << "    cellLevel :"
+            << " min:" << cellLimits.min()
+            << " max:" << cellLimits.max() << nl
+            << "    pointLevel :"
+            << " min:" << pointLimits.min()
+            << " max:" << pointLimits.max() << nl
+            << endl;
+    }
 
 
     // Maintain 2:1 ratio
diff --git a/applications/utilities/mesh/advanced/selectCells/edgeStats.C b/applications/utilities/mesh/advanced/selectCells/edgeStats.C
index 579a989b4d..1f54ceb793 100644
--- a/applications/utilities/mesh/advanced/selectCells/edgeStats.C
+++ b/applications/utilities/mesh/advanced/selectCells/edgeStats.C
@@ -37,7 +37,6 @@ License
 const Foam::scalar Foam::edgeStats::edgeTol_ = 1e-3;
 
 
-
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 Foam::direction Foam::edgeStats::getNormalDir
@@ -45,26 +44,25 @@ Foam::direction Foam::edgeStats::getNormalDir
     const twoDPointCorrector* correct2DPtr
 ) const
 {
-    direction dir = 3;
-
     if (correct2DPtr)
     {
         const vector& normal = correct2DPtr->planeNormal();
 
-        if (mag(normal & vector(1, 0, 0)) > 1-edgeTol_)
+        if (mag(normal.x()) > 1-edgeTol_)
         {
-            dir = 0;
+            return vector::X;
         }
-        else if (mag(normal & vector(0, 1, 0)) > 1-edgeTol_)
+        else if (mag(normal.y()) > 1-edgeTol_)
         {
-            dir = 1;
+            return vector::Y;
         }
-        else if (mag(normal & vector(0, 0, 1)) > 1-edgeTol_)
+        else if (mag(normal.z()) > 1-edgeTol_)
         {
-            dir = 2;
+            return vector::Z;
         }
     }
-    return dir;
+
+    return direction(3);
 }
 
 
@@ -95,12 +93,9 @@ Foam::edgeStats::edgeStats(const polyMesh& mesh)
         {
             Info<< "Correcting for 2D motion" << endl << endl;
 
-            autoPtr<twoDPointCorrector> correct2DPtr
-            (
-                new twoDPointCorrector(mesh)
-            );
+            twoDPointCorrector correct2D(mesh);
 
-            normalDir_ = getNormalDir(&correct2DPtr());
+            normalDir_ = getNormalDir(&correct2D);
         }
     }
 }
@@ -122,24 +117,15 @@ Foam::edgeStats::edgeStats
 
 Foam::scalar Foam::edgeStats::minLen(Ostream& os) const
 {
-    label nX = 0;
-    label nY = 0;
-    label nZ = 0;
+    label nAny(0);
+    label nX(0);
+    label nY(0);
+    label nZ(0);
 
-    scalar minX = GREAT;
-    scalar maxX = -GREAT;
-    vector x(1, 0, 0);
-
-    scalar minY = GREAT;
-    scalar maxY = -GREAT;
-    vector y(0, 1, 0);
-
-    scalar minZ = GREAT;
-    scalar maxZ = -GREAT;
-    vector z(0, 0, 1);
-
-    scalar minOther = GREAT;
-    scalar maxOther = -GREAT;
+    scalarMinMax limitsAny(GREAT, -GREAT);
+    scalarMinMax limitsX(limitsAny);
+    scalarMinMax limitsY(limitsAny);
+    scalarMinMax limitsZ(limitsAny);
 
     const edgeList& edges = mesh_.edges();
 
@@ -151,58 +137,75 @@ Foam::scalar Foam::edgeStats::minLen(Ostream& os) const
 
         eVec /= eMag;
 
-        if (mag(eVec & x) > 1-edgeTol_)
+        if (mag(eVec.x()) > 1-edgeTol_)
         {
-            minX = min(minX, eMag);
-            maxX = max(maxX, eMag);
+            limitsX.add(eMag);
             nX++;
         }
-        else if (mag(eVec & y) > 1-edgeTol_)
+        else if (mag(eVec.y()) > 1-edgeTol_)
         {
-            minY = min(minY, eMag);
-            maxY = max(maxY, eMag);
+            limitsY.add(eMag);
             nY++;
         }
-        else if (mag(eVec & z) > 1-edgeTol_)
+        else if (mag(eVec.z()) > 1-edgeTol_)
         {
-            minZ = min(minZ, eMag);
-            maxZ = max(maxZ, eMag);
+            limitsZ.add(eMag);
             nZ++;
         }
         else
         {
-            minOther = min(minOther, eMag);
-            maxOther = max(maxOther, eMag);
+            limitsAny.add(eMag);
+            nAny++;
         }
     }
 
     os  << "Mesh bounding box:" << boundBox(mesh_.points()) << nl << nl
         << "Mesh edge statistics:" << nl
-        << "    x aligned :  number:" << nX << "\tminLen:" << minX
-        << "\tmaxLen:" << maxX << nl
-        << "    y aligned :  number:" << nY << "\tminLen:" << minY
-        << "\tmaxLen:" << maxY << nl
-        << "    z aligned :  number:" << nZ << "\tminLen:" << minZ
-        << "\tmaxLen:" << maxZ << nl
-        << "    other     :  number:" << mesh_.nEdges() - nX - nY - nZ
-        << "\tminLen:" << minOther
-        << "\tmaxLen:" << maxOther << nl << endl;
+        << "    x aligned :  number:" << nX
+        << "\tminLen:" << limitsX.min() << "\tmaxLen:" << limitsX.max() << nl
+        << "    y aligned :  number:" << nY
+        << "\tminLen:" << limitsY.min() << "\tmaxLen:" << limitsY.max() << nl
+        << "    z aligned :  number:" << nZ
+        << "\tminLen:" << limitsZ.min() << "\tmaxLen:" << limitsZ.max() << nl
+        << "    other     :  number:" << nAny
+        << "\tminLen:" << limitsAny.min()
+        << "\tmaxLen:" << limitsAny.max() << nl << endl;
 
-    if (normalDir_ == 0)
+    if (normalDir_ == vector::X)
     {
-        return min(minY, min(minZ, minOther));
+        return Foam::min
+        (
+            limitsAny.min(),
+            Foam::min(limitsY.min(), limitsZ.min())
+        );
     }
-    else if (normalDir_ == 1)
+    else if (normalDir_ == vector::Y)
     {
-        return min(minX, min(minZ, minOther));
+        return Foam::min
+        (
+            limitsAny.min(),
+            Foam::min(limitsX.min(), limitsZ.min())
+        );
     }
-    else if (normalDir_ == 2)
+    else if (normalDir_ == vector::Z)
     {
-        return min(minX, min(minY, minOther));
+        return Foam::min
+        (
+            limitsAny.min(),
+            Foam::min(limitsX.min(), limitsY.min())
+        );
     }
     else
     {
-        return min(minX, min(minY, min(minZ, minOther)));
+        return Foam::min
+        (
+            limitsAny.min(),
+            Foam::min
+            (
+                limitsX.min(),
+                Foam::min(limitsY.min(), limitsZ.min())
+            )
+        );
     }
 }
 
diff --git a/applications/utilities/mesh/advanced/snappyRefineMesh/snappyRefineMesh.C b/applications/utilities/mesh/advanced/snappyRefineMesh/snappyRefineMesh.C
index 6aa76e6653..2b537081fe 100644
--- a/applications/utilities/mesh/advanced/snappyRefineMesh/snappyRefineMesh.C
+++ b/applications/utilities/mesh/advanced/snappyRefineMesh/snappyRefineMesh.C
@@ -81,26 +81,25 @@ void writeSet(const cellSet& cells, const string& msg)
 
 direction getNormalDir(const twoDPointCorrector* correct2DPtr)
 {
-    direction dir = 3;
-
     if (correct2DPtr)
     {
         const vector& normal = correct2DPtr->planeNormal();
 
-        if (mag(normal & vector(1, 0, 0)) > 1-edgeTol)
+        if (mag(normal.x()) > 1-edgeTol)
         {
-            dir = 0;
+            return vector::X;
         }
-        else if (mag(normal & vector(0, 1, 0)) > 1-edgeTol)
+        else if (mag(normal.y()) > 1-edgeTol)
         {
-            dir = 1;
+            return vector::Y;
         }
-        else if (mag(normal & vector(0, 0, 1)) > 1-edgeTol)
+        else if (mag(normal.z()) > 1-edgeTol)
         {
-            dir = 2;
+            return vector::Z;
         }
     }
-    return dir;
+
+    return direction(3);
 }
 
 
@@ -109,89 +108,95 @@ direction getNormalDir(const twoDPointCorrector* correct2DPtr)
 // directions but exclude component (0=x, 1=y, 2=z, other=none)
 scalar getEdgeStats(const primitiveMesh& mesh, const direction excludeCmpt)
 {
-    label nX = 0;
-    label nY = 0;
-    label nZ = 0;
+    label nAny(0);
+    label nX(0);
+    label nY(0);
+    label nZ(0);
 
-    scalar minX = GREAT;
-    scalar maxX = -GREAT;
-    vector x(1, 0, 0);
-
-    scalar minY = GREAT;
-    scalar maxY = -GREAT;
-    vector y(0, 1, 0);
-
-    scalar minZ = GREAT;
-    scalar maxZ = -GREAT;
-    vector z(0, 0, 1);
-
-    scalar minOther = GREAT;
-    scalar maxOther = -GREAT;
+    scalarMinMax limitsAny(GREAT, -GREAT);
+    scalarMinMax limitsX(limitsAny);
+    scalarMinMax limitsY(limitsAny);
+    scalarMinMax limitsZ(limitsAny);
 
     const edgeList& edges = mesh.edges();
 
-    forAll(edges, edgei)
+    for (const edge& e : edges)
     {
-        const edge& e = edges[edgei];
-
         vector eVec(e.vec(mesh.points()));
 
         scalar eMag = mag(eVec);
 
         eVec /= eMag;
 
-        if (mag(eVec & x) > 1-edgeTol)
+        if (mag(eVec.x()) > 1-edgeTol)
         {
-            minX = min(minX, eMag);
-            maxX = max(maxX, eMag);
+            limitsX.add(eMag);
             nX++;
         }
-        else if (mag(eVec & y) > 1-edgeTol)
+        else if (mag(eVec.y()) > 1-edgeTol)
         {
-            minY = min(minY, eMag);
-            maxY = max(maxY, eMag);
+            limitsY.add(eMag);
             nY++;
         }
-        else if (mag(eVec & z) > 1-edgeTol)
+        else if (mag(eVec.z()) > 1-edgeTol)
         {
-            minZ = min(minZ, eMag);
-            maxZ = max(maxZ, eMag);
+            limitsZ.add(eMag);
             nZ++;
         }
         else
         {
-            minOther = min(minOther, eMag);
-            maxOther = max(maxOther, eMag);
+            limitsAny.add(eMag);
+            nAny++;
         }
     }
 
     Info<< "Mesh bounding box:" << boundBox(mesh.points()) << nl << nl
         << "Mesh edge statistics:" << nl
-        << "    x aligned :  number:" << nX << "\tminLen:" << minX
-        << "\tmaxLen:" << maxX << nl
-        << "    y aligned :  number:" << nY << "\tminLen:" << minY
-        << "\tmaxLen:" << maxY << nl
-        << "    z aligned :  number:" << nZ << "\tminLen:" << minZ
-        << "\tmaxLen:" << maxZ << nl
-        << "    other     :  number:" << mesh.nEdges() - nX - nY - nZ
-        << "\tminLen:" << minOther
-        << "\tmaxLen:" << maxOther << nl << endl;
+        << "    x aligned :  number:" << nX
+        << "\tminLen:" << limitsX.min() << "\tmaxLen:" << limitsX.max() << nl
+        << "    y aligned :  number:" << nY
+        << "\tminLen:" << limitsY.min() << "\tmaxLen:" << limitsY.max() << nl
+        << "    z aligned :  number:" << nZ
+        << "\tminLen:" << limitsZ.min() << "\tmaxLen:" << limitsZ.max() << nl
+        << "    other     :  number:" << nAny
+        << "\tminLen:" << limitsAny.min()
+        << "\tmaxLen:" << limitsAny.max() << nl << endl;
 
-    if (excludeCmpt == 0)
+    if (excludeCmpt == vector::X)
     {
-        return min(minY, min(minZ, minOther));
+        return Foam::min
+        (
+            limitsAny.min(),
+            Foam::min(limitsY.min(), limitsZ.min())
+        );
     }
-    else if (excludeCmpt == 1)
+    else if (excludeCmpt == vector::Y)
     {
-        return min(minX, min(minZ, minOther));
+        return Foam::min
+        (
+            limitsAny.min(),
+            Foam::min(limitsX.min(), limitsZ.min())
+        );
     }
-    else if (excludeCmpt == 2)
+    else if (excludeCmpt == vector::Z)
     {
-        return min(minX, min(minY, minOther));
+        return Foam::min
+        (
+            limitsAny.min(),
+            Foam::min(limitsX.min(), limitsY.min())
+        );
     }
     else
     {
-        return min(minX, min(minY, min(minZ, minOther)));
+        return Foam::min
+        (
+            limitsAny.min(),
+            Foam::min
+            (
+                limitsX.min(),
+                Foam::min(limitsY.min(), limitsZ.min())
+            )
+        );
     }
 }
 
diff --git a/applications/utilities/mesh/manipulation/refineMesh/refineMesh.C b/applications/utilities/mesh/manipulation/refineMesh/refineMesh.C
index 0076597534..3a51efb298 100644
--- a/applications/utilities/mesh/manipulation/refineMesh/refineMesh.C
+++ b/applications/utilities/mesh/manipulation/refineMesh/refineMesh.C
@@ -63,96 +63,72 @@ static const scalar edgeTol = 1e-3;
 // Print edge statistics on mesh.
 void printEdgeStats(const polyMesh& mesh)
 {
-    label nX = 0;
-    label nY = 0;
-    label nZ = 0;
+    label nAny(0);
+    label nX(0);
+    label nY(0);
+    label nZ(0);
 
-    scalar minX = GREAT;
-    scalar maxX = -GREAT;
-    static const vector x(1, 0, 0);
-
-    scalar minY = GREAT;
-    scalar maxY = -GREAT;
-    static const vector y(0, 1, 0);
-
-    scalar minZ = GREAT;
-    scalar maxZ = -GREAT;
-    static const vector z(0, 0, 1);
-
-    scalar minOther = GREAT;
-    scalar maxOther = -GREAT;
+    scalarMinMax limitsAny(GREAT, -GREAT);
+    scalarMinMax limitsX(limitsAny);
+    scalarMinMax limitsY(limitsAny);
+    scalarMinMax limitsZ(limitsAny);
 
     bitSet isMasterEdge(syncTools::getMasterEdges(mesh));
 
     const edgeList& edges = mesh.edges();
 
-    forAll(edges, edgeI)
+    for (const label edgei : isMasterEdge)
     {
-        if (isMasterEdge.test(edgeI))
+        const edge& e = edges[edgei];
+
+        vector eVec(e.vec(mesh.points()));
+
+        scalar eMag = mag(eVec);
+
+        eVec /= eMag;
+
+        if (mag(eVec.x()) > 1-edgeTol)
         {
-            const edge& e = edges[edgeI];
-
-            vector eVec(e.vec(mesh.points()));
-
-            scalar eMag = mag(eVec);
-
-            eVec /= eMag;
-
-            if (mag(eVec & x) > 1-edgeTol)
-            {
-                minX = min(minX, eMag);
-                maxX = max(maxX, eMag);
-                nX++;
-            }
-            else if (mag(eVec & y) > 1-edgeTol)
-            {
-                minY = min(minY, eMag);
-                maxY = max(maxY, eMag);
-                nY++;
-            }
-            else if (mag(eVec & z) > 1-edgeTol)
-            {
-                minZ = min(minZ, eMag);
-                maxZ = max(maxZ, eMag);
-                nZ++;
-            }
-            else
-            {
-                minOther = min(minOther, eMag);
-                maxOther = max(maxOther, eMag);
-            }
+            limitsX.add(eMag);
+            nX++;
+        }
+        else if (mag(eVec.y()) > 1-edgeTol)
+        {
+            limitsY.add(eMag);
+            nY++;
+        }
+        else if (mag(eVec.z()) > 1-edgeTol)
+        {
+            limitsZ.add(eMag);
+            nZ++;
+        }
+        else
+        {
+            limitsAny.add(eMag);
+            nAny++;
         }
     }
 
-    label nEdges = mesh.nEdges();
-    reduce(nEdges, sumOp<label>());
     reduce(nX, sumOp<label>());
     reduce(nY, sumOp<label>());
     reduce(nZ, sumOp<label>());
+    reduce(nAny, sumOp<label>());
 
-    reduce(minX, minOp<scalar>());
-    reduce(maxX, maxOp<scalar>());
-
-    reduce(minY, minOp<scalar>());
-    reduce(maxY, maxOp<scalar>());
-
-    reduce(minZ, minOp<scalar>());
-    reduce(maxZ, maxOp<scalar>());
-
-    reduce(minOther, minOp<scalar>());
-    reduce(maxOther, maxOp<scalar>());
-
+    reduce(limitsX, sumOp<scalarMinMax>());
+    reduce(limitsY, sumOp<scalarMinMax>());
+    reduce(limitsZ, sumOp<scalarMinMax>());
+    reduce(limitsAny, sumOp<scalarMinMax>());
 
     Info<< "Mesh edge statistics:" << nl
-        << "    x aligned :  number:" << nX << "\tminLen:" << minX
-        << "\tmaxLen:" << maxX << nl
-        << "    y aligned :  number:" << nY << "\tminLen:" << minY
-        << "\tmaxLen:" << maxY << nl
-        << "    z aligned :  number:" << nZ << "\tminLen:" << minZ
-        << "\tmaxLen:" << maxZ << nl
-        << "    other     :  number:" << nEdges - nX - nY - nZ
-        << "\tminLen:" << minOther
-        << "\tmaxLen:" << maxOther << nl << endl;
+        << "    x aligned :  number:" << nX
+        << "\tminLen:" << limitsX.min() << "\tmaxLen:" << limitsX.max() << nl
+        << "    y aligned :  number:" << nY
+        << "\tminLen:" << limitsY.min() << "\tmaxLen:" << limitsY.max() << nl
+        << "    z aligned :  number:" << nZ
+        << "\tminLen:" << limitsZ.min() << "\tmaxLen:" << limitsZ.max() << nl
+        << "    other     :  number:" << nAny
+        << "\tminLen:" << limitsAny.min()
+        << "\tmaxLen:" << limitsAny.max() << nl << endl;
 }
 
 
diff --git a/applications/utilities/preProcessing/createBoxTurb/createBoxTurb.C b/applications/utilities/preProcessing/createBoxTurb/createBoxTurb.C
index 5e01b77ada..4c1400fc08 100644
--- a/applications/utilities/preProcessing/createBoxTurb/createBoxTurb.C
+++ b/applications/utilities/preProcessing/createBoxTurb/createBoxTurb.C
@@ -167,18 +167,24 @@ int main(int argc, char *argv[])
         Info<< "Generating kinetic energy field" << endl;
         volScalarField k("k", 0.5*magSqr(U));
         k.write();
+
+        auto limits = gMinMax(k);
+        auto avg = gAverage(k);
+
         Info<< "min/max/average k = "
-            << gMin(k) << ", " << gMax(k) << ", " << gAverage(k)
-            << endl;
+            << limits.min() << ", " << limits.max() << ", " << avg << endl;
     }
 
     {
         Info<< "Generating div(U) field" << endl;
         volScalarField divU(fvc::div(U));
         divU.write();
+
+        auto limits = gMinMax(divU);
+        auto avg = gAverage(divU);
+
         Info<< "min/max/average div(U) = "
-            << gMin(divU) << ", " << gMax(divU) << ", " << gAverage(divU)
-            << endl;
+            << limits.min() << ", " << limits.max() << ", " << avg << endl;
     }
 
     Info<< nl;
diff --git a/applications/utilities/preProcessing/setAlphaField/setAlphaField.C b/applications/utilities/preProcessing/setAlphaField/setAlphaField.C
index 89e10cd444..867e8762b7 100644
--- a/applications/utilities/preProcessing/setAlphaField/setAlphaField.C
+++ b/applications/utilities/preProcessing/setAlphaField/setAlphaField.C
@@ -223,11 +223,16 @@ int main(int argc, char *argv[])
     Info<< "Writing new alpha field " << alpha1.name() << endl;
     alpha1.write();
 
-    const scalarField& alpha = alpha1.internalField();
+    {
+        const auto& alpha = alpha1.primitiveField();
 
-    Info<< "sum(alpha*V):" << gSum(mesh.V()*alpha)
-        << ", 1-max(alpha1): " << 1 - gMax(alpha)
-        << " min(alpha1): " << gMin(alpha) << endl;
+        auto limits = gMinMax(alpha);
+        auto total = gWeightedSum(mesh.V(), alpha);
+
+        Info<< "sum(alpha*V):" << total
+            << ", 1-max(alpha1): " << 1 - limits.max()
+            << " min(alpha1): " << limits.min() << endl;
+    }
 
     Info<< "\nEnd\n" << endl;
 
diff --git a/src/combustionModels/EDC/EDC.C b/src/combustionModels/EDC/EDC.C
index 90fc7e2605..977270aa9d 100644
--- a/src/combustionModels/EDC/EDC.C
+++ b/src/combustionModels/EDC/EDC.C
@@ -180,8 +180,9 @@ void Foam::combustionModels::EDC<ReactionThermo>::correct()
             kappa_.correctBoundaryConditions();
         }
 
-        Info<< "Chemistry time solved max/min : "
-            << gMax(tauStar) << " / " << gMin(tauStar) << endl;
+        auto limits = gMinMax(tauStar);
+        Info<< "Chemistry time solved min/max : "
+            << limits.min() << ", " << limits.max() << endl;
 
         this->chemistryPtr_->solve(tauStar);
     }
diff --git a/src/dynamicFaMesh/interfaceTrackingFvMesh/interfaceTrackingFvMesh.C b/src/dynamicFaMesh/interfaceTrackingFvMesh/interfaceTrackingFvMesh.C
index aa5bf5a4fd..d74fe5b006 100644
--- a/src/dynamicFaMesh/interfaceTrackingFvMesh/interfaceTrackingFvMesh.C
+++ b/src/dynamicFaMesh/interfaceTrackingFvMesh/interfaceTrackingFvMesh.C
@@ -2022,8 +2022,11 @@ bool Foam::interfaceTrackingFvMesh::update()
 
         const scalarField& K = aMesh().faceCurvatures().internalField();
 
-        Info<< "Free surface curvature: min = " << gMin(K)
-            << ", max = " << gMax(K) << ", average = " << gAverage(K) << nl;
+        auto limits = gMinMax(K);
+
+        Info<< "Free surface curvature: min = " << limits.min()
+            << ", max = " << limits.max()
+            << ", average = " << gAverage(K) << nl;
 
         timeIndex_ = mesh().time().timeIndex();
     }
diff --git a/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C b/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C
index 93d9c97a8a..7f209b9ef6 100644
--- a/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C
+++ b/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C
@@ -879,10 +879,13 @@ bool Foam::motionSmootherAlgo::scaleMesh
         vector::zero    // null value
     );
 
-    Info<< "Moving mesh using displacement scaling :"
-        << " min:" << gMin(scale_.primitiveField())
-        << "  max:" << gMax(scale_.primitiveField())
-        << endl;
+    {
+        auto limits = gMinMax(scale_.primitiveField());
+        Info<< "Moving mesh using displacement scaling :"
+            << " min:" << limits.min()
+            << " max:" << limits.max()
+            << endl;
+    }
 
     // Get points using current displacement and scale. Optionally 2D corrected.
     pointField newPoints(curPoints());
@@ -1018,9 +1021,11 @@ bool Foam::motionSmootherAlgo::scaleMesh
 
         if (debug)
         {
+            auto limits = gMinMax(scale_);
+
             Pout<< "scale_ after smoothing :"
-                << " min:" << Foam::gMin(scale_)
-                << " max:" << Foam::gMax(scale_)
+                << " min:" << limits.min()
+                << " max:" << limits.max()
                 << endl;
         }
 
diff --git a/src/dynamicMesh/motionSolvers/componentDisplacement/componentDisplacementMotionSolver.C b/src/dynamicMesh/motionSolvers/componentDisplacement/componentDisplacementMotionSolver.C
index b0a98da199..01b0642f87 100644
--- a/src/dynamicMesh/motionSolvers/componentDisplacement/componentDisplacementMotionSolver.C
+++ b/src/dynamicMesh/motionSolvers/componentDisplacement/componentDisplacementMotionSolver.C
@@ -166,9 +166,7 @@ void Foam::componentDisplacementMotionSolver::updateMesh(const mapPolyMesh& mpm)
     );
 
     // Get extents of points0 and points and determine scale
-    const scalar scale =
-        (gMax(points0_)-gMin(points0_))
-       /(gMax(points)-gMin(points));
+    const scalar scale = gMinMax(points0_).span() / gMinMax(points).span();
 
     scalarField newPoints0(mpm.pointMap().size());
 
diff --git a/src/dynamicMesh/motionSolvers/displacement/interpolation/displacementInterpolationMotionSolver.C b/src/dynamicMesh/motionSolvers/displacement/interpolation/displacementInterpolationMotionSolver.C
index 2026e4e09e..24b07f1718 100644
--- a/src/dynamicMesh/motionSolvers/displacement/interpolation/displacementInterpolationMotionSolver.C
+++ b/src/dynamicMesh/motionSolvers/displacement/interpolation/displacementInterpolationMotionSolver.C
@@ -138,16 +138,18 @@ void Foam::displacementInterpolationMotionSolver::calcInterpolation()
             scalar minCoord = VGREAT;
             scalar maxCoord = -VGREAT;
 
-            forAll(fz().meshPoints(), localI)
+            scalarMinMax limits;
+
+            for (const label pointi : fz().meshPoints())
             {
-                label pointi = fz().meshPoints()[localI];
                 const scalar coord = points0()[pointi][dir];
-                minCoord = min(minCoord, coord);
-                maxCoord = max(maxCoord, coord);
+                limits.add(coord);
             }
 
-            zoneCoordinates[2*i] = returnReduce(minCoord, minOp<scalar>());
-            zoneCoordinates[2*i+1] = returnReduce(maxCoord, maxOp<scalar>());
+            reduce(limits, sumOp<scalarMinMax>());
+
+            zoneCoordinates[2*i] = limits.min();
+            zoneCoordinates[2*i+1] = limits.max();
 
             if (debug)
             {
@@ -167,8 +169,13 @@ void Foam::displacementInterpolationMotionSolver::calcInterpolation()
         // Check if we have static min and max mesh bounds
         const scalarField meshCoords(points0().component(dir));
 
-        scalar minCoord = gMin(meshCoords);
-        scalar maxCoord = gMax(meshCoords);
+        scalar minCoord, maxCoord;
+        {
+            auto limits = gMinMax(meshCoords);
+
+            minCoord = limits.min();
+            maxCoord = limits.max();
+        }
 
         if (debug)
         {
diff --git a/src/dynamicMesh/motionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.C b/src/dynamicMesh/motionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.C
index 645be3cea2..14b8ed41c7 100644
--- a/src/dynamicMesh/motionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.C
+++ b/src/dynamicMesh/motionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.C
@@ -386,15 +386,21 @@ void Foam::displacementLayeredMotionMotionSolver::cellZoneSolve
             patchi
         );
 
-        DebugInfo
-            << "For cellZone:" << cellZoneI
-            << " for faceZone:" << fz.name()
-            << " nPoints:" << tseed().size()
-            << " have patchField:"
-            << " max:" << gMax(tseed())
-            << " min:" << gMin(tseed())
-            << " avg:" << gAverage(tseed())
-            << endl;
+        if (debug)
+        {
+            auto limits = gMinMax(tseed());
+            auto avg = gAverage(tseed());
+
+            Info
+                << "For cellZone:" << cellZoneI
+                << " for faceZone:" << fz.name()
+                << " nPoints:" << tseed().size()
+                << " have patchField:"
+                << " min:" << limits.min()
+                << " max:" << limits.max()
+                << " avg:" << avg
+                << endl;
+        }
 
 
         // Set distance and transported value
diff --git a/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/edgeLimitedFaGrad/edgeLimitedFaGrads.C b/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/edgeLimitedFaGrad/edgeLimitedFaGrads.C
index a68d4f63ec..2e11e21291 100644
--- a/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/edgeLimitedFaGrad/edgeLimitedFaGrads.C
+++ b/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/edgeLimitedFaGrad/edgeLimitedFaGrads.C
@@ -175,10 +175,13 @@ tmp<areaVectorField> edgeLimitedGrad<scalar>::calcGrad
 
     if (fa::debug)
     {
+        auto limits = gMinMax(limiter);
+        auto avg = gAverage(limiter);
+
         Info<< "gradient limiter for: " << vsf.name()
-            << " max = " << gMax(limiter)
-            << " min = " << gMin(limiter)
-            << " average: " << gAverage(limiter) << endl;
+            << " min = " << limits.min()
+            << " max = " << limits.max()
+            << " average: " << avg << endl;
     }
 
     g.primitiveFieldRef() *= limiter;
@@ -312,10 +315,13 @@ tmp<areaTensorField> edgeLimitedGrad<vector>::calcGrad
 
     if (fa::debug)
     {
+        auto limits = gMinMax(limiter);
+        auto avg = gAverage(limiter);
+
         Info<< "gradient limiter for: " << vvf.name()
-            << " max = " << gMax(limiter)
-            << " min = " << gMin(limiter)
-            << " average: " << gAverage(limiter) << endl;
+            << " min = " << limits.min()
+            << " max = " << limits.max()
+            << " average: " << avg << endl;
     }
 
     g.primitiveFieldRef() *= limiter;
diff --git a/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/faceLimitedFaGrad/faceLimitedFaGrads.C b/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/faceLimitedFaGrad/faceLimitedFaGrads.C
index eb0cc9bcbd..80a2a60a6a 100644
--- a/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/faceLimitedFaGrad/faceLimitedFaGrads.C
+++ b/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/faceLimitedFaGrad/faceLimitedFaGrads.C
@@ -221,10 +221,13 @@ tmp<areaVectorField> faceLimitedGrad<scalar>::calcGrad
 
     if (fa::debug)
     {
+        auto limits = gMinMax(limiter);
+        auto avg = gAverage(limiter);
+
         Info<< "gradient limiter for: " << vsf.name()
-            << " max = " << gMax(limiter)
-            << " min = " << gMin(limiter)
-            << " average: " << gAverage(limiter) << endl;
+            << " min = " << limits.min()
+            << " max = " << limits.max()
+            << " average: " << avg << endl;
     }
 
     g.primitiveFieldRef() *= limiter;
@@ -370,10 +373,13 @@ tmp<areaTensorField> faceLimitedGrad<vector>::calcGrad
 
     if (fa::debug)
     {
+        auto limits = gMinMax(limiter);
+        auto avg = gAverage(limiter);
+
         Info<< "gradient limiter for: " << vsf.name()
-            << " max = " << gMax(limiter)
-            << " min = " << gMin(limiter)
-            << " average: " << gAverage(limiter) << endl;
+            << " min = " << limits.min()
+            << " max = " << limits.max()
+            << " average: " << avg << endl;
     }
 
     tensorField& gIf = g.primitiveFieldRef();
diff --git a/src/finiteVolume/cfdTools/general/include/setDeltaT.H b/src/finiteVolume/cfdTools/general/include/setDeltaT.H
index 25ee3713dd..aa80893203 100644
--- a/src/finiteVolume/cfdTools/general/include/setDeltaT.H
+++ b/src/finiteVolume/cfdTools/general/include/setDeltaT.H
@@ -36,11 +36,13 @@ Description
 if (adjustTimeStep)
 {
     scalar maxDeltaTFact = maxCo/(CoNum + SMALL);
-    scalar deltaTFact = min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
+
+    const scalar deltaTFact =
+        Foam::min(Foam::min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
 
     runTime.setDeltaT
     (
-        min
+        Foam::min
         (
             deltaTFact*runTime.deltaTValue(),
             maxDeltaT
diff --git a/src/finiteVolume/cfdTools/general/include/setInitialDeltaT.H b/src/finiteVolume/cfdTools/general/include/setInitialDeltaT.H
index 55f23f13e1..1f45166220 100644
--- a/src/finiteVolume/cfdTools/general/include/setInitialDeltaT.H
+++ b/src/finiteVolume/cfdTools/general/include/setInitialDeltaT.H
@@ -38,10 +38,10 @@ if (adjustTimeStep)
     {
         runTime.setDeltaT
         (
-            min
+            Foam::min
             (
                 maxCo*runTime.deltaTValue()/CoNum,
-                min(runTime.deltaTValue(), maxDeltaT)
+                Foam::min(runTime.deltaTValue(), maxDeltaT)
             )
         );
     }
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/electrostaticDeposition/electrostaticDepositionFvPatchScalarField.C b/src/finiteVolume/fields/fvPatchFields/derived/electrostaticDeposition/electrostaticDepositionFvPatchScalarField.C
index 2ffd8892c1..61e2959eac 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/electrostaticDeposition/electrostaticDepositionFvPatchScalarField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/electrostaticDeposition/electrostaticDepositionFvPatchScalarField.C
@@ -496,16 +496,15 @@ void Foam::electrostaticDepositionFvPatchScalarField::updateCoeffs()
     timei_ = db().time().timeIndex();
 
     {
-        const scalar hMin = gMin(h_);
-        const scalar hMax = gMax(h_);
-        const scalar hAvg = gAverage(h_);
+        auto limits = gMinMax(h_);
+        auto avg = gAverage(h_);
 
-        if (Pstream::master())
+        if (UPstream::master())
         {
             Info<< "    patch: " << patch().name()
-                << ", h: min = " << hMin
-                << ", max = " << hMax
-                << ", average = " << hAvg << nl
+                << ", h: min = " << limits.min()
+                << ", max = " << limits.max()
+                << ", average = " << avg << nl
                 << endl;
         }
     }
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/mappedField/mappedFieldFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/mappedField/mappedFieldFvPatchField.C
index cd3da65504..3ad7284506 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/mappedField/mappedFieldFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/mappedField/mappedFieldFvPatchField.C
@@ -175,11 +175,14 @@ void Foam::mappedFieldFvPatchField<Type>::updateCoeffs()
 
     if (debug)
     {
+        auto limits = gMinMax(*this);
+        auto avg = gAverage(*this);
+
         Info<< "operating on field:" << this->internalField().name()
             << " patch:" << this->patch().name()
-            << "  avg:" << gAverage(*this)
-            << "  min:" << gMin(*this)
-            << "  max:" << gMax(*this)
+            << "  avg:" << avg
+            << "  min:" << limits.min()
+            << "  max:" << limits.max()
             << endl;
     }
 
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/mappedField/mappedMixedFieldFvPatchField/mappedMixedFieldFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/mappedField/mappedMixedFieldFvPatchField/mappedMixedFieldFvPatchField.C
index d0a225f0e4..643b495772 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/mappedField/mappedMixedFieldFvPatchField/mappedMixedFieldFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/mappedField/mappedMixedFieldFvPatchField/mappedMixedFieldFvPatchField.C
@@ -173,6 +173,9 @@ void Foam::mappedMixedFieldFvPatchField<Type>::updateCoeffs()
 
     if (debug)
     {
+        auto limits = gMinMax(*this);
+        auto avg = gAverage(*this);
+
         Info<< this->patch().boundaryMesh().mesh().name() << ':'
             << this->patch().name() << ':'
             << this->internalField().name() << " <- "
@@ -180,9 +183,9 @@ void Foam::mappedMixedFieldFvPatchField<Type>::updateCoeffs()
             << this->mapper_.samplePatch() << ':'
             << this->fieldName_ << " :"
             << " value "
-            << " min:" << gMin(*this)
-            << " max:" << gMax(*this)
-            << " avg:" << gAverage(*this)
+            << " min:" << limits.min()
+            << " max:" << limits.max()
+            << " avg:" << avg
             << endl;
     }
 }
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/mappedFixedValue/mappedFixedValueFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/mappedFixedValue/mappedFixedValueFvPatchField.C
index ec190df0f1..58970edd6e 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/mappedFixedValue/mappedFixedValueFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/mappedFixedValue/mappedFixedValueFvPatchField.C
@@ -122,12 +122,15 @@ void Foam::mappedFixedValueFvPatchField<Type>::updateCoeffs()
 
     if (debug)
     {
+        auto limits = gMinMax(*this);
+        auto avg = gAverage(*this);
+
         Info<< "mapped on field:"
             << this->internalField().name()
             << " patch:" << this->patch().name()
-            << "  avg:" << gAverage(*this)
-            << "  min:" << gMin(*this)
-            << "  max:" << gMax(*this)
+            << "  avg:" << avg
+            << "  min:" << limits.min()
+            << "  max:" << limits.max()
             << endl;
     }
 
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/mappedMixed/mappedMixedFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/mappedMixed/mappedMixedFvPatchField.C
index f3d9870e0d..95d553d508 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/mappedMixed/mappedMixedFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/mappedMixed/mappedMixedFvPatchField.C
@@ -215,6 +215,9 @@ void Foam::mappedMixedFvPatchField<Type>::updateCoeffs()
 
     if (debug)
     {
+        auto limits = gMinMax(*this);
+        auto avg = gAverage(*this);
+
         Info<< this->patch().boundaryMesh().mesh().name() << ':'
             << this->patch().name() << ':'
             << this->internalField().name() << " <- "
@@ -222,9 +225,9 @@ void Foam::mappedMixedFvPatchField<Type>::updateCoeffs()
             << this->mapper_.samplePatch() << ':'
             << this->fieldName_ << " :"
             << " value "
-            << " min:" << gMin(*this)
-            << " max:" << gMax(*this)
-            << " avg:" << gAverage(*this)
+            << " min:" << limits.min()
+            << " max:" << limits.max()
+            << " avg:" << avg
             << endl;
     }
 }
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C
index 40dc0780cc..f5a3a61818 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C
@@ -181,9 +181,12 @@ void Foam::timeVaryingMappedFixedValueFvPatchField<Type>::updateCoeffs()
 
     if (debug)
     {
-        Pout<< "updateCoeffs : set fixedValue to min:" << gMin(*this)
-            << " max:" << gMax(*this)
-            << " avg:" << gAverage(*this) << endl;
+        auto limits = gMinMax(*this);
+        auto avg = gAverage(*this);
+
+        Pout<< "updateCoeffs : set fixedValue to min:" << limits.min()
+            << " max:" << limits.max()
+            << " avg:" << avg << endl;
     }
 
     fixedValueFvPatchField<Type>::updateCoeffs();
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C
index 888e5bc6c0..24114e2055 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C
@@ -962,6 +962,8 @@ void Foam::turbulentDFSEMInletFvPatchVectorField::updateCoeffs()
 
         if (debug)
         {
+            auto limits = gMinMax(*this);
+
             Info<< "Magnitude of bulk velocity: " << UBulk << endl;
 
             Info<< "Number of eddies: "
@@ -969,7 +971,7 @@ void Foam::turbulentDFSEMInletFvPatchVectorField::updateCoeffs()
                 << endl;
 
             Info<< "Patch:" << patch().patch().name()
-                << " min/max(U):" << gMin(U) << ", " << gMax(U)
+                << " min/max(U):" << limits.min() << ", " << limits.max()
                 << endl;
 
             if (db().time().writeTime())
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/waveSurfacePressure/waveSurfacePressureFvPatchScalarField.C b/src/finiteVolume/fields/fvPatchFields/derived/waveSurfacePressure/waveSurfacePressureFvPatchScalarField.C
index 0eb44f87fa..e98f7f20b0 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/waveSurfacePressure/waveSurfacePressureFvPatchScalarField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/waveSurfacePressure/waveSurfacePressureFvPatchScalarField.C
@@ -210,8 +210,12 @@ void Foam::waveSurfacePressureFvPatchScalarField::updateCoeffs()
     }
 
 
-    Info<< "min/max zetap = " << gMin(zetap & nf()) << ", "
-        << gMax(zetap & nf()) << endl;
+    {
+        auto limits = gMinMax(zetap & nf());
+
+        Info<< "min/max zetap = "
+            << limits.min() << ", " << limits.max() << endl;
+    }
 
     // Update the surface pressure
     const uniformDimensionedVectorField& g =
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.C
index 06209c1334..9f339aecc5 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.C
@@ -186,12 +186,12 @@ tmp<surfaceScalarField> ddtScheme<Type>::fvcDdtPhiCoeff
 
     if (debug > 1)
     {
+        auto limits = gMinMax(ddtCouplingCoeff.primitiveField());
+        auto avg = gAverage(ddtCouplingCoeff.primitiveField());
+
         InfoInFunction
             << "ddtCouplingCoeff mean max min = "
-            << gAverage(ddtCouplingCoeff.primitiveField())
-            << " " << gMax(ddtCouplingCoeff.primitiveField())
-            << " " << gMin(ddtCouplingCoeff.primitiveField())
-            << endl;
+            << avg << ' ' << limits.max() << ' ' << limits.min() << endl;
     }
 
     return tddtCouplingCoeff;
@@ -267,12 +267,12 @@ tmp<surfaceScalarField> ddtScheme<Type>::fvcDdtPhiCoeffExperimental
 
     if (debug > 1)
     {
+        auto limits = gMinMax(ddtCouplingCoeff.primitiveField());
+        auto avg = gAverage(ddtCouplingCoeff.primitiveField());
+
         InfoInFunction
             << "ddtCouplingCoeff mean max min = "
-            << gAverage(ddtCouplingCoeff.primitiveField())
-            << " " << gMax(ddtCouplingCoeff.primitiveField())
-            << " " << gMin(ddtCouplingCoeff.primitiveField())
-            << endl;
+            << avg << ' ' << limits.max() << ' ' << limits.min() << endl;
     }
 
     return tddtCouplingCoeff;
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrad.C b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrad.C
index f259794f1f..0d6605cbfc 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrad.C
+++ b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrad.C
@@ -215,10 +215,13 @@ Foam::fv::cellLimitedGrad<Type, Limiter>::calcGrad
 
     if (fv::debug)
     {
+        auto limits = gMinMax(limiter);
+        auto avg = gAverage(limiter);
+
         Info<< "gradient limiter for: " << vsf.name()
-            << " max = " << gMax(limiter)
-            << " min = " << gMin(limiter)
-            << " average: " << gAverage(limiter) << endl;
+            << " min = " << limits.min()
+            << " max = " << limits.max()
+            << " average: " << avg << endl;
     }
 
     limitGradient(limiter, g);
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceLimitedGrad/faceLimitedGrads.C b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceLimitedGrad/faceLimitedGrads.C
index 5ef4594ef7..a6591c366e 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceLimitedGrad/faceLimitedGrads.C
+++ b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceLimitedGrad/faceLimitedGrads.C
@@ -162,10 +162,13 @@ Foam::fv::faceLimitedGrad<Foam::scalar>::calcGrad
 
     if (fv::debug)
     {
+        auto limits = gMinMax(limiter);
+        auto avg = gAverage(limiter);
+
         Info<< "gradient limiter for: " << vsf.name()
-            << " max = " << gMax(limiter)
-            << " min = " << gMin(limiter)
-            << " average: " << gAverage(limiter) << endl;
+            << " min = " << limits.min()
+            << " max = " << limits.max()
+            << " average: " << avg << endl;
     }
 
     g.primitiveFieldRef() *= limiter;
@@ -323,10 +326,13 @@ Foam::fv::faceLimitedGrad<Foam::vector>::calcGrad
 
     if (fv::debug)
     {
+        auto limits = gMinMax(limiter);
+        auto avg = gAverage(limiter);
+
         Info<< "gradient limiter for: " << vvf.name()
-            << " max = " << gMax(limiter)
-            << " min = " << gMin(limiter)
-            << " average: " << gAverage(limiter) << endl;
+            << " min = " << limits.min()
+            << " max = " << limits.max()
+            << " average: " << avg << endl;
     }
 
     g.primitiveFieldRef() *= limiter;
diff --git a/src/finiteVolume/fvMesh/fvGeometryScheme/averageNeighbour/averageNeighbourFvGeometryScheme.C b/src/finiteVolume/fvMesh/fvGeometryScheme/averageNeighbour/averageNeighbourFvGeometryScheme.C
index 6cfef0e944..feb35be095 100644
--- a/src/finiteVolume/fvMesh/fvGeometryScheme/averageNeighbour/averageNeighbourFvGeometryScheme.C
+++ b/src/finiteVolume/fvMesh/fvGeometryScheme/averageNeighbour/averageNeighbourFvGeometryScheme.C
@@ -814,6 +814,7 @@ void Foam::averageNeighbourFvGeometryScheme::movePoints()
                         Foam::acos(faceOrthogonality)
                     )
                 );
+
                 Pout<< "    iter:" << iter
                     << " nClipped:" << nClipped
                     << " average displacement:" << gAverage(magCorrection)
@@ -845,10 +846,13 @@ void Foam::averageNeighbourFvGeometryScheme::movePoints()
 
         if (debug)
         {
+            auto limits = gMinMax(cellWeight);
+            auto avg = gAverage(cellWeight);
+
             Pout<< "averageNeighbourFvGeometryScheme::movePoints() :"
                 << " averageNeighbour weight"
-                << " max:" << gMax(cellWeight) << " min:" << gMin(cellWeight)
-                << " average:" << gAverage(cellWeight) << endl;
+                << " min:" << limits.min() << " max:" << limits.max()
+                << " average:" << avg << endl;
 
             // Dump lines from old to new location
             const fileName tp(mesh_.time().timePath());
diff --git a/src/finiteVolume/fvMesh/fvGeometryScheme/highAspectRatio/cellAspectRatio.C b/src/finiteVolume/fvMesh/fvGeometryScheme/highAspectRatio/cellAspectRatio.C
index efa5a01fd9..ce87146904 100644
--- a/src/finiteVolume/fvMesh/fvGeometryScheme/highAspectRatio/cellAspectRatio.C
+++ b/src/finiteVolume/fvMesh/fvGeometryScheme/highAspectRatio/cellAspectRatio.C
@@ -120,8 +120,11 @@ void Foam::cellAspectRatio::calcAspectRatio()
 
     if (debug)
     {
-        InfoInFunction << "Calculated cell aspect ratio min:" << gMin(aRatio)
-            << " max:" << gMax(aRatio) << " average:" << gAverage(aRatio)
+        auto limits = gMinMax(aRatio);
+        auto avg = gAverage(aRatio);
+
+        InfoInFunction << "Calculated cell aspect ratio min:" << limits.min()
+            << " max:" << limits.max() << " average:" << avg
             << endl;
     }
 }
diff --git a/src/finiteVolume/fvMesh/fvGeometryScheme/highAspectRatio/highAspectRatioFvGeometryScheme.C b/src/finiteVolume/fvMesh/fvGeometryScheme/highAspectRatio/highAspectRatioFvGeometryScheme.C
index 0f64298abc..6e60a2eac3 100644
--- a/src/finiteVolume/fvMesh/fvGeometryScheme/highAspectRatio/highAspectRatioFvGeometryScheme.C
+++ b/src/finiteVolume/fvMesh/fvGeometryScheme/highAspectRatio/highAspectRatioFvGeometryScheme.C
@@ -434,10 +434,13 @@ void Foam::highAspectRatioFvGeometryScheme::movePoints()
 
         if (debug)
         {
+            auto limits = gMinMax(cellWeight);
+            auto avg = gAverage(cellWeight);
+
             Pout<< "highAspectRatioFvGeometryScheme::movePoints() :"
                 << " highAspectRatio weight"
-                << " max:" << gMax(cellWeight) << " min:" << gMin(cellWeight)
-                << " average:" << gAverage(cellWeight) << endl;
+                << " min:" << limits.max() << " max:" << limits.max()
+                << " average:" << avg << endl;
         }
 
         vectorField faceAreas(mesh_.faceAreas());
diff --git a/src/functionObjects/field/AMIWeights/AMIWeights.C b/src/functionObjects/field/AMIWeights/AMIWeights.C
index 7ec095d55d..a3a49b20aa 100644
--- a/src/functionObjects/field/AMIWeights/AMIWeights.C
+++ b/src/functionObjects/field/AMIWeights/AMIWeights.C
@@ -88,52 +88,52 @@ void Foam::functionObjects::AMIWeights::reportPatch
     const Switch distributed = pp.AMI().distributed();
 
     const scalarField& srcWeightsSum = pp.AMI().srcWeightsSum();
-    const scalar srcMinWeight = gMin(srcWeightsSum);
-    const scalar srcMaxWeight = gMax(srcWeightsSum);
-    const scalar srcAveWeight = gAverage(srcWeightsSum);
+    const auto srcWeightLimits = gMinMax(srcWeightsSum);
+    const auto srcWeightAvg = gAverage(srcWeightsSum);
 
     const labelListList& srcAddress = pp.AMI().srcAddress();
-    label srcMinNbr = labelMax;
-    label srcMaxNbr = labelMin;
-    scalar srcAveNbr = 0;
+
+    labelMinMax srcNbrLimits(labelMax, labelMin);
+    scalar srcNbrAvg(0);
     for (const labelList& srcFace : srcAddress)
     {
         const label n = srcFace.size();
-        srcAveNbr += n;
-        srcMinNbr = min(srcMinNbr, n);
-        srcMaxNbr = max(srcMaxNbr, n);
+
+        srcNbrAvg += n;
+        srcNbrLimits.add(n);
     }
 
-    reduce(srcMinNbr, minOp<label>());
-    reduce(srcMaxNbr, maxOp<label>());
+    {
+        reduce(srcNbrLimits, sumOp<labelMinMax>());
 
-    srcAveNbr =
-        returnReduce(srcAveNbr, sumOp<scalar>())
-       /(returnReduce(srcAddress.size(), sumOp<scalar>()) + ROOTVSMALL);
+        label count = srcAddress.size();
+        sumReduce(srcNbrAvg, count);
+        srcNbrAvg /= (count + ROOTVSMALL);
+    }
 
     const scalarField& tgtWeightsSum = pp.AMI().tgtWeightsSum();
-    const scalar tgtMinWeight = gMin(tgtWeightsSum);
-    const scalar tgtMaxWeight = gMax(tgtWeightsSum);
-    const scalar tgtAveWeight = gAverage(tgtWeightsSum);
+    const auto tgtWeightLimits = gMinMax(tgtWeightsSum);
+    const auto tgtWeightAvg = gAverage(tgtWeightsSum);
 
     const labelListList& tgtAddress = pp.AMI().tgtAddress();
-    label tgtMinNbr = labelMax;
-    label tgtMaxNbr = labelMin;
-    scalar tgtAveNbr = 0;
+
+    labelMinMax tgtNbrLimits(labelMax, labelMin);
+    scalar tgtNbrAvg(0);
     for (const labelList& tgtFace : tgtAddress)
     {
         const label n = tgtFace.size();
-        tgtAveNbr += n;
-        tgtMinNbr = min(tgtMinNbr, n);
-        tgtMaxNbr = max(tgtMaxNbr, n);
+
+        tgtNbrAvg += n;
+        tgtNbrLimits.add(n);
     }
 
-    reduce(tgtMinNbr, minOp<label>());
-    reduce(tgtMaxNbr, maxOp<label>());
+    {
+        reduce(tgtNbrLimits, sumOp<labelMinMax>());
 
-    tgtAveNbr =
-        returnReduce(tgtAveNbr, sumOp<scalar>())
-       /(returnReduce(tgtAddress.size(), sumOp<scalar>()) + ROOTVSMALL);
+        label count = tgtAddress.size();
+        sumReduce(tgtNbrAvg, count);
+        tgtNbrAvg /= (count + ROOTVSMALL);
+    }
 
     file()
         << mesh_.time().timeName() << tab
@@ -147,18 +147,18 @@ void Foam::functionObjects::AMIWeights::reportPatch
     }
 
     file()
-        << srcMinWeight << tab
-        << srcMaxWeight << tab
-        << srcAveWeight << tab
-        << srcMinNbr << tab
-        << srcMaxNbr << tab
-        << srcAveNbr << tab
-        << tgtMinWeight << tab
-        << tgtMaxWeight << tab
-        << tgtAveWeight << tab
-        << tgtMinNbr << tab
-        << tgtMaxNbr << tab
-        << tgtAveNbr << tab
+        << srcWeightLimits.min() << tab
+        << srcWeightLimits.max() << tab
+        << srcWeightAvg << tab
+        << srcNbrLimits.min() << tab
+        << srcNbrLimits.max() << tab
+        << srcNbrAvg << tab
+        << tgtWeightLimits.min() << tab
+        << tgtWeightLimits.max() << tab
+        << tgtWeightAvg << tab
+        << tgtNbrLimits.min() << tab
+        << tgtNbrLimits.max() << tab
+        << tgtNbrAvg
         << endl;
 
     Log << "    Patches: " << nl
@@ -176,34 +176,34 @@ void Foam::functionObjects::AMIWeights::reportPatch
 
     Log << "                     | " << setw(w) << pp.name()
         << " | " << setw(w) << nbrPatchName << " | " << nl
-        << "        min(weight)  | " << setw(w) << srcMinWeight
-        << " | " << setw(w) << tgtMinWeight << " | " << nl
-        << "        max(weight)  | " << setw(w) << srcMaxWeight
-        << " | " << setw(w) << tgtMaxWeight << " | " << nl
-        << "        ave(weight)  | " << setw(w) << srcAveWeight
-        << " | " << setw(w) << tgtAveWeight << " | " << nl
-        << "        min(address) | " << setw(w) << srcMinNbr
-        << " | " << setw(w) << tgtMinNbr << " | " << nl
-        << "        max(address) | " << setw(w) << srcMaxNbr
-        << " | " << setw(w) << tgtMaxNbr << " | " << nl
-        << "        ave(address) | " << setw(w) << srcAveNbr
-        << " | " << setw(w) << tgtAveNbr << " | " << nl
+        << "        min(weight)  | " << setw(w) << srcWeightLimits.min()
+        << " | " << setw(w) << tgtWeightLimits.min() << " | " << nl
+        << "        max(weight)  | " << setw(w) << srcWeightLimits.max()
+        << " | " << setw(w) << tgtWeightLimits.max() << " | " << nl
+        << "        ave(weight)  | " << setw(w) << srcWeightAvg
+        << " | " << setw(w) << tgtWeightAvg << " | " << nl
+        << "        min(address) | " << setw(w) << srcNbrLimits.min()
+        << " | " << setw(w) << tgtNbrLimits.min() << " | " << nl
+        << "        max(address) | " << setw(w) << srcNbrLimits.max()
+        << " | " << setw(w) << tgtNbrLimits.max() << " | " << nl
+        << "        ave(address) | " << setw(w) << srcNbrAvg
+        << " | " << setw(w) << tgtNbrAvg << " | " << nl
         << endl;
 
     setResult(pp.name() + ":src", pp.name());
     setResult(pp.name() + ":tgt", nbrPatchName);
-    setResult(pp.name() + ":src:min(weight)", srcMinWeight);
-    setResult(pp.name() + ":src:max(weight)", srcMaxWeight);
-    setResult(pp.name() + ":src:ave(weight)", srcAveWeight);
-    setResult(pp.name() + ":src:min(address)", srcMinNbr);
-    setResult(pp.name() + ":src:max(address)", srcMaxNbr);
-    setResult(pp.name() + ":src:ave(address)", srcAveNbr);
-    setResult(pp.name() + ":tgt:min(weight)", tgtMinWeight);
-    setResult(pp.name() + ":tgt:max(weight)", tgtMaxWeight);
-    setResult(pp.name() + ":tgt:ave(weight)", tgtAveWeight);
-    setResult(pp.name() + ":tgt:min(address)", tgtMinNbr);
-    setResult(pp.name() + ":tgt:max(address)", tgtMaxNbr);
-    setResult(pp.name() + ":tgt:ave(address)", tgtAveNbr);
+    setResult(pp.name() + ":src:min(weight)", srcWeightLimits.min());
+    setResult(pp.name() + ":src:max(weight)", srcWeightLimits.max());
+    setResult(pp.name() + ":src:ave(weight)", srcWeightAvg);
+    setResult(pp.name() + ":src:min(address)", srcNbrLimits.min());
+    setResult(pp.name() + ":src:max(address)", srcNbrLimits.max());
+    setResult(pp.name() + ":src:ave(address)", srcNbrAvg);
+    setResult(pp.name() + ":tgt:min(weight)", tgtWeightLimits.min());
+    setResult(pp.name() + ":tgt:max(weight)", tgtWeightLimits.max());
+    setResult(pp.name() + ":tgt:ave(weight)", tgtWeightAvg);
+    setResult(pp.name() + ":tgt:min(address)", tgtNbrLimits.min());
+    setResult(pp.name() + ":tgt:max(address)", tgtNbrLimits.max());
+    setResult(pp.name() + ":tgt:ave(address)", tgtNbrAvg);
 }
 
 
diff --git a/src/functionObjects/field/yPlus/yPlus.C b/src/functionObjects/field/yPlus/yPlus.C
index f9b5f3a3b5..f767e00f5b 100644
--- a/src/functionObjects/field/yPlus/yPlus.C
+++ b/src/functionObjects/field/yPlus/yPlus.C
@@ -211,25 +211,24 @@ bool Foam::functionObjects::yPlus::write()
         {
             const scalarField& yPlusp = yPlusBf[patchi];
 
-            const scalar minYplus = gMin(yPlusp);
-            const scalar maxYplus = gMax(yPlusp);
-            const scalar avgYplus = gAverage(yPlusp);
+            auto limits = gMinMax(yPlusp);
+            auto avg = gAverage(yPlusp);
 
             if (UPstream::master())
             {
                 writeCurrentTime(file());
                 file()
                     << token::TAB << patch.name()
-                    << token::TAB << minYplus
-                    << token::TAB << maxYplus
-                    << token::TAB << avgYplus
+                    << token::TAB << limits.min()
+                    << token::TAB << limits.max()
+                    << token::TAB << avg
                     << endl;
             }
 
             Log << "    patch " << patch.name()
-                << " y+ : min = " << minYplus
-                << ", max = " << maxYplus
-                << ", average = " << avgYplus << endl;
+                << " y+ : min = " << limits.min()
+                << ", max = " << limits.max()
+                << ", average = " << avg << endl;
         }
     }
 
diff --git a/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.C b/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.C
index 73f8923439..d57ddb09a9 100644
--- a/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.C
+++ b/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.C
@@ -552,9 +552,12 @@ void Foam::timeVaryingMappedFixedValuePointPatchField<Type>::updateCoeffs()
 
     if (debug)
     {
-        Pout<< "updateCoeffs : set fixedValue to min:" << gMin(*this)
-            << " max:" << gMax(*this)
-            << " avg:" << gAverage(*this) << endl;
+        auto limits = gMinMax(*this);
+        auto avg = gAverage(*this);
+
+        Pout<< "updateCoeffs : set fixedValue to min:" << limits.min()
+            << " max:" << limits.max()
+            << " avg:" << avg << endl;
     }
 
     fixedValuePointPatchField<Type>::updateCoeffs();
diff --git a/src/fvOptions/sources/derived/heatExchangerSource/heatExchangerModels/effectivenessTable/effectivenessTable.C b/src/fvOptions/sources/derived/heatExchangerSource/heatExchangerModels/effectivenessTable/effectivenessTable.C
index 08af91be26..b5e42d1af5 100644
--- a/src/fvOptions/sources/derived/heatExchangerSource/heatExchangerModels/effectivenessTable/effectivenessTable.C
+++ b/src/fvOptions/sources/derived/heatExchangerSource/heatExchangerModels/effectivenessTable/effectivenessTable.C
@@ -195,23 +195,23 @@ Foam::heatExchangerModels::effectivenessTable::energyDensity
         secondaryInletT_ += targetQdotRelax_*dT;
     }
 
-    const scalarField TCells(T, cells);
-    scalarField deltaTCells(cells.size(), Zero);
+    // start with a copy
+    scalarField deltaTCells(T, cells);
     Tref_ = 0;
     if (Qt_ > 0)
     {
-        Tref_ = gMax(TCells);
-        forAll(deltaTCells, i)
+        Tref_ = gMax(deltaTCells);
+        for (scalar& delta : deltaTCells)
         {
-            deltaTCells[i] = max(Tref_ - TCells[i], scalar(0));
+            delta = Foam::max(Tref_ - delta, scalar(0));
         }
     }
     else
     {
-        Tref_ = gMin(TCells);
-        forAll(deltaTCells, i)
+        Tref_ = gMin(deltaTCells);
+        for (scalar& delta : deltaTCells)
         {
-            deltaTCells[i] = max(TCells[i] - Tref_, scalar(0));
+            delta = Foam::max(delta - Tref_, scalar(0));
         }
     }
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
index 776d3af811..a513fb2ffb 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
@@ -950,13 +950,12 @@ void Foam::KinematicCloud<CloudType>::info()
             alpha().write();
         }
 
-        const scalar alphaMin = gMin(alpha().primitiveField());
-        const scalar alphaMax = gMax(alpha().primitiveField());
+        auto limits = gMinMax(alpha().primitiveField());
 
-        Log_<< "    Min cell volume fraction        = " << alphaMin << nl
-            << "    Max cell volume fraction        = " << alphaMax << endl;
+        Log_<< "    Min cell volume fraction        = " << limits.min() << nl
+            << "    Max cell volume fraction        = " << limits.max() << endl;
 
-        if (alphaMax < SMALL)
+        if (limits.max() < SMALL)
         {
             return;
         }
@@ -969,7 +968,7 @@ void Foam::KinematicCloud<CloudType>::info()
 
             if (n > 0)
             {
-                const scalar nPack = n*alphaMax/alpha()[celli];
+                const scalar nPack = n*limits.max()/alpha()[celli];
 
                 if (nPack < nMin)
                 {
diff --git a/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloud.C b/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloud.C
index 14dcb633ed..3e7bb38560 100644
--- a/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloud.C
@@ -283,13 +283,12 @@ void Foam::MPPICCloud<CloudType>::info()
 
     tmp<volScalarField> alpha = this->theta();
 
-    const scalar alphaMin = gMin(alpha().primitiveField());
-    const scalar alphaMax = gMax(alpha().primitiveField());
+    auto limits = gMinMax(alpha().primitiveField());
 
-    Log_ << "    Min cell volume fraction        = " << alphaMin << nl
-         << "    Max cell volume fraction        = " << alphaMax << endl;
+    Log_ << "    Min cell volume fraction        = " << limits.min() << nl
+         << "    Max cell volume fraction        = " << limits.max() << endl;
 
-    if (alphaMax < SMALL)
+    if (limits.max() < SMALL)
     {
         return;
     }
@@ -302,7 +301,7 @@ void Foam::MPPICCloud<CloudType>::info()
 
         if (n > 0)
         {
-            const scalar nPack = n*alphaMax/alpha()[celli];
+            const scalar nPack = n*limits.max()/alpha()[celli];
 
             if (nPack < nMin)
             {
diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.C
index d362fabd28..181a588f87 100644
--- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.C
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.C
@@ -1676,8 +1676,7 @@ void Foam::snappyLayerDriver::calculateLayerThickness
         GREAT               // null value
     );
 
-    //Info<< "calculateLayerThickness : min:" << gMin(thickness)
-    //    << " max:" << gMax(thickness) << endl;
+    //Info<< "calculateLayerThickness : " << gMinMax(thickness) << endl;
 
     // Print a bit
     {
diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.C
index a37f5de341..d4372068a6 100644
--- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.C
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.C
@@ -2087,10 +2087,12 @@ Foam::vectorField Foam::snappySnapDriver::calcNearestSurface
 
             scalarField magDisp(mag(patchDisp));
 
+            auto limits = gMinMax(magDisp);
+
             Info<< "Wanted displacement : average:"
                 <<  meshRefinement::gAverage(isPatchMasterPoint, magDisp)
-                << " min:" << gMin(magDisp)
-                << " max:" << gMax(magDisp) << endl;
+                << " min:" << limits.min()
+                << " max:" << limits.max() << endl;
         }
     }
 
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
index 4d2ab16688..472aaff449 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
@@ -279,37 +279,31 @@ void Foam::AMIInterpolation::normaliseWeights
         }
     }
 
-    if (output && comm != -1)
+    if (output && comm != -1 && returnReduceOr(wght.size(), comm))
     {
-        // Note: change global communicator since gMin,gAverage etc don't
-        // support user communicator
-        const label oldWorldComm(UPstream::worldComm);
-        UPstream::worldComm = comm;
+        auto limits = gMinMax(wghtSum, comm);
+        auto avg = gAverage(wghtSum, comm);
 
-        if (returnReduceOr(wght.size()))
+        label nLow =
+            returnReduce(nLowWeight, sumOp<label>(), UPstream::msgType(), comm);
+
+        Info.masterStream(comm)
+            << indent
+            << "AMI: Patch " << patchName
+            << " sum(weights)"
+            << " min:" << limits.min()
+            << " max:" << limits.max()
+            << " average:" << avg << nl;
+
+        if (nLow)
         {
             Info.masterStream(comm)
                 << indent
                 << "AMI: Patch " << patchName
-                << " sum(weights)"
-                << " min:" << gMin(wghtSum)
-                << " max:" << gMax(wghtSum)
-                << " average:" << gAverage(wghtSum) << nl;
-
-            const label nLow = returnReduce(nLowWeight, sumOp<label>());
-
-            if (nLow)
-            {
-                Info.masterStream(comm)
-                    << indent
-                    << "AMI: Patch " << patchName
-                    << " identified " << nLow
-                    << " faces with weights less than " << lowWeightTol
-                    << endl;
-            }
+                << " identified " << nLow
+                << " faces with weights less than " << lowWeightTol
+                << endl;
         }
-
-        UPstream::worldComm = oldWorldComm;
     }
 }
 
diff --git a/src/meshTools/AMIInterpolation/patches/cyclicPeriodicAMI/cyclicPeriodicAMIPolyPatch/cyclicPeriodicAMIPolyPatch.C b/src/meshTools/AMIInterpolation/patches/cyclicPeriodicAMI/cyclicPeriodicAMIPolyPatch/cyclicPeriodicAMIPolyPatch.C
index 410a797230..6aba9c7b53 100644
--- a/src/meshTools/AMIInterpolation/patches/cyclicPeriodicAMI/cyclicPeriodicAMIPolyPatch/cyclicPeriodicAMIPolyPatch.C
+++ b/src/meshTools/AMIInterpolation/patches/cyclicPeriodicAMI/cyclicPeriodicAMIPolyPatch/cyclicPeriodicAMIPolyPatch.C
@@ -513,18 +513,29 @@ void Foam::cyclicPeriodicAMIPolyPatch::resetAMI() const
                 tgtWghtSum[faceI] = sum(AMIPtr_->tgtWeights()[faceI]);
             }
 
-            Info<< indent
-                << "AMI: Patch " << name()
-                << " sum(weights)"
-                << " min:" << gMin(srcWghtSum)
-                << " max:" << gMax(srcWghtSum)
-                << " average:" << gAverage(srcWghtSum) << nl;
-            Info<< indent
-                << "AMI: Patch " << neighbPatch().name()
-                << " sum(weights)"
-                << " min:" << gMin(tgtWghtSum)
-                << " max:" << gMax(tgtWghtSum)
-                << " average:" << gAverage(tgtWghtSum) << nl;
+            {
+                auto limits = gMinMax(srcWghtSum);
+                auto avg = gAverage(srcWghtSum);
+
+                Info<< indent
+                    << "AMI: Patch " << name()
+                    << " sum(weights)"
+                    << " min:" << limits.min()
+                    << " max:" << limits.max()
+                    << " average:" << avg << nl;
+            }
+
+            {
+                auto limits = gMinMax(tgtWghtSum);
+                auto avg = gAverage(tgtWghtSum);
+
+                Info<< indent
+                    << "AMI: Patch " << neighbPatch().name()
+                    << " sum(weights)"
+                    << " min:" << limits.min()
+                    << " max:" << limits.max()
+                    << " average:" << avg << nl;
+            }
         }
     }
 }
diff --git a/src/meshTools/PatchFunction1/MappedFile/MappedFile.C b/src/meshTools/PatchFunction1/MappedFile/MappedFile.C
index ef0856fe7b..2e6adb35f5 100644
--- a/src/meshTools/PatchFunction1/MappedFile/MappedFile.C
+++ b/src/meshTools/PatchFunction1/MappedFile/MappedFile.C
@@ -747,9 +747,13 @@ Foam::PatchFunction1Types::MappedFile<Type>::value
 
     if (debug)
     {
-        Pout<< "MappedFile<Type>::value : set fixedValue to min:" << gMin(fld)
-            << " max:" << gMax(fld)
-            << " avg:" << gAverage(fld) << endl;
+        auto limits = gMinMax(fld);
+        auto avg = gAverage(fld);
+
+        Pout<< "MappedFile<Type>::value : set fixedValue to min:"
+            << limits.min()
+            << " max:" << limits.max()
+            << " avg:" << avg << endl;
     }
 
     return this->transform(tfld);
diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/regularisation/regularisationPDE/Helmoltz/Helmholtz.C b/src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/regularisation/regularisationPDE/Helmoltz/Helmholtz.C
index 39980eece8..ced36121cb 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/regularisation/regularisationPDE/Helmoltz/Helmholtz.C
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/regularisation/regularisationPDE/Helmoltz/Helmholtz.C
@@ -231,7 +231,7 @@ void Foam::Helmholtz::regularise
             );
             // Map result back to original field
             result.rmap(resultSub, cellMap);
-            Info<< "min max " << gMin(result) << " " << gMax(result) << endl;
+            Info<< "min max " << gMinMax(result) << endl;
             return;
         }
     }
diff --git a/src/phaseSystemModels/reactingEuler/multiphaseSystem/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C b/src/phaseSystemModels/reactingEuler/multiphaseSystem/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
index e3a2be4b9a..ac4ee7ca4c 100644
--- a/src/phaseSystemModels/reactingEuler/multiphaseSystem/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
+++ b/src/phaseSystemModels/reactingEuler/multiphaseSystem/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
@@ -430,21 +430,29 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctInterfaceThermo()
 
         Tf = (H1*T1 + H2*T2 + iDmdtNew*L)/(H1 + H2);
 
-        Info<< "Tf." << pair.name()
-            << ": min = " << gMin(Tf.primitiveField())
-            << ", mean = " << gAverage(Tf.primitiveField())
-            << ", max = " << gMax(Tf.primitiveField())
-            << endl;
+        {
+            auto limits = gMinMax(Tf.primitiveField());
+            auto avg = gAverage(Tf.primitiveField());
+
+            Info<< "Tf." << pair.name()
+                << ": min = " << limits.min()
+                << ", mean = " << avg
+                << ", max = " << limits.max()
+                << endl;
+        }
 
         scalar iDmdtRelax(this->mesh().fieldRelaxationFactor("iDmdt"));
         iDmdt = (1 - iDmdtRelax)*iDmdt + iDmdtRelax*iDmdtNew;
 
         if (phaseChange_)
         {
+            auto limits = gMinMax(iDmdt.primitiveField());
+            auto avg = gAverage(iDmdt.primitiveField());
+
             Info<< "iDmdt." << pair.name()
-                << ": min = " << gMin(iDmdt.primitiveField())
-                << ", mean = " << gAverage(iDmdt.primitiveField())
-                << ", max = " << gMax(iDmdt.primitiveField())
+                << ": min = " << limits.min()
+                << ", mean = " << avg
+                << ", max = " << limits.max()
                 << ", integral = " << fvc::domainIntegrate(iDmdt).value()
                 << endl;
         }
@@ -525,10 +533,13 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctInterfaceThermo()
 
         if (wallBoilingActive)
         {
+            auto limits = gMinMax(wDmdt.primitiveField());
+            auto avg = gAverage(wDmdt.primitiveField());
+
             Info<< "wDmdt." << pair.name()
-                << ": min = " << gMin(wDmdt.primitiveField())
-                << ", mean = " << gAverage(wDmdt.primitiveField())
-                << ", max = " << gMax(wDmdt.primitiveField())
+                << ": min = " << limits.min()
+                << ", mean = " << avg
+                << ", max = " << limits.max()
                 << ", integral = " << fvc::domainIntegrate(wDmdt).value()
                 << endl;
         }
diff --git a/src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C b/src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C
index 56cdee16f8..efca2b963c 100644
--- a/src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C
+++ b/src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C
@@ -505,13 +505,12 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
             {
                 Info<< "alphat for vapour : " << nl << endl;
 
-                Info<< "  alphatEffv: " << gMin(vaporw*(*this + alphaw))
-                    << " - " << gMax(vaporw*(*this + alphaw)) << endl;
+                Info<< "  alphatEffv: " << gMinMax(vaporw*(*this + alphaw))
+                    << endl;
 
                 const scalarField qEff(vaporw*(*this + alphaw)*hewv.snGrad());
 
-                 Info<< "  qEffVap: " << gMin(qEff) << " - "
-                     << gMax(qEff) << endl;
+                Info<< "  qEffVap: " << gMinMax(qEff) << endl;
 
                 scalar Qeff = gWeightedSum(patch().magSf(), qEff);
                 Info<< " Effective heat transfer rate to vapor:" << Qeff
@@ -1096,24 +1095,16 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                     fLiquid*liquidw*(*this + alphaw)*hew.snGrad()
                 );
 
-                Info<< "alphat for liquid:  " <<  nl << endl;
+                Info<< "alphat for liquid:  " <<  nl << nl;
+                Info<< "  qEffLiq: " << gMinMax(qEff) << nl;
+                Info<< "  alphatl: " << gMinMax(*this) << nl;
+                Info<< "  dmdt: " << gMinMax(dmdt_) << nl;
+                Info<< "  alphatlEff: "
+                    << gMinMax(liquidw*(*this + alphaw)) << nl;
 
-                Info<< "  qEffLiq: " << gMin(qEff) << " - "
-                    << gMax(qEff) << endl;
-
-
-                Info<< "  alphatl: " << gMin((*this)) << " - "
-                    << gMax((*this)) << endl;
-
-                Info<< "  dmdt: " << gMin((dmdt_)) << " - "
-                    << gMax((dmdt_)) << endl;
-
-                Info<< "  alphatlEff: " << gMin(liquidw*(*this + alphaw))
-                    << " - " << gMax(liquidw*(*this + alphaw)) << endl;
-
-                scalar Qeff = gWeightedSum(patch().magSf(), qEff);
-                Info<< " Effective heat transfer rate to liquid: " << Qeff
-                    << endl << nl;
+                Info<< " Effective heat transfer rate to liquid: "
+                    << gWeightedSum(patch().magSf(), qEff)
+                    << nl << endl;
 
                 if (debug == 2)
                 {
diff --git a/src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/fixedMultiPhaseHeatFlux/fixedMultiPhaseHeatFluxFvPatchScalarField.C b/src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/fixedMultiPhaseHeatFlux/fixedMultiPhaseHeatFluxFvPatchScalarField.C
index 3f32741033..9139ea3a54 100644
--- a/src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/fixedMultiPhaseHeatFlux/fixedMultiPhaseHeatFluxFvPatchScalarField.C
+++ b/src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/fixedMultiPhaseHeatFlux/fixedMultiPhaseHeatFluxFvPatchScalarField.C
@@ -147,8 +147,11 @@ void Foam::fixedMultiPhaseHeatFluxFvPatchScalarField::updateCoeffs()
             const scalarField q0(T.snGrad()*alpha*kappaEff);
             Q += q0;
 
+            auto limits = gMinMax(q0);
+
             Info<< patch().name() << " " << phase.name()
-                << ": Heat flux " << gMin(q0) << " - " << gMax(q0) << endl;
+                << ": Heat flux "
+                << limits.min() << " - " << limits.max() << endl;
         }
 
         A += T.patchInternalField()*alpha*kappaEff*patch().deltaCoeffs();
@@ -157,7 +160,7 @@ void Foam::fixedMultiPhaseHeatFluxFvPatchScalarField::updateCoeffs()
 
     if (debug)
     {
-        MinMax<scalar> limits = gMinMax(Q);
+        auto limits = gMinMax(Q);
 
         Info<< patch().name() << " " << ": overall heat flux "
             << limits.min() << " - " << limits.max() << " W/m2, power: "
diff --git a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C
index acadb10b20..2c81e9ac7c 100644
--- a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C
+++ b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C
@@ -694,11 +694,10 @@ void reactingOneDim::evolveRegion()
 
     solidThermo_->correct();
 
+    auto limits = gMinMax(solidThermo_->T().primitiveField());
+
     Info<< "pyrolysis min/max(T) = "
-        << gMin(solidThermo_->T().primitiveField())
-        << ", "
-        << gMax(solidThermo_->T().primitiveField())
-        << endl;
+        << limits.min() << ", " << limits.max() << endl;
 }
 
 
diff --git a/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C b/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C
index ba9bf0b0cb..b62a18da60 100644
--- a/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C
+++ b/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C
@@ -697,10 +697,11 @@ void thermoSingleLayer::info()
 
     const scalarField& Tinternal = T_;
 
+    auto limits = gMinMax(Tinternal);
+    auto avg = gAverage(Tinternal);
+
     Info<< indent << "min/mean/max(T)    = "
-        << gMin(Tinternal) << ", "
-        << gAverage(Tinternal) << ", "
-        << gMax(Tinternal) << nl;
+        << limits.min() << ", " << avg << ", " << limits.max() << nl;
 
     phaseChange_->info(Info);
 }
diff --git a/src/thermoTools/derivedFvPatchFields/lumpedMassWallTemperature/lumpedMassWallTemperatureFvPatchScalarField.C b/src/thermoTools/derivedFvPatchFields/lumpedMassWallTemperature/lumpedMassWallTemperatureFvPatchScalarField.C
index df66f3288c..215d17698f 100644
--- a/src/thermoTools/derivedFvPatchFields/lumpedMassWallTemperature/lumpedMassWallTemperatureFvPatchScalarField.C
+++ b/src/thermoTools/derivedFvPatchFields/lumpedMassWallTemperature/lumpedMassWallTemperatureFvPatchScalarField.C
@@ -195,14 +195,17 @@ void Foam::lumpedMassWallTemperatureFvPatchScalarField::updateCoeffs()
             }
         }
 
+        auto limits = gMinMax(*this);
+        auto avg = gAverage(*this);
+
         Info<< patch().boundaryMesh().mesh().name() << ':'
             << patch().name() << ':'
             << this->internalField().name() << " :"
             << " heat transfer rate:" << Q
             << " wall temperature "
-            << " min:" << gMin(*this)
-            << " max:" << gMax(*this)
-            << " avg:" << gAverage(*this)
+            << " min:" << limits.min()
+            << " max:" << limits.max()
+            << " avg:" << avg
             << " Qin [W]:" << Qin
             << " Qout [W]:" << Qout
             << endl;
diff --git a/src/thermoTools/derivedFvPatchFields/thermalBaffle1D/thermalBaffle1DFvPatchScalarField.C b/src/thermoTools/derivedFvPatchFields/thermalBaffle1D/thermalBaffle1DFvPatchScalarField.C
index bb3532b43d..c73f38cd62 100644
--- a/src/thermoTools/derivedFvPatchFields/thermalBaffle1D/thermalBaffle1DFvPatchScalarField.C
+++ b/src/thermoTools/derivedFvPatchFields/thermalBaffle1D/thermalBaffle1DFvPatchScalarField.C
@@ -400,6 +400,10 @@ void thermalBaffle1DFvPatchScalarField<solidType>::updateCoeffs()
         if (debug)
         {
             scalar Q = gAverage(kappaw*snGrad());
+
+            auto limits = gMinMax(*this);
+            auto avg = gAverage(*this);
+
             Info<< patch().boundaryMesh().mesh().name() << ':'
                 << patch().name() << ':'
                 << this->internalField().name() << " <- "
@@ -407,9 +411,9 @@ void thermalBaffle1DFvPatchScalarField<solidType>::updateCoeffs()
                 << this->internalField().name() << " :"
                 << " heat[W]:" << Q
                 << " walltemperature "
-                << " min:" << gMin(*this)
-                << " max:" << gMax(*this)
-                << " avg:" << gAverage(*this)
+                << " min:" << limits.min()
+                << " max:" << limits.max()
+                << " avg:" << avg
                 << endl;
         }
     }
diff --git a/src/thermoTools/derivedFvPatchFields/turbulentTemperatureCoupledBaffleMixed/turbulentTemperatureCoupledBaffleMixedFvPatchScalarField.C b/src/thermoTools/derivedFvPatchFields/turbulentTemperatureCoupledBaffleMixed/turbulentTemperatureCoupledBaffleMixedFvPatchScalarField.C
index d2066bec6d..dcd7821636 100644
--- a/src/thermoTools/derivedFvPatchFields/turbulentTemperatureCoupledBaffleMixed/turbulentTemperatureCoupledBaffleMixedFvPatchScalarField.C
+++ b/src/thermoTools/derivedFvPatchFields/turbulentTemperatureCoupledBaffleMixed/turbulentTemperatureCoupledBaffleMixedFvPatchScalarField.C
@@ -385,6 +385,9 @@ void turbulentTemperatureCoupledBaffleMixedFvPatchScalarField::updateCoeffs()
     {
         scalar Q = gSum(kappaTp*patch().magSf()*snGrad());
 
+        auto limits = gMinMax(*this);
+        auto avg = gAverage(*this);
+
         Info<< patch().boundaryMesh().mesh().name() << ':'
             << patch().name() << ':'
             << this->internalField().name() << " <- "
@@ -393,9 +396,9 @@ void turbulentTemperatureCoupledBaffleMixedFvPatchScalarField::updateCoeffs()
             << this->internalField().name() << " :"
             << " heat transfer rate:" << Q
             << " walltemperature "
-            << " min:" << gMin(*this)
-            << " max:" << gMax(*this)
-            << " avg:" << gAverage(*this)
+            << " min:" << limits.min()
+            << " max:" << limits.max()
+            << " avg:" << avg
             << endl;
     }
 
diff --git a/src/thermophysicalModels/reactionThermo/derivedFvPatchFields/enthalpySorption/enthalpySorptionFvPatchScalarField.C b/src/thermophysicalModels/reactionThermo/derivedFvPatchFields/enthalpySorption/enthalpySorptionFvPatchScalarField.C
index ca9bee458f..3a8d0004ee 100644
--- a/src/thermophysicalModels/reactionThermo/derivedFvPatchFields/enthalpySorption/enthalpySorptionFvPatchScalarField.C
+++ b/src/thermophysicalModels/reactionThermo/derivedFvPatchFields/enthalpySorption/enthalpySorptionFvPatchScalarField.C
@@ -234,8 +234,10 @@ patchSource() const
 
     if (debug)
     {
+        auto limits = gMinMax(dhdt);
+
         Info<< " Patch enthalpy rate min/max [J/m3/sec]: "
-            << gMin(dhdt) << " - " << gMax(dhdt) << endl;
+            << limits.min() << " - " << limits.max() << endl;
     }
 
     return tmp<scalarField>::New(dhdt);
@@ -282,8 +284,10 @@ void Foam::enthalpySorptionFvPatchScalarField::updateCoeffs()
 
     if (debug)
     {
+        auto limits = gMinMax(dhdt_);
+
         Info<< "  Enthalpy change min/max [J/kg]: "
-            << gMin(dhdt_) << " - " << gMax(dhdt_) << endl;
+            << limits.min() << " - " << limits.max() << endl;
     }
 
     zeroGradientFvPatchScalarField::updateCoeffs();
diff --git a/src/thermophysicalModels/reactionThermo/derivedFvPatchFields/speciesSorption/speciesSorptionFvPatchScalarField.C b/src/thermophysicalModels/reactionThermo/derivedFvPatchFields/speciesSorption/speciesSorptionFvPatchScalarField.C
index 26a4a27fd6..2318c26550 100644
--- a/src/thermophysicalModels/reactionThermo/derivedFvPatchFields/speciesSorption/speciesSorptionFvPatchScalarField.C
+++ b/src/thermophysicalModels/reactionThermo/derivedFvPatchFields/speciesSorption/speciesSorptionFvPatchScalarField.C
@@ -321,8 +321,10 @@ patchSource() const
 
     if (debug)
     {
+        auto limits = gMinMax(dfldp);
+
         Info<< " Patch mass rate min/max [kg/m3/sec]: "
-            << gMin(dfldp) << " - " << gMax(dfldp) << endl;
+            << limits.min() << " - " << limits.max() << endl;
     }
 
     return tmp<scalarField>::New(dfldp);
@@ -391,8 +393,10 @@ void Foam::speciesSorptionFvPatchScalarField::updateCoeffs()
 
     if (debug)
     {
+        auto limits = gMinMax(dfldp_);
+
         Info<< "  Absorption rate min/max [mol/kg/sec]: "
-            << gMin(dfldp_) << " - " << gMax(dfldp_) << endl;
+            << limits.min() << " - " << limits.max() << endl;
     }
 
     zeroGradientFvPatchScalarField::updateCoeffs();
diff --git a/src/topoChangerFvMesh/movingConeTopoFvMesh/movingConeTopoFvMesh.C b/src/topoChangerFvMesh/movingConeTopoFvMesh/movingConeTopoFvMesh.C
index 0b0d1e8bef..be8b4aab60 100644
--- a/src/topoChangerFvMesh/movingConeTopoFvMesh/movingConeTopoFvMesh.C
+++ b/src/topoChangerFvMesh/movingConeTopoFvMesh/movingConeTopoFvMesh.C
@@ -134,7 +134,7 @@ void Foam::movingConeTopoFvMesh::addZonesAndModifiers()
          && fc[facei].x() < -0.003499
         )
         {
-            if ((fa[facei] & vector(1, 0, 0)) < 0)
+            if (fa[facei].x() < 0)
             {
                 flipZone1[nZoneFaces1] = true;
             }
@@ -152,7 +152,7 @@ void Foam::movingConeTopoFvMesh::addZonesAndModifiers()
         {
             zone2[nZoneFaces2] = facei;
 
-            if ((fa[facei] & vector(1, 0, 0)) > 0)
+            if (fa[facei].x() > 0)
             {
                 flipZone2[nZoneFaces2] = true;
             }
diff --git a/src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvection.C b/src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvection.C
index 6ecf2ad017..d34604b256 100644
--- a/src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvection.C
+++ b/src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvection.C
@@ -160,11 +160,9 @@ Foam::isoAdvection::isoAdvection
 
         if (porosityPtr_)
         {
-            if
-            (
-                gMin(porosityPtr_->primitiveField()) <= 0
-             || gMax(porosityPtr_->primitiveField()) > 1 + SMALL
-            )
+            auto limits = gMinMax(porosityPtr_->primitiveField());
+
+            if (limits.min() <= 0 || limits.max() > 1 + SMALL)
             {
                 FatalErrorInFunction
                     << "Porosity field has values <= 0 or > 1"
diff --git a/src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvectionTemplates.C b/src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvectionTemplates.C
index 422d86c896..b0dc0607d7 100644
--- a/src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvectionTemplates.C
+++ b/src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvectionTemplates.C
@@ -121,9 +121,11 @@ void Foam::isoAdvection::limitFluxes
     addProfilingInFunction(geometricVoF);
     DebugInFunction << endl;
 
+    auto alpha1Limits = gMinMax(alpha1_);
+
     const scalar aTol = 100*SMALL;          // Note: tolerances
-    scalar maxAlphaMinus1 = gMax(alpha1_) - 1;      // max(alphaNew - 1);
-    scalar minAlpha = gMin(alpha1_);           // min(alphaNew);
+    scalar maxAlphaMinus1 = alpha1Limits.max() - 1;  // max(alphaNew - 1);
+    scalar minAlpha = alpha1Limits.min();            // min(alphaNew);
     const label nOvershoots = 20;         // sum(pos0(alphaNew - 1 - aTol));
 
     const labelList& owner = mesh_.faceOwner();
@@ -232,8 +234,10 @@ void Foam::isoAdvection::limitFluxes
             break;
         }
 
-        maxAlphaMinus1 = gMax(alpha1_) - 1;     // max(alphaNew - 1);
-        minAlpha = gMin(alpha1_);               // min(alphaNew);
+        alpha1Limits = gMinMax(alpha1_);
+
+        maxAlphaMinus1 = alpha1Limits.max() - 1;  // max(alphaNew - 1);
+        minAlpha = alpha1Limits.min();            // min(alphaNew);
 
         if (debug)
         {
@@ -486,11 +490,13 @@ void Foam::isoAdvection::advect(const SpType& Sp, const SuType& Su)
     // Adjust dVf for unbounded cells
     limitFluxes(Sp, Su);
 
-    scalar maxAlphaMinus1 = gMax(alpha1In_) - 1;
-    scalar minAlpha = gMin(alpha1In_);
+    {
+        auto limits = gMinMax(alpha1In_);
 
-    Info<< "isoAdvection: After conservative bounding: min(alpha) = "
-        << minAlpha << ", max(alpha) = 1 + " << maxAlphaMinus1 << endl;
+        Info<< "isoAdvection: After conservative bounding:"
+            << " min(alpha) = " << limits.min()
+            << ", max(alpha) = 1 + " << (limits.max()-1) << endl;
+    }
 
     // Apply non-conservative bounding mechanisms (clipping and snapping)
     // Note: We should be able to write out alpha before this is done!