diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/laserDTRM/laserDTRM.C b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/laserDTRM/laserDTRM.C
index 40c906b57e..a106c8e5aa 100644
--- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/laserDTRM/laserDTRM.C
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/laserDTRM/laserDTRM.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2023 OpenCFD Ltd.
+    Copyright (C) 2017-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -679,7 +679,7 @@ void Foam::radiation::laserDTRM::calculate()
             }
         }
 
-        scalar totalQ = gSum(Q_.primitiveFieldRef()*mesh_.V());
+        scalar totalQ = gWeightedSum(mesh_.V(), Q_.primitiveField());
         Info << "Total energy absorbed [W]: " << totalQ << endl;
 
         if (mesh_.time().writeTime())
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/cylindricalInletVelocity/cylindricalInletVelocityFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/cylindricalInletVelocity/cylindricalInletVelocityFvPatchVectorField.C
index cc824309b6..5a33ae09a9 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/cylindricalInletVelocity/cylindricalInletVelocityFvPatchVectorField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/cylindricalInletVelocity/cylindricalInletVelocityFvPatchVectorField.C
@@ -79,7 +79,7 @@ cylindricalInletVelocityFvPatchVectorField
 :
     fixedValueFvPatchField<vector>(p, iF, dict),
     origin_(dict.getCompat<vector>("origin", {{"centre", 1712}})),
-    axis_(dict.lookup("axis")),
+    axis_(dict.get<vector>("axis")),
     axialVelocity_(Function1<scalar>::New("axialVelocity", dict, &db())),
     radialVelocity_
     (
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchField.C
index d6d8a6eb24..6e40b3b97c 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchField.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2016 OpenFOAM Foundation
-    Copyright (C) 2017-2021 OpenCFD Ltd.
+    Copyright (C) 2017-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -109,8 +109,7 @@ void Foam::fixedMeanFvPatchField<Type>::updateCoeffs()
     Field<Type> newValues(this->patchInternalField());
 
     Type meanValuePsi =
-        gSum(this->patch().magSf()*newValues)
-       /gSum(this->patch().magSf());
+        gWeightedAverage(this->patch().magSf(), newValues);
 
     if (mag(meanValue) > SMALL && mag(meanValuePsi) > 0.5*mag(meanValue))
     {
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fixedMeanOutletInlet/fixedMeanOutletInletFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/fixedMeanOutletInlet/fixedMeanOutletInletFvPatchField.C
index ff40e480f1..3e12360249 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/fixedMeanOutletInlet/fixedMeanOutletInletFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fixedMeanOutletInlet/fixedMeanOutletInletFvPatchField.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2018 OpenFOAM Foundation
-    Copyright (C) 2020-2021 OpenCFD Ltd.
+    Copyright (C) 2020-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -117,9 +117,7 @@ void Foam::fixedMeanOutletInletFvPatchField<Type>::updateCoeffs()
 
     Field<Type> newValues(this->patchInternalField());
 
-    Type meanValuePsi =
-        gSum(this->patch().magSf()*newValues)
-       /gSum(this->patch().magSf());
+    Type meanValuePsi = gWeightedAverage(this->patch().magSf(), newValues);
 
     if (mag(meanValue) > SMALL && mag(meanValuePsi) > 0.5*mag(meanValue))
     {
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/mappedField/Sampled/Sampled.C b/src/finiteVolume/fields/fvPatchFields/derived/mappedField/Sampled/Sampled.C
index 7cc3c3048d..56423e5480 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/mappedField/Sampled/Sampled.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/mappedField/Sampled/Sampled.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2022 OpenCFD Ltd.
+    Copyright (C) 2018-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -285,7 +285,7 @@ Foam::PatchFunction1Types::Sampled<Type>::value
             (
                 mag(this->mappedPatchBase::patch_.faceAreas())
             );
-            averagePsi = gSum(magSf*newValues)/gSum(magSf);
+            averagePsi = gWeightedAverage(magSf, newValues);
         }
         else
         {
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/mappedField/mappedPatchFieldBase.C b/src/finiteVolume/fields/fvPatchFields/derived/mappedField/mappedPatchFieldBase.C
index 786ab2511b..16573052aa 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/mappedField/mappedPatchFieldBase.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/mappedField/mappedPatchFieldBase.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2016 OpenFOAM Foundation
-    Copyright (C) 2018-2022 OpenCFD Ltd.
+    Copyright (C) 2018-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -788,8 +788,7 @@ Foam::mappedPatchFieldBase<Type>::mappedField
     if (setAverage_)
     {
         Type averagePsi =
-            gSum(patchField_.patch().magSf()*newValues)
-           /gSum(patchField_.patch().magSf());
+            gWeightedAverage(patchField_.patch().magSf(), newValues);
 
         if (mag(averagePsi) > 0.5*mag(average_))
         {
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/outletMappedUniformInlet/outletMappedUniformInletFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/outletMappedUniformInlet/outletMappedUniformInletFvPatchField.C
index 72ae806e5a..bbb48fa41f 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/outletMappedUniformInlet/outletMappedUniformInletFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/outletMappedUniformInlet/outletMappedUniformInletFvPatchField.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2020-2022 OpenCFD Ltd.
+    Copyright (C) 2020-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -369,7 +369,7 @@ void Foam::outletMappedUniformInletFvPatchField<Type>::updateCoeffs()
 
                 mapField.append
                 (
-                    gSum(outlet.magSf()*outletFld)/gSum(outlet.magSf())
+                    gWeightedAverage(outlet.magSf(), outletFld)
                 );
             }
         }
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/outletPhaseMeanVelocity/outletPhaseMeanVelocityFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/outletPhaseMeanVelocity/outletPhaseMeanVelocityFvPatchVectorField.C
index 6d90bf382b..01eee817ae 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/outletPhaseMeanVelocity/outletPhaseMeanVelocityFvPatchVectorField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/outletPhaseMeanVelocity/outletPhaseMeanVelocityFvPatchVectorField.C
@@ -130,8 +130,8 @@ void Foam::outletPhaseMeanVelocityFvPatchVectorField::updateCoeffs()
         patch().lookupPatchField<volScalarField>(alphaName_)
     );
 
-    alphap = max(alphap, scalar(0));
-    alphap = min(alphap, scalar(1));
+    // Clamp to 0-1 range
+    alphap.clamp_range(Foam::zero_one{});
 
     // Get the patchInternalField (zero-gradient field)
     vectorField Uzg(patchInternalField());
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/plenumPressure/plenumPressureFvPatchScalarField.C b/src/finiteVolume/fields/fvPatchFields/derived/plenumPressure/plenumPressureFvPatchScalarField.C
index 4b5dc3559f..2dc83bd329 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/plenumPressure/plenumPressureFvPatchScalarField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/plenumPressure/plenumPressureFvPatchScalarField.C
@@ -211,7 +211,7 @@ void Foam::plenumPressureFvPatchScalarField::updateCoeffs()
     {
         if (hasRho_)
         {
-            massFlowRate = - gSum(rho_*phi);
+            massFlowRate = - rho_*gSum(phi);
         }
         else
         {
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/pressurePIDControlInletVelocity/pressurePIDControlInletVelocityFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/pressurePIDControlInletVelocity/pressurePIDControlInletVelocityFvPatchVectorField.C
index 84c265f990..4ca7a6c842 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/pressurePIDControlInletVelocity/pressurePIDControlInletVelocityFvPatchVectorField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/pressurePIDControlInletVelocity/pressurePIDControlInletVelocityFvPatchVectorField.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2015-2020 OpenCFD Ltd.
+    Copyright (C) 2015-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -243,7 +243,7 @@ void Foam::pressurePIDControlInletVelocityFvPatchVectorField::updateCoeffs()
     {
         const auto& rhop = patch().lookupPatchField<volScalarField>(rhoName_);
 
-        rho = gSum(rhop*patch().magSf())/gSum(patch().magSf());
+        rho = gWeightedAverage(patch().magSf(), rhop);
     }
     else
     {
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/rotatingWallVelocity/rotatingWallVelocityFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/rotatingWallVelocity/rotatingWallVelocityFvPatchVectorField.C
index 095c96367f..ca57fe5ad7 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/rotatingWallVelocity/rotatingWallVelocityFvPatchVectorField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/rotatingWallVelocity/rotatingWallVelocityFvPatchVectorField.C
@@ -41,7 +41,7 @@ rotatingWallVelocityFvPatchVectorField
 )
 :
     fixedValueFvPatchField<vector>(p, iF),
-    origin_(),
+    origin_(Zero),
     axis_(Zero),
     omega_(nullptr)
 {}
@@ -56,8 +56,8 @@ rotatingWallVelocityFvPatchVectorField
 )
 :
     fixedValueFvPatchField<vector>(p, iF, dict, IOobjectOption::NO_READ),
-    origin_(dict.lookup("origin")),
-    axis_(dict.lookup("axis")),
+    origin_(dict.get<vector>("origin")),
+    axis_(dict.get<vector>("axis")),
     omega_(Function1<scalar>::New("omega", dict, &db()))
 {
     if (!this->readValueEntry(dict))
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/swirlFanVelocity/swirlFanVelocityFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/swirlFanVelocity/swirlFanVelocityFvPatchField.C
index 86e5291206..5e75e74855 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/swirlFanVelocity/swirlFanVelocityFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/swirlFanVelocity/swirlFanVelocityFvPatchField.C
@@ -57,8 +57,7 @@ void Foam::swirlFanVelocityFvPatchField::calcFanJump()
             deltaP /= patch().lookupPatchField<volScalarField>(rhoName_);
         }
 
-        const vector axisHat =
-            gSum(patch().nf()*patch().magSf())/gSum(patch().magSf());
+        const vector axisHat = gSum(patch().Sf())/gSum(patch().magSf());
 
         vectorField tanDir
         (
@@ -131,7 +130,7 @@ Foam::swirlFanVelocityFvPatchField::swirlFanVelocityFvPatchField
     phiName_("phi"),
     pName_("p"),
     rhoName_("rho"),
-    origin_(),
+    origin_(Zero),
     rpm_(nullptr),
     fanEff_(1),
     rEff_(0),
@@ -158,7 +157,7 @@ Foam::swirlFanVelocityFvPatchField::swirlFanVelocityFvPatchField
         (
             "origin",
             returnReduceOr(patch().size())
-          ? gSum(patch().Cf()*patch().magSf())/gSum(patch().magSf())
+          ? gWeightedAverage(patch().magSf(), patch().Cf())
           : Zero
         )
     ),
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/swirlFlowRateInletVelocity/swirlFlowRateInletVelocityFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/swirlFlowRateInletVelocity/swirlFlowRateInletVelocityFvPatchVectorField.C
index b4e4a38435..beb80502b3 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/swirlFlowRateInletVelocity/swirlFlowRateInletVelocityFvPatchVectorField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/swirlFlowRateInletVelocity/swirlFlowRateInletVelocityFvPatchVectorField.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2018-2022 OpenCFD Ltd.
+    Copyright (C) 2018-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -45,7 +45,7 @@ swirlFlowRateInletVelocityFvPatchVectorField
     fixedValueFvPatchField<vector>(p, iF),
     phiName_("phi"),
     rhoName_("rho"),
-    origin_(),
+    origin_(Zero),
     axis_(Zero),
     flowRate_(),
     rpm_()
@@ -69,7 +69,7 @@ swirlFlowRateInletVelocityFvPatchVectorField
         (
             "origin",
             returnReduceOr(patch().size())
-          ? gSum(patch().Cf()*patch().magSf())/gSum(patch().magSf())
+          ? gWeightedAverage(patch().magSf(), patch().Cf())
           : Zero
         )
     ),
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/swirlInletVelocity/swirlInletVelocityFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/swirlInletVelocity/swirlInletVelocityFvPatchVectorField.C
index d3957f09ae..76ab66f0b1 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/swirlInletVelocity/swirlInletVelocityFvPatchVectorField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/swirlInletVelocity/swirlInletVelocityFvPatchVectorField.C
@@ -74,8 +74,8 @@ swirlInletVelocityFvPatchVectorField
 )
 :
     fixedValueFvPatchField<vector>(p, iF, dict),
-    origin_(dict.lookup("origin")),
-    axis_(dict.lookup("axis")),
+    origin_(dict.get<vector>("origin")),
+    axis_(dict.get<vector>("axis")),
     axialVelocity_(Function1<scalar>::New("axialVelocity", dict, &db())),
     radialVelocity_
     (
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C
index a59ecba23f..888e5bc6c0 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C
@@ -895,8 +895,7 @@ void Foam::turbulentDFSEMInletFvPatchVectorField::updateCoeffs()
         // (PCR:p. 522)
         const vector UBulk
         (
-            gSum(UMean()*patch().magSf())
-           /(gSum(patch().magSf()) + ROOTVSMALL)
+            gWeightedAverage(patch().magSf(), UMean())
         );
 
         // Move eddies using bulk velocity
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/turbulentDigitalFilterInletFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/turbulentDigitalFilterInletFvPatchField.C
index f947bfb24f..770dfcc135 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/turbulentDigitalFilterInletFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/turbulentDigitalFilterInletFvPatchField.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2022 OpenCFD Ltd.
+    Copyright (C) 2019-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -174,11 +174,7 @@ void Foam::turbulentDigitalFilterInletFvPatchField<Type>::mapMean
     const vectorField& mean = tmean.cref();
 
     // Calculate flow-rate correction factor for vector fields (KCX:Eq. 8)
-    const vector bulk
-    (
-        gSum(mean*this->patch().magSf())
-       /(gSum(this->patch().magSf()) + ROOTVSMALL)
-    );
+    const vector bulk = gWeightedAverage(this->patch().magSf(), mean);
 
     const scalar correct
     (
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/variableHeightFlowRateInletVelocity/variableHeightFlowRateInletVelocityFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/variableHeightFlowRateInletVelocity/variableHeightFlowRateInletVelocityFvPatchVectorField.C
index cc709fe383..067c93d67f 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/variableHeightFlowRateInletVelocity/variableHeightFlowRateInletVelocityFvPatchVectorField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/variableHeightFlowRateInletVelocity/variableHeightFlowRateInletVelocityFvPatchVectorField.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2016 OpenFOAM Foundation
-    Copyright (C) 2021 OpenCFD Ltd.
+    Copyright (C) 2021-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -125,14 +125,14 @@ void Foam::variableHeightFlowRateInletVelocityFvPatchVectorField::updateCoeffs()
 
     scalarField alphap(alpha.boundaryField()[patch().index()]);
 
-    alphap = max(alphap, scalar(0));
-    alphap = min(alphap, scalar(1));
+    // Clamp to 0-1 range
+    alphap.clamp_range(Foam::zero_one{});
 
     const scalar t = db().time().timeOutputValue();
     scalar flowRate = flowRate_->value(t);
 
     // a simpler way of doing this would be nice
-    scalar avgU = -flowRate/gSum(patch().magSf()*alphap);
+    const scalar avgU = -flowRate/gWeightedSum(patch().magSf(), alphap);
 
     vectorField n(patch().nf());
 
diff --git a/src/functionObjects/field/DESModelRegions/DESModelRegions.C b/src/functionObjects/field/DESModelRegions/DESModelRegions.C
index ce2e343942..d026e4530b 100644
--- a/src/functionObjects/field/DESModelRegions/DESModelRegions.C
+++ b/src/functionObjects/field/DESModelRegions/DESModelRegions.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2015 OpenFOAM Foundation
-    Copyright (C) 2015-2020 OpenCFD Ltd.
+    Copyright (C) 2015-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -110,23 +110,18 @@ bool Foam::functionObjects::DESModelRegions::execute()
 {
     Log << type() << " " << name() <<  " execute:" << nl;
 
-    volScalarField& DESModelRegions =
-        lookupObjectRef<volScalarField>(resultName_);
+    auto& regions = lookupObjectRef<volScalarField>(resultName_);
 
-
-    if (foundObject<DESModelBase>(turbulenceModel::propertiesName))
+    if
+    (
+        const auto* modelp
+      = cfindObject<DESModelBase>(turbulenceModel::propertiesName)
+    )
     {
-        const DESModelBase& model =
-            lookupObject<DESModelBase>
-            (
-                turbulenceModel::propertiesName
-            );
-
-        DESModelRegions == model.LESRegion();
+        regions == (*modelp).LESRegion();
 
         const scalar prc =
-            gSum(DESModelRegions.primitiveField()*mesh_.V())
-           /gSum(mesh_.V())*100.0;
+            gWeightedAverage(mesh_.V(), regions.primitiveField())*100.0;
 
         file() << time_.value()
             << token::TAB << prc
@@ -149,14 +144,13 @@ bool Foam::functionObjects::DESModelRegions::execute()
 
 bool Foam::functionObjects::DESModelRegions::write()
 {
-    const volScalarField& DESModelRegions =
-        lookupObject<volScalarField>(resultName_);
+    const auto& regions = lookupObject<volScalarField>(resultName_);
 
     Log << type() << " " << name() <<  " output:" << nl
-        << "    writing field " << DESModelRegions.name() << nl
+        << "    writing field " << regions.name() << nl
         << endl;
 
-    DESModelRegions.write();
+    regions.write();
 
     return true;
 }
diff --git a/src/functionObjects/field/comfort/comfort.C b/src/functionObjects/field/comfort/comfort.C
index ed0d465e5c..9eca338fe8 100644
--- a/src/functionObjects/field/comfort/comfort.C
+++ b/src/functionObjects/field/comfort/comfort.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2019 OpenFOAM Foundation
-    Copyright (C) 2021-2023 OpenCFD Ltd.
+    Copyright (C) 2021-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -89,12 +89,15 @@ Foam::dimensionedScalar Foam::functionObjects::comfort::Trad() const
 
             if (isType<wallFvPatch>(pTBf))
             {
-                areaIntegral += gSum(pSf);
-                TareaIntegral += gSum(pSf*pT);
+                areaIntegral += sum(pSf);
+                TareaIntegral += weightedSum(pSf, pT);
             }
         }
 
-        Trad.value() = TareaIntegral/areaIntegral;
+        reduce(TareaIntegral, sumOp<scalar>());
+        reduce(areaIntegral, sumOp<scalar>());
+
+        Trad.value() = TareaIntegral/(areaIntegral + VSMALL);
     }
 
     // Bounds based on EN ISO 7730
diff --git a/src/functionObjects/field/wallHeatFlux/wallHeatFlux.C b/src/functionObjects/field/wallHeatFlux/wallHeatFlux.C
index 70618231d0..4112083a3e 100644
--- a/src/functionObjects/field/wallHeatFlux/wallHeatFlux.C
+++ b/src/functionObjects/field/wallHeatFlux/wallHeatFlux.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 OpenFOAM Foundation
-    Copyright (C) 2016-2024 OpenCFD Ltd.
+    Copyright (C) 2016-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -272,7 +272,7 @@ bool Foam::functionObjects::wallHeatFlux::execute()
         const scalarField& hfp = wallHeatFlux.boundaryField()[patchi];
 
         const MinMax<scalar> limits = gMinMax(hfp);
-        const scalar integralHfp = gSum(magSf[patchi]*hfp);
+        const scalar integralHfp = gWeightedSum(magSf[patchi], hfp);
 
         if (Pstream::master())
         {
diff --git a/src/fvMotionSolver/pointPatchFields/derived/angularOscillatingDisplacement/angularOscillatingDisplacementPointPatchVectorField.C b/src/fvMotionSolver/pointPatchFields/derived/angularOscillatingDisplacement/angularOscillatingDisplacementPointPatchVectorField.C
index 7ff26948cb..0594512e12 100644
--- a/src/fvMotionSolver/pointPatchFields/derived/angularOscillatingDisplacement/angularOscillatingDisplacementPointPatchVectorField.C
+++ b/src/fvMotionSolver/pointPatchFields/derived/angularOscillatingDisplacement/angularOscillatingDisplacementPointPatchVectorField.C
@@ -64,8 +64,8 @@ angularOscillatingDisplacementPointPatchVectorField
 )
 :
     fixedValuePointPatchField<vector>(p, iF, dict),
-    axis_(dict.lookup("axis")),
-    origin_(dict.lookup("origin")),
+    axis_(dict.get<vector>("axis")),
+    origin_(dict.get<vector>("origin")),
     angle0_(dict.get<scalar>("angle0")),
     amplitude_(dict.get<scalar>("amplitude")),
     omega_(dict.get<scalar>("omega"))
diff --git a/src/fvMotionSolver/pointPatchFields/derived/angularOscillatingVelocity/angularOscillatingVelocityPointPatchVectorField.C b/src/fvMotionSolver/pointPatchFields/derived/angularOscillatingVelocity/angularOscillatingVelocityPointPatchVectorField.C
index 792daffdcf..ff3ca22233 100644
--- a/src/fvMotionSolver/pointPatchFields/derived/angularOscillatingVelocity/angularOscillatingVelocityPointPatchVectorField.C
+++ b/src/fvMotionSolver/pointPatchFields/derived/angularOscillatingVelocity/angularOscillatingVelocityPointPatchVectorField.C
@@ -64,8 +64,8 @@ angularOscillatingVelocityPointPatchVectorField
 )
 :
     fixedValuePointPatchField<vector>(p, iF, dict),
-    axis_(dict.lookup("axis")),
-    origin_(dict.lookup("origin")),
+    axis_(dict.get<vector>("axis")),
+    origin_(dict.get<vector>("origin")),
     angle0_(dict.get<scalar>("angle0")),
     amplitude_(dict.get<scalar>("amplitude")),
     omega_(dict.get<scalar>("omega"))
diff --git a/src/fvOptions/sources/derived/directionalPressureGradientExplicitSource/directionalPressureGradientExplicitSource.C b/src/fvOptions/sources/derived/directionalPressureGradientExplicitSource/directionalPressureGradientExplicitSource.C
index c2d6ea05c1..7a78483534 100644
--- a/src/fvOptions/sources/derived/directionalPressureGradientExplicitSource/directionalPressureGradientExplicitSource.C
+++ b/src/fvOptions/sources/derived/directionalPressureGradientExplicitSource/directionalPressureGradientExplicitSource.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2015-2023 OpenCFD Ltd.
+    Copyright (C) 2015-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -337,7 +337,7 @@ void Foam::fv::directionalPressureGradientExplicitSource::correct
                     );
                 const scalarField rho(turbModel.rho(),cells_);
                 const scalarField cv(mesh_.V(), cells_);
-                scalar rhoAve = gSumProd(rho, cv)/gSum(cv);
+                scalar rhoAve = gWeightedAverage(cv, rho);
                 volFlowRate = mag(totalphi)/rhoAve;
             }
 
diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C
index 6381b3733c..ddb5f7fd0d 100644
--- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C
@@ -2687,7 +2687,7 @@ Foam::label Foam::snappyRefineDriver::directionalSmooth
 
                     const vectorField residual(baseNewPoints - baseMeshPoints);
                     {
-                        scalar res(gSum(mag(residual)));
+                        scalar res(gSumMag(residual));
 
                         if (iter == 0)
                         {
diff --git a/src/meshTools/PatchFunction1/MappedFile/MappedFile.C b/src/meshTools/PatchFunction1/MappedFile/MappedFile.C
index 16d7fbc850..ef0856fe7b 100644
--- a/src/meshTools/PatchFunction1/MappedFile/MappedFile.C
+++ b/src/meshTools/PatchFunction1/MappedFile/MappedFile.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2022 OpenCFD Ltd.
+    Copyright (C) 2018-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -700,8 +700,7 @@ Foam::PatchFunction1Types::MappedFile<Type>::value
 
         if (this->faceValues())
         {
-            const scalarField magSf(mag(this->patch_.faceAreas()));
-            averagePsi = gSum(magSf*fld)/gSum(magSf);
+            averagePsi = gWeightedAverage(this->patch_.magFaceAreas(), fld);
         }
         else
         {
diff --git a/src/optimisation/adjointOptimisation/adjoint/adjointBoundaryConditions/adjointRotatingWallVelocity/adjointRotatingWallVelocityFvPatchVectorField.C b/src/optimisation/adjointOptimisation/adjoint/adjointBoundaryConditions/adjointRotatingWallVelocity/adjointRotatingWallVelocityFvPatchVectorField.C
index ad2e887809..c6f0996b86 100644
--- a/src/optimisation/adjointOptimisation/adjoint/adjointBoundaryConditions/adjointRotatingWallVelocity/adjointRotatingWallVelocityFvPatchVectorField.C
+++ b/src/optimisation/adjointOptimisation/adjoint/adjointBoundaryConditions/adjointRotatingWallVelocity/adjointRotatingWallVelocityFvPatchVectorField.C
@@ -39,7 +39,7 @@ adjointRotatingWallVelocityFvPatchVectorField
 )
 :
     adjointWallVelocityFvPatchVectorField(p, iF),
-    origin_(),
+    origin_(Zero),
     axis_(Zero),
     omega_(nullptr)
 {}
diff --git a/src/optimisation/adjointOptimisation/adjoint/objectives/geometric/objectiveTopOSolidVolume/objectiveTopOSolidVolume.C b/src/optimisation/adjointOptimisation/adjoint/objectives/geometric/objectiveTopOSolidVolume/objectiveTopOSolidVolume.C
index b315d85985..7324bfc065 100644
--- a/src/optimisation/adjointOptimisation/adjoint/objectives/geometric/objectiveTopOSolidVolume/objectiveTopOSolidVolume.C
+++ b/src/optimisation/adjointOptimisation/adjoint/objectives/geometric/objectiveTopOSolidVolume/objectiveTopOSolidVolume.C
@@ -71,17 +71,17 @@ objectiveTopOSolidVolume::objectiveTopOSolidVolume
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-
 scalar objectiveTopOSolidVolume::J()
 {
     J_ = Zero;
-    if (mesh_.foundObject<volScalarField>("beta"))
+
+    if (const auto* betap = mesh_.cfindObject<volScalarField>("beta"))
     {
-        const volScalarField& beta = mesh_.lookupObject<volScalarField>("beta");
-        const DimensionedField<scalar, volMesh>& V = mesh_.V();
+        const auto& beta = *betap;
         const scalar time = mesh_.time().timeOutputValue();
+
         J_ =
-            gSum(beta.primitiveField()*V)/gSum(V)
+            gWeightedAverage(mesh_.V(), beta.primitiveField())
           - targetPercentage_->value(time);
         if (percentInDenom_)
         {
@@ -101,6 +101,7 @@ scalar objectiveTopOSolidVolume::J()
 void objectiveTopOSolidVolume::update_dJdb()
 {
     const scalar time = mesh_.time().timeOutputValue();
+
     dJdbPtr_().primitiveFieldRef() = scalar(1)/gSum(mesh_.V());
     if (percentInDenom_)
     {
diff --git a/src/optimisation/adjointOptimisation/adjoint/objectives/geometric/objectiveTopOVolume/objectiveTopOVolume.C b/src/optimisation/adjointOptimisation/adjoint/objectives/geometric/objectiveTopOVolume/objectiveTopOVolume.C
index d15b68e3df..9d9b4f5da8 100644
--- a/src/optimisation/adjointOptimisation/adjoint/objectives/geometric/objectiveTopOVolume/objectiveTopOVolume.C
+++ b/src/optimisation/adjointOptimisation/adjoint/objectives/geometric/objectiveTopOVolume/objectiveTopOVolume.C
@@ -71,17 +71,17 @@ objectiveTopOVolume::objectiveTopOVolume
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-
 scalar objectiveTopOVolume::J()
 {
     J_ = Zero;
-    if (mesh_.foundObject<volScalarField>("beta"))
+
+    if (const auto* betap = mesh_.cfindObject<volScalarField>("beta"))
     {
-        const volScalarField& beta = mesh_.lookupObject<volScalarField>("beta");
-        const DimensionedField<scalar, volMesh>& V = mesh_.V();
+        const auto& beta = *betap;
         const scalar time = mesh_.time().timeOutputValue();
+
         J_ =
-            scalar(1) - gSum(beta.primitiveField()*V)/gSum(V)
+            scalar(1) - gWeightedAverage(mesh_.V(), beta.primitiveField())
           - targetPercentage_->value(time);
         if (percentInDenom_)
         {
@@ -101,6 +101,7 @@ scalar objectiveTopOVolume::J()
 void objectiveTopOVolume::update_dJdb()
 {
     const scalar time = mesh_.time().timeOutputValue();
+
     dJdbPtr_().primitiveFieldRef() = -scalar(1)/gSum(mesh_.V());
     if (percentInDenom_)
     {
diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/betaMax/betaMax/betaMax.C b/src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/betaMax/betaMax/betaMax.C
index 99a1c78eb3..d26f95171d 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/betaMax/betaMax/betaMax.C
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/betaMax/betaMax/betaMax.C
@@ -62,8 +62,9 @@ Foam::scalar Foam::betaMax::computeLength(const dictionary& dict) const
             for (const label pI : inletPatches)
             {
                 const fvPatch& patch = mesh_.boundary()[pI];
-                length += gSum(patch.magSf());
+                length += sum(patch.magSf());
             }
+            reduce(length, sumOp<scalar>());
 
             // Divide with the span in the empty direction
             const Vector<label>& geometricD = mesh_.geometricD();
diff --git a/src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C b/src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C
index 7adb2c0d55..56cdee16f8 100644
--- a/src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C
+++ b/src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C
@@ -513,7 +513,7 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                  Info<< "  qEffVap: " << gMin(qEff) << " - "
                      << gMax(qEff) << endl;
 
-                scalar Qeff = gSum(qEff*patch().magSf());
+                scalar Qeff = gWeightedSum(patch().magSf(), qEff);
                 Info<< " Effective heat transfer rate to vapor:" << Qeff
                     << nl << endl;
             }
@@ -1111,7 +1111,7 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                 Info<< "  alphatlEff: " << gMin(liquidw*(*this + alphaw))
                     << " - " << gMax(liquidw*(*this + alphaw)) << endl;
 
-                scalar Qeff = gSum(qEff*patch().magSf());
+                scalar Qeff = gWeightedSum(patch().magSf(), qEff);
                 Info<< " Effective heat transfer rate to liquid: " << Qeff
                     << endl << nl;
 
@@ -1166,7 +1166,7 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                         *hew.snGrad()
                     );
 
-                    scalar Qc = gSum(qc*patch().magSf());
+                    scalar Qc = gWeightedSum(patch().magSf(), qc);
                     Info<< " Convective heat transfer: " << Qc << endl;
 
                     const scalarField qFilm
@@ -1174,7 +1174,7 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                         relax*fLiquid*nFilms*htcFilmBoiling*(Tw - Tsatw)
                     );
 
-                    scalar QFilm = gSum(qFilm*patch().magSf());
+                    scalar QFilm = gWeightedSum(patch().magSf(), qFilm);
                     Info<< " Film boiling heat transfer: " << QFilm << endl;
 
                     Info<< " Htc Film Boiling coeff: "
@@ -1199,7 +1199,7 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                         )
                     );
 
-                    scalar QsubCool = gSum(qSubCool*patch().magSf());
+                    scalar QsubCool = gWeightedSum(patch().magSf(), qSubCool);
 
                     Info<< " Sub Cool boiling heat transfer: " << QsubCool
                         << endl;
diff --git a/src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/fixedMultiPhaseHeatFlux/fixedMultiPhaseHeatFluxFvPatchScalarField.C b/src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/fixedMultiPhaseHeatFlux/fixedMultiPhaseHeatFluxFvPatchScalarField.C
index cce56b32a4..3f32741033 100644
--- a/src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/fixedMultiPhaseHeatFlux/fixedMultiPhaseHeatFluxFvPatchScalarField.C
+++ b/src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/fixedMultiPhaseHeatFlux/fixedMultiPhaseHeatFluxFvPatchScalarField.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015-2020 OpenFOAM Foundation
-    Copyright (C) 2020-2021 OpenCFD Ltd.
+    Copyright (C) 2020-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -157,9 +157,11 @@ void Foam::fixedMultiPhaseHeatFluxFvPatchScalarField::updateCoeffs()
 
     if (debug)
     {
+        MinMax<scalar> limits = gMinMax(Q);
+
         Info<< patch().name() << " " << ": overall heat flux "
-            << gMin(Q) << " - " << gMax(Q) << " W/m2, power: "
-            << gSum(patch().magSf()*Q) << " W" << endl;
+            << limits.min() << " - " << limits.max() << " W/m2, power: "
+            << gWeightedSum(patch().magSf(), Q) << " W" << endl;
     }
 
     operator==((scalar(1) - relax_)*Tp + relax_*max(Tmin_,(q_ + A)/(B)));
diff --git a/src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisRadiativeCoupledMixed/filmPyrolysisRadiativeCoupledMixedFvPatchScalarField.C b/src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisRadiativeCoupledMixed/filmPyrolysisRadiativeCoupledMixedFvPatchScalarField.C
index 6bc8bbed29..7a7499e571 100644
--- a/src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisRadiativeCoupledMixed/filmPyrolysisRadiativeCoupledMixedFvPatchScalarField.C
+++ b/src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisRadiativeCoupledMixed/filmPyrolysisRadiativeCoupledMixedFvPatchScalarField.C
@@ -402,10 +402,12 @@ void filmPyrolysisRadiativeCoupledMixedFvPatchScalarField::updateCoeffs()
 
     if (debug)
     {
-        scalar Qc = gSum(qConv*patch().magSf());
-        scalar qr = gSum(qRad*patch().magSf());
+        scalar Qc = gWeightedSum(patch().magSf(), qConv);
+        scalar qr = gWeightedSum(patch().magSf(), qRad);
         scalar Qt = gSum((qConv + qRad)*patch().magSf());
 
+        MinMax<scalar> limits = gMinMax(*this);
+
         Info<< mesh.name() << ':'
             << patch().name() << ':'
             << this->internalField().name() << " <- "
@@ -416,8 +418,8 @@ void filmPyrolysisRadiativeCoupledMixedFvPatchScalarField::updateCoeffs()
             << "     radiative heat [W] : " << qr << nl
             << "     total heat     [W] : " << Qt << nl
             << "     wall temperature "
-            << " min:" << gMin(*this)
-            << " max:" << gMax(*this)
+            << " min:" << limits.min()
+            << " max:" << limits.max()
             << " avg:" << gAverage(*this)
             << endl;
     }
diff --git a/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C b/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C
index a35106bc1f..4e3ca4b73e 100644
--- a/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C
+++ b/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2019-2023 OpenCFD Ltd.
+    Copyright (C) 2019-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -1048,15 +1048,20 @@ void kinematicSingleLayer::info()
     outputProperties().readIfPresent("addedMassTotal", addedMassTotal);
     addedMassTotal += returnReduce(addedMassTotal_, sumOp<scalar>());
 
+    MinMax<scalar> limits_U = gMinMaxMag(Uinternal);
+    MinMax<scalar> limits_delta = gMinMax(deltaInternal);
+
     Info<< indent << "added mass         = " << addedMassTotal << nl
         << indent << "current mass       = "
-        << gSum((deltaRho_*magSf())()) << nl
-        << indent << "min/max(mag(U))    = " << gMin(mag(Uinternal)) << ", "
-        << gMax(mag(Uinternal)) << nl
-        << indent << "min/max(delta)     = " << gMin(deltaInternal) << ", "
-        << gMax(deltaInternal) << nl
+        << gWeightedSum(magSf(), deltaRho_) << nl
+        << indent << "min/max(mag(U))    = "
+        << limits_U.min() << ", "
+        << limits_U.max() << nl
+        << indent << "min/max(delta)     = "
+        << limits_delta.min() << ", "
+        << limits_delta.max() << nl
         << indent << "coverage           = "
-        << gSum(alpha_.primitiveField()*magSf())/gSum(magSf()) <<  nl;
+        << gWeightedAverage(magSf(), alpha_.primitiveField()) <<  nl;
 
     injection_.info(Info);
     transfer_.info(Info);
diff --git a/src/thermoTools/derivedFvPatchFields/fixedIncidentRadiation/fixedIncidentRadiationFvPatchScalarField.C b/src/thermoTools/derivedFvPatchFields/fixedIncidentRadiation/fixedIncidentRadiationFvPatchScalarField.C
index d005c768df..04e1438bc8 100644
--- a/src/thermoTools/derivedFvPatchFields/fixedIncidentRadiation/fixedIncidentRadiationFvPatchScalarField.C
+++ b/src/thermoTools/derivedFvPatchFields/fixedIncidentRadiation/fixedIncidentRadiationFvPatchScalarField.C
@@ -170,14 +170,16 @@ void Foam::radiation::fixedIncidentRadiationFvPatchScalarField::updateCoeffs()
 
     if (debug)
     {
+        MinMax<scalar> limits = gMinMax(*this);
+
         scalar qr = gSum(kappa(*this)*gradient()*patch().magSf());
         Info<< patch().boundaryMesh().mesh().name() << ':'
             << patch().name() << ':'
             << this->internalField().name() << " -> "
             << " radiativeFlux:" << qr
             << " walltemperature "
-            << " min:" << gMin(*this)
-            << " max:" << gMax(*this)
+            << " min:" << limits.min()
+            << " max:" << limits.max()
             << " avg:" << gAverage(*this)
             << endl;
     }
diff --git a/src/thermoTools/derivedFvPatchFields/humidityTemperatureCoupledMixed/humidityTemperatureCoupledMixedFvPatchScalarField.C b/src/thermoTools/derivedFvPatchFields/humidityTemperatureCoupledMixed/humidityTemperatureCoupledMixedFvPatchScalarField.C
index d59610e575..abc00d8c4d 100644
--- a/src/thermoTools/derivedFvPatchFields/humidityTemperatureCoupledMixed/humidityTemperatureCoupledMixedFvPatchScalarField.C
+++ b/src/thermoTools/derivedFvPatchFields/humidityTemperatureCoupledMixed/humidityTemperatureCoupledMixedFvPatchScalarField.C
@@ -733,11 +733,13 @@ void Foam::humidityTemperatureCoupledMixedFvPatchScalarField::updateCoeffs()
 
     if (debug && fluid_)
     {
-        scalar Qdm = gSum(dm*magSf);
+        scalar Qdm = gWeightedSum(magSf, dm);
         scalar QMass = gSum(mass_);
         scalar Qt = gSum(myKDelta_*(Tp - Tin)*magSf);
         scalar QtSolid = gSum(KDeltaNbr*(Tp - nbrIntFld)*magSf);
 
+        MinMax<scalar> limits = gMinMax(*this);
+
         Info<< mesh.name() << ':'
             << patch().name() << ':'
             << internalField().name() << " <- "
@@ -751,8 +753,8 @@ void Foam::humidityTemperatureCoupledMixedFvPatchScalarField::updateCoeffs()
             << "     Total heat (>0 leaving the wall to the solid) [W] : "
             << QtSolid << nl
             << "     wall temperature "
-            << " min:" << gMin(*this)
-            << " max:" << gMax(*this)
+            << " min:" << limits.min()
+            << " max:" << limits.max()
             << " avg:" << gAverage(*this)
             << endl;
     }
diff --git a/src/thermoTools/derivedFvPatchFields/lumpedMassWallTemperature/lumpedMassWallTemperatureFvPatchScalarField.C b/src/thermoTools/derivedFvPatchFields/lumpedMassWallTemperature/lumpedMassWallTemperatureFvPatchScalarField.C
index d2ddbbe2e2..df66f3288c 100644
--- a/src/thermoTools/derivedFvPatchFields/lumpedMassWallTemperature/lumpedMassWallTemperatureFvPatchScalarField.C
+++ b/src/thermoTools/derivedFvPatchFields/lumpedMassWallTemperature/lumpedMassWallTemperatureFvPatchScalarField.C
@@ -168,7 +168,7 @@ void Foam::lumpedMassWallTemperatureFvPatchScalarField::updateCoeffs()
     const scalarField q(tkappa.ref()*snGrad());
 
     // Total heat in or out of the wall
-    const scalar Q = gSum(q*magSf);
+    const scalar Q = gWeightedSum(magSf, q);
 
     Tp += -(Q/mass_/Cp_)*deltaT;
 
diff --git a/src/thermoTools/derivedFvPatchFields/outletMachNumberPressure/outletMachNumberPressureFvPatchScalarField.C b/src/thermoTools/derivedFvPatchFields/outletMachNumberPressure/outletMachNumberPressureFvPatchScalarField.C
index 8c17484d03..a45417038c 100644
--- a/src/thermoTools/derivedFvPatchFields/outletMachNumberPressure/outletMachNumberPressureFvPatchScalarField.C
+++ b/src/thermoTools/derivedFvPatchFields/outletMachNumberPressure/outletMachNumberPressureFvPatchScalarField.C
@@ -212,7 +212,7 @@ void Foam::outletMachNumberPressureFvPatchScalarField::updateCoeffs()
         }
 
         const scalarField r(pBack_/ptot);
-        const scalar area = gSum(mag(patch().Sf()));
+        const scalar area = gSum(patch().magSf());
         M =
             A1_/(c1_*area)
            *sqrt(2/(gamma-1)*(pow(r, 2/gamma) - pow(r, (gamma+1)/gamma)));
diff --git a/src/thermoTools/derivedFvPatchFields/outletMappedUniformInletHeatAddition/outletMappedUniformInletHeatAdditionFvPatchField.C b/src/thermoTools/derivedFvPatchFields/outletMappedUniformInletHeatAddition/outletMappedUniformInletHeatAdditionFvPatchField.C
index 05d0f08b11..bbede04dbc 100644
--- a/src/thermoTools/derivedFvPatchFields/outletMappedUniformInletHeatAddition/outletMappedUniformInletHeatAdditionFvPatchField.C
+++ b/src/thermoTools/derivedFvPatchFields/outletMappedUniformInletHeatAddition/outletMappedUniformInletHeatAdditionFvPatchField.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2016-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -153,7 +153,7 @@ void Foam::outletMappedUniformInletHeatAdditionFvPatchField::updateCoeffs()
         db().lookupObject<surfaceScalarField>(phiName_);
 
     const scalarField& outletPatchPhi = phi.boundaryField()[outletPatchID];
-    scalar sumOutletPatchPhi = gSum(outletPatchPhi);
+    const scalar sumOutletPatchPhi = gSum(outletPatchPhi);
 
     if (sumOutletPatchPhi > SMALL)
     {
@@ -164,22 +164,21 @@ void Foam::outletMappedUniformInletHeatAdditionFvPatchField::updateCoeffs()
         const scalarField& pT = thermo.T().boundaryField()[outletPatchID];
 
         scalar averageOutletField =
-            gSum(outletPatchPhi*outletPatchField)/sumOutletPatchPhi;
+            gWeightedSum(outletPatchPhi, outletPatchField)/sumOutletPatchPhi;
 
         // Calculate Q as a function of average outlet temperature
         const scalar Q = Qptr_->value(averageOutletField);
 
         const scalarField Cpf(thermo.Cp(pp, pT, outletPatchID));
 
-        scalar totalPhiCp = gSum(outletPatchPhi)*gAverage(Cpf);
+        scalar totalPhiCp = sumOutletPatchPhi*gAverage(Cpf);
 
         operator==(clamp(averageOutletField + Q/totalPhiCp, TMin_, TMax_));
     }
     else
     {
         scalar averageOutletField =
-            gSum(outletPatch.magSf()*outletPatchField)
-           /gSum(outletPatch.magSf());
+            gWeightedAverage(outletPatch.magSf(), outletPatchField);
 
         operator==(averageOutletField);
     }
diff --git a/src/thermoTools/derivedFvPatchFields/totalFlowRateAdvectiveDiffusive/totalFlowRateAdvectiveDiffusiveFvPatchScalarField.C b/src/thermoTools/derivedFvPatchFields/totalFlowRateAdvectiveDiffusive/totalFlowRateAdvectiveDiffusiveFvPatchScalarField.C
index e65ce1728a..2e2a70eec5 100644
--- a/src/thermoTools/derivedFvPatchFields/totalFlowRateAdvectiveDiffusive/totalFlowRateAdvectiveDiffusiveFvPatchScalarField.C
+++ b/src/thermoTools/derivedFvPatchFields/totalFlowRateAdvectiveDiffusive/totalFlowRateAdvectiveDiffusiveFvPatchScalarField.C
@@ -179,7 +179,7 @@ void Foam::totalFlowRateAdvectiveDiffusiveFvPatchScalarField::updateCoeffs()
 
     if (debug)
     {
-        scalar phi = gSum(-phip*(*this));
+        scalar phi = -gWeightedSum(phip, *this);
 
         Info<< patch().boundaryMesh().mesh().name() << ':'
             << patch().name() << ':'
diff --git a/src/thermoTools/derivedFvPatchFields/turbulentTemperatureRadCoupledMixed/turbulentTemperatureRadCoupledMixedFvPatchScalarField.C b/src/thermoTools/derivedFvPatchFields/turbulentTemperatureRadCoupledMixed/turbulentTemperatureRadCoupledMixedFvPatchScalarField.C
index 81b6f09f68..94789c6abb 100644
--- a/src/thermoTools/derivedFvPatchFields/turbulentTemperatureRadCoupledMixed/turbulentTemperatureRadCoupledMixedFvPatchScalarField.C
+++ b/src/thermoTools/derivedFvPatchFields/turbulentTemperatureRadCoupledMixed/turbulentTemperatureRadCoupledMixedFvPatchScalarField.C
@@ -690,7 +690,7 @@ void turbulentTemperatureRadCoupledMixedFvPatchScalarField::updateCoeffs()
                 htc[i] = qField[i]/deltaT[i];
             }
         }
-        const scalar aveHtc = gSum(htc*patch().magSf())/max(magSf, SMALL);
+        const scalar aveHtc = gWeightedAverage(patch().magSf(), htc);
 
 
         // Calculate heat-transfer coeff based on the second definition
@@ -704,8 +704,7 @@ void turbulentTemperatureRadCoupledMixedFvPatchScalarField::updateCoeffs()
                 htcPatch[i] = qField[i]/deltaTPatch[i];
             }
         }
-        const scalar aveHtcPatch =
-            gSum(htcPatch*patch().magSf())/max(magSf, SMALL);
+        const scalar aveHtcPatch = gWeightedAverage(patch().magSf(), htcPatch);
 
         // Calculate various averages of temperature
         const scalarMinMax TpMinMax = gMinMax(Tp);
diff --git a/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveViewFactor/greyDiffusiveViewFactorFixedValueFvPatchScalarField.C b/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveViewFactor/greyDiffusiveViewFactorFixedValueFvPatchScalarField.C
index c5d3d5ff61..471036865d 100644
--- a/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveViewFactor/greyDiffusiveViewFactorFixedValueFvPatchScalarField.C
+++ b/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveViewFactor/greyDiffusiveViewFactorFixedValueFvPatchScalarField.C
@@ -140,15 +140,17 @@ updateCoeffs()
 
     if (debug)
     {
-        scalar Q = gSum((*this)*patch().magSf());
+        scalar Q = gWeightedSum(patch().magSf(), *this);
+
+        MinMax<scalar> limits = gMinMax(*this);
 
         Info<< patch().boundaryMesh().mesh().name() << ':'
             << patch().name() << ':'
             << this->internalField().name() << " <- "
             << " heat transfer rate:" << Q
             << " wall radiative heat flux "
-            << " min:" << gMin(*this)
-            << " max:" << gMax(*this)
+            << " min:" << limits.min()
+            << " max:" << limits.max()
             << " avg:" << gAverage(*this)
             << endl;
     }
diff --git a/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C b/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C
index 79267e4c80..c1f56cf05e 100644
--- a/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C
+++ b/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C
@@ -1118,7 +1118,7 @@ void Foam::radiation::viewFactor::calculate()
         {
             const scalarField& qrp = qrBf[patchID];
             const scalarField& magSf = mesh_.magSf().boundaryField()[patchID];
-            const scalar heatFlux = gSum(qrp*magSf);
+            const scalar heatFlux = gWeightedSum(magSf, qrp);
 
             Info<< "Total heat transfer rate at patch: "
                 << patchID << " "