diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFPatchTransfer/VoFPatchTransfer.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFPatchTransfer/VoFPatchTransfer.C
index 63cc806366..f3af8611d3 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFPatchTransfer/VoFPatchTransfer.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFPatchTransfer/VoFPatchTransfer.C
@@ -86,6 +86,7 @@ VoFPatchTransfer::VoFPatchTransfer
     wordRes patchNames;
     if (coeffDict_.readIfPresent("patches", patchNames))
     {
+        // Can also use pbm.indices(), but no warnings...
         patchIDs_ = pbm.patchSet(patchNames).sortedToc();
 
         Info<< "        applying to " << patchIDs_.size() << " patches:" << nl;
diff --git a/src/functionObjects/field/Curle/Curle.C b/src/functionObjects/field/Curle/Curle.C
index 9cdcbd3c35..6b53be0723 100644
--- a/src/functionObjects/field/Curle/Curle.C
+++ b/src/functionObjects/field/Curle/Curle.C
@@ -63,7 +63,6 @@ Foam::functionObjects::Curle::Curle
     fvMeshFunctionObject(name, runTime, dict),
     writeFile(mesh_, name),
     pName_("p"),
-    patchSet_(),
     observerPositions_(),
     c0_(0),
     rawFilePtrs_(),
@@ -78,6 +77,8 @@ Foam::functionObjects::Curle::Curle
 
 bool Foam::functionObjects::Curle::read(const dictionary& dict)
 {
+    const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
+
     if (!(fvMeshFunctionObject::read(dict) && writeFile::read(dict)))
     {
         return false;
@@ -85,14 +86,12 @@ bool Foam::functionObjects::Curle::read(const dictionary& dict)
 
     dict.readIfPresent("p", pName_);
 
-    patchSet_ = mesh_.boundaryMesh().patchSet(dict.get<wordRes>("patches"));
+    patchIDs_ = pbm.patchSet(dict.get<wordRes>("patches")).sortedToc();
 
-    if (patchSet_.empty())
+    if (patchIDs_.empty())
     {
         WarningInFunction
-            << "No patches defined"
-            << endl;
-
+            << "No patches defined" << endl;
         return false;
     }
 
@@ -117,8 +116,7 @@ bool Foam::functionObjects::Curle::read(const dictionary& dict)
     if (observerPositions_.empty())
     {
         WarningInFunction
-            << "No observer positions defined"
-            << endl;
+            << "No observer positions defined" << endl;
 
         return false;
     }
@@ -205,7 +203,7 @@ bool Foam::functionObjects::Curle::execute()
 
     scalarField pDash(observerPositions_.size(), 0);
 
-    for (auto patchi : patchSet_)
+    for (const label patchi : patchIDs_)
     {
         const scalarField& pp = pBf[patchi];
         const scalarField& dpdtp = dpdtBf[patchi];
@@ -228,7 +226,7 @@ bool Foam::functionObjects::Curle::execute()
 
     if (surfaceWriterPtr_)
     {
-        if (Pstream::master())
+        if (UPstream::master())
         {
             // Time-aware, with time spliced into the output path
             surfaceWriterPtr_->beginTime(time_);
diff --git a/src/functionObjects/field/Curle/Curle.H b/src/functionObjects/field/Curle/Curle.H
index 19afa8e3de..06beac04f6 100644
--- a/src/functionObjects/field/Curle/Curle.H
+++ b/src/functionObjects/field/Curle/Curle.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2020 OpenCFD Ltd.
+    Copyright (C) 2017-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -124,8 +124,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef functionObjects_Curle_H
-#define functionObjects_Curle_H
+#ifndef Foam_functionObjects_Curle_H
+#define Foam_functionObjects_Curle_H
 
 #include "fvMeshFunctionObject.H"
 #include "writeFile.H"
@@ -163,8 +163,8 @@ class Curle
         //- Name of pressure field; default = p
         word pName_;
 
-        //- Patches to integrate forces over
-        labelHashSet patchSet_;
+        //- List of patches to process
+        labelList patchIDs_;
 
         //- Observer positions
         List<point> observerPositions_;
diff --git a/src/functionObjects/field/binField/binModels/binModel/binModel.C b/src/functionObjects/field/binField/binModels/binModel/binModel.C
index fc31d3e860..e610bada7f 100644
--- a/src/functionObjects/field/binField/binModels/binModel/binModel.C
+++ b/src/functionObjects/field/binField/binModels/binModel/binModel.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021-2022 OpenCFD Ltd.
+    Copyright (C) 2021-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -122,12 +122,8 @@ Foam::binModel::binModel
     mesh_(mesh),
     decomposePatchValues_(false),
     cumulative_(false),
-    coordSysPtr_(),
-    nBin_(1),
-    patchSet_(),
-    fieldNames_(),
-    cellZoneIDs_(),
-    filePtrs_()
+    coordSysPtr_(nullptr),
+    nBin_(1)
 {}
 
 
@@ -135,19 +131,23 @@ Foam::binModel::binModel
 
 bool Foam::binModel::read(const dictionary& dict)
 {
+    const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
+
     if (!functionObjects::writeFile::read(dict))
     {
         return false;
     }
 
-    patchSet_ = mesh_.boundaryMesh().patchSet(dict.get<wordRes>("patches"));
+    // Can also use pbm.indices(), but no warnings...
+    patchIDs_ = pbm.patchSet(dict.get<wordRes>("patches")).sortedToc();
     fieldNames_ = dict.get<wordHashSet>("fields").sortedToc();
 
-    if (dict.found("cellZones"))
+    wordRes zoneNames;
+    if (dict.readIfPresent("cellZones", zoneNames))
     {
         DynamicList<label> zoneIDs;
         DynamicList<wordRe> czUnmatched;
-        for (const auto& cz : dict.get<wordRes>("cellZones"))
+        for (const auto& cz : zoneNames)
         {
             const labelList czi(mesh_.cellZones().indices(cz));
 
diff --git a/src/functionObjects/field/binField/binModels/binModel/binModel.H b/src/functionObjects/field/binField/binModels/binModel/binModel.H
index 3583a0a973..1531507e67 100644
--- a/src/functionObjects/field/binField/binModels/binModel/binModel.H
+++ b/src/functionObjects/field/binField/binModels/binModel/binModel.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021-2022 OpenCFD Ltd.
+    Copyright (C) 2021-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -84,7 +84,7 @@ protected:
         label nBin_;
 
         //- Indices of operand patches
-        labelHashSet patchSet_;
+        labelList patchIDs_;
 
         //- Names of operand fields
         wordList fieldNames_;
diff --git a/src/functionObjects/field/binField/binModels/singleDirectionUniformBin/singleDirectionUniformBin.C b/src/functionObjects/field/binField/binModels/singleDirectionUniformBin/singleDirectionUniformBin.C
index 2140504cc0..e111827617 100644
--- a/src/functionObjects/field/binField/binModels/singleDirectionUniformBin/singleDirectionUniformBin.C
+++ b/src/functionObjects/field/binField/binModels/singleDirectionUniformBin/singleDirectionUniformBin.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021-2022 OpenCFD Ltd.
+    Copyright (C) 2021-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -49,7 +49,7 @@ void Foam::binModels::singleDirectionUniformBin::initialise()
     // Determine extents of patches in a given direction
     scalar geomMin = GREAT;
     scalar geomMax = -GREAT;
-    for (const label patchi : patchSet_)
+    for (const label patchi : patchIDs_)
     {
         const polyPatch& pp = pbm[patchi];
         const scalarField d(pp.faceCentres() & binDir_);
diff --git a/src/functionObjects/field/binField/binModels/singleDirectionUniformBin/singleDirectionUniformBinTemplates.C b/src/functionObjects/field/binField/binModels/singleDirectionUniformBin/singleDirectionUniformBinTemplates.C
index ad5024263c..64ebc66d58 100644
--- a/src/functionObjects/field/binField/binModels/singleDirectionUniformBin/singleDirectionUniformBinTemplates.C
+++ b/src/functionObjects/field/binField/binModels/singleDirectionUniformBin/singleDirectionUniformBinTemplates.C
@@ -154,9 +154,8 @@ bool Foam::binModels::singleDirectionUniformBin::processField
         }
     }
 
-    forAllIters(patchSet_, iter)
+    for (const label patchi : patchIDs_)
     {
-        const label patchi = iter();
         const polyPatch& pp = mesh_.boundaryMesh()[patchi];
         const vectorField np(mesh_.boundary()[patchi].nf());
 
diff --git a/src/functionObjects/field/binField/binModels/uniformBin/uniformBin.C b/src/functionObjects/field/binField/binModels/uniformBin/uniformBin.C
index 77836bf690..064e654eb8 100644
--- a/src/functionObjects/field/binField/binModels/uniformBin/uniformBin.C
+++ b/src/functionObjects/field/binField/binModels/uniformBin/uniformBin.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021-2022 OpenCFD Ltd.
+    Copyright (C) 2021-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -49,7 +49,7 @@ void Foam::binModels::uniformBin::initialise()
     vector geomMin(GREAT, GREAT, GREAT);
     vector geomMax(-GREAT, -GREAT, -GREAT);
 
-    for (const label patchi : patchSet_)
+    for (const label patchi : patchIDs_)
     {
         const polyPatch& pp = pbm[patchi];
         const vectorField ppcs(coordSysPtr_->localPosition(pp.faceCentres()));
@@ -159,19 +159,19 @@ Foam::labelList Foam::binModels::uniformBin::binAddr(const vectorField& d) const
 
 void Foam::binModels::uniformBin::setBinsAddressing()
 {
-    faceToBin_.setSize(mesh_.nBoundaryFaces());
+    faceToBin_.resize_nocopy(mesh_.nBoundaryFaces());
     faceToBin_ = -1;
 
-    forAllIters(patchSet_, iter)
+    for (const label patchi : patchIDs_)
     {
-        const polyPatch& pp = mesh_.boundaryMesh()[iter()];
+        const polyPatch& pp = mesh_.boundaryMesh()[patchi];
         const label i0 = pp.start() - mesh_.nInternalFaces();
 
         SubList<label>(faceToBin_, pp.size(), i0) =
             binAddr(coordSysPtr_->localPosition(pp.faceCentres()));
     }
 
-    cellToBin_.setSize(mesh_.nCells());
+    cellToBin_.resize_nocopy(mesh_.nCells());
     cellToBin_ = -1;
 
     for (const label zonei : cellZoneIDs_)
diff --git a/src/functionObjects/field/binField/binModels/uniformBin/uniformBinTemplates.C b/src/functionObjects/field/binField/binModels/uniformBin/uniformBinTemplates.C
index c9edd55a13..b4dd98f225 100644
--- a/src/functionObjects/field/binField/binModels/uniformBin/uniformBinTemplates.C
+++ b/src/functionObjects/field/binField/binModels/uniformBin/uniformBinTemplates.C
@@ -142,9 +142,8 @@ bool Foam::binModels::uniformBin::processField(const label fieldi)
         }
     }
 
-    forAllIters(patchSet_, iter)
+    for (const label patchi : patchIDs_)
     {
-        const label patchi = iter();
         const polyPatch& pp = mesh_.boundaryMesh()[patchi];
         const vectorField np(mesh_.boundary()[patchi].nf());
 
diff --git a/src/functionObjects/field/fieldExtents/fieldExtents.C b/src/functionObjects/field/fieldExtents/fieldExtents.C
index db42a4ac55..d2adc5f4c2 100644
--- a/src/functionObjects/field/fieldExtents/fieldExtents.C
+++ b/src/functionObjects/field/fieldExtents/fieldExtents.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -111,8 +111,7 @@ Foam::functionObjects::fieldExtents::fieldExtents
     internalField_(true),
     threshold_(0),
     C0_(Zero),
-    fieldSet_(mesh_),
-    patchIDs_()
+    fieldSet_(mesh_)
 {
     read(dict);
 
@@ -125,6 +124,8 @@ Foam::functionObjects::fieldExtents::fieldExtents
 
 bool Foam::functionObjects::fieldExtents::read(const dictionary& dict)
 {
+    const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
+
     if (fvMeshFunctionObject::read(dict) && writeFile::read(dict))
     {
         dict.readIfPresent<bool>("internalField", internalField_);
@@ -133,28 +134,31 @@ bool Foam::functionObjects::fieldExtents::read(const dictionary& dict)
 
         dict.readIfPresent<vector>("referencePosition", C0_);
 
-        patchIDs_.clear();
-        const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
-
         wordRes patchNames;
         if (dict.readIfPresent("patches", patchNames))
         {
-            for (const wordRe& name : patchNames)
-            {
-                patchIDs_.insert(pbm.indices(name));
-            }
+            patchIDs_ = pbm.indices(patchNames);
         }
         else
         {
-            // Add all non-processor and non-empty patches
+            labelHashSet patchSet(2*pbm.size());
+
+            // All non-processor and non-empty patches
             forAll(pbm, patchi)
             {
                 const polyPatch& pp = pbm[patchi];
-                if (!isA<processorPolyPatch>(pp) && !isA<emptyPolyPatch>(pp))
+
+                if
+                (
+                    !isA<processorPolyPatch>(pp)
+                 && !isA<emptyPolyPatch>(pp)
+                )
                 {
-                    patchIDs_.insert(patchi);
+                    patchSet.insert(patchi);
                 }
             }
+
+            patchIDs_ = patchSet.sortedToc();
         }
 
         if (!internalField_ && patchIDs_.empty())
diff --git a/src/functionObjects/field/fieldExtents/fieldExtents.H b/src/functionObjects/field/fieldExtents/fieldExtents.H
index d042870fa3..8b8f0cd39a 100644
--- a/src/functionObjects/field/fieldExtents/fieldExtents.H
+++ b/src/functionObjects/field/fieldExtents/fieldExtents.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -99,8 +99,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef functionObjects_fieldExtents_H
-#define functionObjects_fieldExtents_H
+#ifndef Foam_functionObjects_fieldExtents_H
+#define Foam_functionObjects_fieldExtents_H
 
 #include "fvMeshFunctionObject.H"
 #include "writeFile.H"
@@ -122,7 +122,6 @@ class fieldExtents
     public fvMeshFunctionObject,
     public writeFile
 {
-
 protected:
 
     // Protected Data
@@ -140,7 +139,7 @@ protected:
         volFieldSelection fieldSet_;
 
         //- Patches to assess
-        labelHashSet patchIDs_;
+        labelList patchIDs_;
 
 
     // Protected Member Functions
diff --git a/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/ReynoldsAnalogy/ReynoldsAnalogy.C b/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/ReynoldsAnalogy/ReynoldsAnalogy.C
index a0e5a63b67..5f363692a7 100644
--- a/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/ReynoldsAnalogy/ReynoldsAnalogy.C
+++ b/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/ReynoldsAnalogy/ReynoldsAnalogy.C
@@ -174,7 +174,7 @@ Foam::heatTransferCoeffModels::ReynoldsAnalogy::Cf() const
     const volSymmTensorField R(devReff());
     const volSymmTensorField::Boundary& Rbf = R.boundaryField();
 
-    for (const label patchi : patchSet_)
+    for (const label patchi : patchIDs_)
     {
         const fvPatchVectorField& Up = Ubf[patchi];
 
@@ -202,7 +202,7 @@ void Foam::heatTransferCoeffModels::ReynoldsAnalogy::htc
 
     volScalarField::Boundary& htcBf = htc.boundaryFieldRef();
 
-    for (const label patchi : patchSet_)
+    for (const label patchi : patchIDs_)
     {
         tmp<scalarField> trhop = rho(patchi);
         tmp<scalarField> tCpp = Cp(patchi);
diff --git a/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/faceZoneReferenceTemperature/faceZoneReferenceTemperature.C b/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/faceZoneReferenceTemperature/faceZoneReferenceTemperature.C
index 66f6a390a6..7e75ccb304 100644
--- a/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/faceZoneReferenceTemperature/faceZoneReferenceTemperature.C
+++ b/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/faceZoneReferenceTemperature/faceZoneReferenceTemperature.C
@@ -190,7 +190,7 @@ void Foam::heatTransferCoeffModels::faceZoneReferenceTemperature::htc
     const scalar Tref = faceZoneAverageTemperature();
 
     // Calculate heat-transfer coefficient boundary fields for current region
-    for (const label patchi : patchSet_)
+    for (const label patchi : patchIDs_)
     {
         htcBf[patchi] = q[patchi]/(Tref - Tbf[patchi] + ROOTVSMALL);
     }
diff --git a/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/fixedReferenceTemperature/fixedReferenceTemperature.C b/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/fixedReferenceTemperature/fixedReferenceTemperature.C
index f37ce5b5ff..65c10017d5 100644
--- a/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/fixedReferenceTemperature/fixedReferenceTemperature.C
+++ b/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/fixedReferenceTemperature/fixedReferenceTemperature.C
@@ -58,7 +58,7 @@ void Foam::heatTransferCoeffModels::fixedReferenceTemperature::htc
     const scalar eps = ROOTVSMALL;
 
     volScalarField::Boundary& htcBf = htc.boundaryFieldRef();
-    for (const label patchi : patchSet_)
+    for (const label patchi : patchIDs_)
     {
         htcBf[patchi] = q[patchi]/(TRef_ - Tbf[patchi] + eps);
     }
diff --git a/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/heatTransferCoeffModel/heatTransferCoeffModel.C b/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/heatTransferCoeffModel/heatTransferCoeffModel.C
index b5d342465c..125fccea7c 100644
--- a/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/heatTransferCoeffModel/heatTransferCoeffModel.C
+++ b/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/heatTransferCoeffModel/heatTransferCoeffModel.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2020 OpenCFD Ltd.
+    Copyright (C) 2017-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -50,7 +50,6 @@ Foam::heatTransferCoeffModel::heatTransferCoeffModel
 )
 :
     mesh_(mesh),
-    patchSet_(),
     TName_(TName),
     qrName_("qr")
 {}
@@ -86,7 +85,7 @@ Foam::heatTransferCoeffModel::q() const
         const volScalarField alphaEff(turb.alphaEff());
         const volScalarField::Boundary& alphaEffbf = alphaEff.boundaryField();
 
-        for (const label patchi : patchSet_)
+        for (const label patchi : patchIDs_)
         {
             q[patchi] = alphaEffbf[patchi]*hebf[patchi].snGrad();
         }
@@ -103,7 +102,7 @@ Foam::heatTransferCoeffModel::q() const
         const volScalarField& alpha(thermo.alpha());
         const volScalarField::Boundary& alphabf = alpha.boundaryField();
 
-        for (const label patchi : patchSet_)
+        for (const label patchi : patchIDs_)
         {
             q[patchi] = alphabf[patchi]*hebf[patchi].snGrad();
         }
@@ -124,7 +123,7 @@ Foam::heatTransferCoeffModel::q() const
     {
         const volScalarField::Boundary& qrbf = qrPtr->boundaryField();
 
-        for (const label patchi : patchSet_)
+        for (const label patchi : patchIDs_)
         {
             q[patchi] += qrbf[patchi];
         }
@@ -148,10 +147,12 @@ bool Foam::heatTransferCoeffModel::calc
 
 bool Foam::heatTransferCoeffModel::read(const dictionary& dict)
 {
-    patchSet_ = mesh_.boundaryMesh().patchSet(dict.get<wordRes>("patches"));
+    const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
 
     dict.readIfPresent("qr", qrName_);
 
+    patchIDs_ = pbm.patchSet(dict.get<wordRes>("patches")).sortedToc();
+
     return true;
 }
 
diff --git a/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/heatTransferCoeffModel/heatTransferCoeffModel.H b/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/heatTransferCoeffModel/heatTransferCoeffModel.H
index a3abbb91ce..9d44e94dbb 100644
--- a/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/heatTransferCoeffModel/heatTransferCoeffModel.H
+++ b/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/heatTransferCoeffModel/heatTransferCoeffModel.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2022 OpenCFD Ltd.
+    Copyright (C) 2017-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -67,8 +67,8 @@ protected:
         //- Const reference to the mesh
         const fvMesh& mesh_;
 
-        //- List of (wall) patches to process
-        labelHashSet patchSet_;
+        //- List of (wall) patches to process (selected by name)
+        labelList patchIDs_;
 
         //- Name of temperature field
         const word TName_;
@@ -152,9 +152,9 @@ public:
         }
 
         //- Return const reference to wall patches to process
-        const labelHashSet& patchSet() const noexcept
+        const labelList& patchIDs() const noexcept
         {
-            return patchSet_;
+            return patchIDs_;
         }
 
         //- Return const reference to name of temperature field
diff --git a/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/localReferenceTemperature/localReferenceTemperature.C b/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/localReferenceTemperature/localReferenceTemperature.C
index 5810a00083..cd9dfbb6f3 100644
--- a/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/localReferenceTemperature/localReferenceTemperature.C
+++ b/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/localReferenceTemperature/localReferenceTemperature.C
@@ -59,7 +59,7 @@ void Foam::heatTransferCoeffModels::localReferenceTemperature::htc
 
     volScalarField::Boundary& htcBf = htc.boundaryFieldRef();
 
-    for (const label patchi : patchSet_)
+    for (const label patchi : patchIDs_)
     {
         tmp<scalarField> tTc = Tbf[patchi].patchInternalField();
         htcBf[patchi] = q[patchi]/(tTc - Tbf[patchi] + eps);
diff --git a/src/functionObjects/field/heatTransferCoeff/multiphaseInterHtcModel/multiphaseInterHtcModel.C b/src/functionObjects/field/heatTransferCoeff/multiphaseInterHtcModel/multiphaseInterHtcModel.C
index 0dc05eede6..206827e72a 100644
--- a/src/functionObjects/field/heatTransferCoeff/multiphaseInterHtcModel/multiphaseInterHtcModel.C
+++ b/src/functionObjects/field/heatTransferCoeff/multiphaseInterHtcModel/multiphaseInterHtcModel.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2022 OpenCFD Ltd.
+    Copyright (C) 2022-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -89,7 +89,7 @@ Foam::functionObjects::multiphaseInterHtcModel::q() const
 
     const multiphaseInterSystem& fluid = *fluidPtr;
 
-    for (const label patchi : htcModelPtr_->patchSet())
+    for (const label patchi : htcModelPtr_->patchIDs())
     {
         q[patchi] += fluid.kappaEff(patchi)()*Tbf[patchi].snGrad();
     }
@@ -103,7 +103,7 @@ Foam::functionObjects::multiphaseInterHtcModel::q() const
     {
         const volScalarField::Boundary& qrbf = qrPtr->boundaryField();
 
-        for (const label patchi : htcModelPtr_->patchSet())
+        for (const label patchi : htcModelPtr_->patchIDs())
         {
             q[patchi] += qrbf[patchi];
         }
diff --git a/src/functionObjects/field/heatTransferCoeff/reactingEulerHtcModel/reactingEulerHtcModel.C b/src/functionObjects/field/heatTransferCoeff/reactingEulerHtcModel/reactingEulerHtcModel.C
index a4ebe4b2df..fa8108e9b3 100644
--- a/src/functionObjects/field/heatTransferCoeff/reactingEulerHtcModel/reactingEulerHtcModel.C
+++ b/src/functionObjects/field/heatTransferCoeff/reactingEulerHtcModel/reactingEulerHtcModel.C
@@ -90,7 +90,7 @@ Foam::functionObjects::reactingEulerHtcModel::q() const
 
     const phaseSystem& fluid = *fluidPtr;
 
-    for (const label patchi : htcModelPtr_->patchSet())
+    for (const label patchi : htcModelPtr_->patchIDs())
     {
         for (const phaseModel& phase : fluid.phases())
         {
@@ -112,7 +112,7 @@ Foam::functionObjects::reactingEulerHtcModel::q() const
     {
         const volScalarField::Boundary& qrbf = qrPtr->boundaryField();
 
-        for (const label patchi : htcModelPtr_->patchSet())
+        for (const label patchi : htcModelPtr_->patchIDs())
         {
             q[patchi] += qrbf[patchi];
         }
diff --git a/src/functionObjects/field/nearWallFields/nearWallFields.C b/src/functionObjects/field/nearWallFields/nearWallFields.C
index a6e32cdb8e..5aeaf99c93 100644
--- a/src/functionObjects/field/nearWallFields/nearWallFields.C
+++ b/src/functionObjects/field/nearWallFields/nearWallFields.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2015-2022 OpenCFD Ltd.
+    Copyright (C) 2015-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -51,13 +51,13 @@ void Foam::functionObjects::nearWallFields::calcAddressing()
 {
     // Count number of faces
     label nPatchFaces = 0;
-    for (const label patchi : patchSet_)
+    for (const label patchi : patchIDs_)
     {
         nPatchFaces += mesh_.boundary()[patchi].size();
     }
 
     // Global indexing
-    globalIndex globalWalls(nPatchFaces);
+    const globalIndex globalWalls(nPatchFaces);
 
     DebugInFunction << "nPatchFaces: " << globalWalls.totalSize() << endl;
 
@@ -72,7 +72,7 @@ void Foam::functionObjects::nearWallFields::calcAddressing()
     // Add particles to track to sample locations
     nPatchFaces = 0;
 
-    for (const label patchi : patchSet_)
+    for (const label patchi : patchIDs_)
     {
         const fvPatch& patch = mesh_.boundary()[patchi];
 
@@ -260,16 +260,15 @@ Foam::functionObjects::nearWallFields::nearWallFields
 
 bool Foam::functionObjects::nearWallFields::read(const dictionary& dict)
 {
+    const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
+
     fvMeshFunctionObject::read(dict);
 
     dict.readEntry("fields", fieldSet_);
     dict.readEntry("distance", distance_);
 
-    patchSet_ =
-        mesh_.boundaryMesh().patchSet
-        (
-            dict.get<wordRes>("patches")
-        );
+    // Can also use pbm.indices(), but no warnings...
+    patchIDs_ = pbm.patchSet(dict.get<wordRes>("patches")).sortedToc();
 
 
     // Clear out any previously loaded fields
diff --git a/src/functionObjects/field/nearWallFields/nearWallFields.H b/src/functionObjects/field/nearWallFields/nearWallFields.H
index 2988672477..587f601d1b 100644
--- a/src/functionObjects/field/nearWallFields/nearWallFields.H
+++ b/src/functionObjects/field/nearWallFields/nearWallFields.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2020 OpenCFD Ltd.
+    Copyright (C) 2015-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -90,8 +90,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef functionObjects_nearWallFields_H
-#define functionObjects_nearWallFields_H
+#ifndef Foam_functionObjects_nearWallFields_H
+#define Foam_functionObjects_nearWallFields_H
 
 #include "fvMeshFunctionObject.H"
 #include "volFields.H"
@@ -123,7 +123,7 @@ protected:
             List<Tuple2<word, word>> fieldSet_;
 
             //- Patches to sample
-            labelHashSet patchSet_;
+            labelList patchIDs_;
 
             //- Distance away from wall
             scalar distance_;
diff --git a/src/functionObjects/field/nearWallFields/nearWallFieldsTemplates.C b/src/functionObjects/field/nearWallFields/nearWallFieldsTemplates.C
index aee5fecc8d..e8c20782f1 100644
--- a/src/functionObjects/field/nearWallFields/nearWallFieldsTemplates.C
+++ b/src/functionObjects/field/nearWallFields/nearWallFieldsTemplates.C
@@ -74,7 +74,7 @@ void Foam::functionObjects::nearWallFields::createFields
                     (
                         io,
                         fld,
-                        patchSet_.toc(),
+                        patchIDs_,
                         fvPatchFieldBase::calculatedType()
                     )
                 );
@@ -120,7 +120,7 @@ void Foam::functionObjects::nearWallFields::sampleBoundaryField
 
     // Pick up data
     label nPatchFaces = 0;
-    for (const label patchi : patchSet_)
+    for (const label patchi : patchIDs_)
     {
         fvPatchField<Type>& pfld = fldBf[patchi];
 
diff --git a/src/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C b/src/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C
index ebdf67cf78..bb60583fc5 100644
--- a/src/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C
+++ b/src/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C
@@ -61,7 +61,7 @@ template<class Type>
 static Map<Type> regionSum(const regionSplit& regions, const Field<Type>& fld)
 {
     // Per region the sum of fld
-    Map<Type> regionToSum(regions.nRegions()/Pstream::nProcs());
+    Map<Type> regionToSum(regions.nRegions()/UPstream::nProcs());
 
     forAll(fld, celli)
     {
@@ -75,6 +75,23 @@ static Map<Type> regionSum(const regionSplit& regions, const Field<Type>& fld)
 }
 
 
+static Map<label> regionSum(const regionSplit& regions, const label nCells)
+{
+    // Per region the sum of fld
+    Map<label> regionToSum(regions.nRegions()/UPstream::nProcs());
+
+    for (label celli = 0; celli < nCells; ++celli)
+    {
+        const label regioni = regions[celli];
+        ++regionToSum(regioni);
+    }
+
+    Pstream::mapCombineReduce(regionToSum, plusEqOp<label>());
+
+    return regionToSum;
+}
+
+
 template<class Type>
 static List<Type> extractData(const labelUList& keys, const Map<Type>& regionData)
 {
@@ -95,7 +112,7 @@ static List<Type> extractData(const labelUList& keys, const Map<Type>& regionDat
 void Foam::functionObjects::regionSizeDistribution::writeAlphaFields
 (
     const regionSplit& regions,
-    const Map<label>& patchRegions,
+    const labelHashSet& keepRegions,
     const Map<scalar>& regionVolume,
     const volScalarField& alpha
 ) const
@@ -138,11 +155,11 @@ void Foam::functionObjects::regionSizeDistribution::writeAlphaFields
     );
 
 
-    // Knock out any cell not in patchRegions
+    // Knock out any cell not in keepRegions (patch regions)
     forAll(liquidCore, celli)
     {
         const label regioni = regions[celli];
-        if (patchRegions.found(regioni))
+        if (keepRegions.found(regioni))
         {
             backgroundAlpha[celli] = 0;
         }
@@ -178,7 +195,7 @@ void Foam::functionObjects::regionSizeDistribution::writeAlphaFields
 }
 
 
-Foam::Map<Foam::label>
+Foam::labelHashSet
 Foam::functionObjects::regionSizeDistribution::findPatchRegions
 (
     const regionSplit& regions
@@ -187,37 +204,23 @@ Foam::functionObjects::regionSizeDistribution::findPatchRegions
     // Mark all regions starting at patches
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    // Count number of patch faces (just for initial sizing)
-    const labelHashSet patchIDs(mesh_.boundaryMesh().patchSet(patchNames_));
+    labelHashSet patchRegions(2*regions.nRegions());
 
-    label nPatchFaces = 0;
-    for (const label patchi : patchIDs)
-    {
-        nPatchFaces += mesh_.boundaryMesh()[patchi].size();
-    }
+    labelHashSet patchSet(mesh_.boundaryMesh().patchSet(patchNames_));
 
-
-    Map<label> patchRegions(nPatchFaces);
-    for (const label patchi : patchIDs)
+    for (const label patchi : patchSet)
     {
         const polyPatch& pp = mesh_.boundaryMesh()[patchi];
 
-        // Collect all regions on the patch
-        const labelList& faceCells = pp.faceCells();
-
-        for (const label celli : faceCells)
+        // All regions connected to the patch
+        for (const label celli : pp.faceCells())
         {
-            patchRegions.insert
-            (
-                regions[celli],
-                Pstream::myProcNo()     // dummy value
-            );
+            patchRegions.insert(regions[celli]);
         }
     }
 
-
-    // Make sure all the processors have the same set of regions
-    Pstream::mapCombineReduce(patchRegions, minEqOp<label>());
+    // Ensure all processors have the same set of regions
+    Pstream::combineReduce(patchRegions, plusEqOp<labelHashSet>());
 
     return patchRegions;
 }
@@ -230,19 +233,15 @@ Foam::functionObjects::regionSizeDistribution::divide
     const scalarField& denom
 )
 {
-    auto tresult = tmp<scalarField>::New(num.size());
+    auto tresult = tmp<scalarField>::New(num.size(), Zero);
     auto& result = tresult.ref();
 
     forAll(denom, i)
     {
-        if (denom[i] != 0)
+        if (ROOTVSMALL < Foam::mag(denom[i]))
         {
             result[i] = num[i]/denom[i];
         }
-        else
-        {
-            result[i] = 0;
-        }
     }
     return tresult;
 }
@@ -258,7 +257,7 @@ void Foam::functionObjects::regionSizeDistribution::writeGraphs
     const coordSet& coords              // graph data for bins
 ) const
 {
-    if (Pstream::master())
+    if (UPstream::master())
     {
         // Calculate per-bin average
         scalarField binSum(nBins_, Zero);
@@ -449,8 +448,8 @@ bool Foam::functionObjects::regionSizeDistribution::write()
                     alphaName_,
                     mesh_.time().timeName(),
                     mesh_,
-                    IOobject::MUST_READ,
-                    IOobject::NO_WRITE
+                    IOobjectOption::MUST_READ,
+                    IOobjectOption::NO_WRITE
                 ),
                 mesh_
             )
@@ -541,14 +540,15 @@ bool Foam::functionObjects::regionSizeDistribution::write()
                 "region",
                 mesh_.time().timeName(),
                 mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
+                IOobjectOption::NO_READ,
+                IOobjectOption::NO_WRITE,
+                IOobjectOption::NO_REGISTER
             ),
             mesh_,
             dimensionedScalar(dimless, Zero)
         );
-        Info<< "    Dumping region as volScalarField to " << region.name()
-            << endl;
+        Info<< "    Dumping region as volScalarField to "
+            << region.name() << endl;
 
         forAll(regions, celli)
         {
@@ -560,21 +560,13 @@ bool Foam::functionObjects::regionSizeDistribution::write()
 
 
     // Determine regions connected to supplied patches
-    Map<label> patchRegions(findPatchRegions(regions));
-
+    const labelHashSet patchRegions(findPatchRegions(regions));
 
     // Sum all regions
     const scalarField alphaVol(alpha.primitiveField()*mesh_.V());
     Map<scalar> allRegionVolume(regionSum(regions, mesh_.V()));
     Map<scalar> allRegionAlphaVolume(regionSum(regions, alphaVol));
-    Map<label> allRegionNumCells
-    (
-        regionSum
-        (
-            regions,
-            labelField(mesh_.nCells(), 1.0)
-        )
-    );
+    Map<label> allRegionNumCells(regionSum(regions, mesh_.nCells()));
 
     if (debug)
     {
@@ -623,9 +615,8 @@ bool Foam::functionObjects::regionSizeDistribution::write()
             << token::TAB << "Volume(" << alpha.name() << "):"
             << nl;
 
-        forAllConstIters(patchRegions, iter)
+        for (const label regioni : patchRegions.sortedToc())
         {
-            const label regioni = iter.key();
             Info<< "    " << token::TAB << regioni
                 << token::TAB << allRegionVolume[regioni]
                 << token::TAB << allRegionAlphaVolume[regioni] << nl;
@@ -786,7 +777,7 @@ bool Foam::functionObjects::regionSizeDistribution::write()
             }
 
             // Write
-            if (Pstream::master())
+            if (UPstream::master())
             {
                 // Construct mids of bins for plotting
                 pointField xBin(nDownstreamBins_, Zero);
@@ -863,7 +854,7 @@ bool Foam::functionObjects::regionSizeDistribution::write()
         }
 
         // Write counts
-        if (Pstream::master())
+        if (UPstream::master())
         {
             auto& writer = formatterPtr_();
             writer.nFields(1);
diff --git a/src/functionObjects/field/regionSizeDistribution/regionSizeDistribution.H b/src/functionObjects/field/regionSizeDistribution/regionSizeDistribution.H
index 0bdc20cbcf..b9cd7f96de 100644
--- a/src/functionObjects/field/regionSizeDistribution/regionSizeDistribution.H
+++ b/src/functionObjects/field/regionSizeDistribution/regionSizeDistribution.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2022 OpenCFD Ltd.
+    Copyright (C) 2016-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -257,13 +257,13 @@ class regionSizeDistribution
         void writeAlphaFields
         (
             const regionSplit& regions,
-            const Map<label>& keepRegions,
+            const labelHashSet& keepRegions,
             const Map<scalar>& regionVolume,
             const volScalarField& alpha
         ) const;
 
         //- Mark all regions starting at patches
-        Map<label> findPatchRegions(const regionSplit&) const;
+        labelHashSet findPatchRegions(const regionSplit&) const;
 
         //- Helper: divide if denom != 0
         static tmp<scalarField> divide(const scalarField&, const scalarField&);
diff --git a/src/functionObjects/field/wallHeatFlux/wallHeatFlux.C b/src/functionObjects/field/wallHeatFlux/wallHeatFlux.C
index ca638e238f..a86d780a0f 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-2022 OpenCFD Ltd.
+    Copyright (C) 2016-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -75,7 +75,7 @@ void Foam::functionObjects::wallHeatFlux::calcHeatFlux
 
     const volScalarField::Boundary& alphaBf = alpha.boundaryField();
 
-    for (const label patchi : patchSet_)
+    for (const label patchi : patchIDs_)
     {
         wallHeatFluxBf[patchi] = alphaBf[patchi]*heBf[patchi].snGrad();
     }
@@ -87,7 +87,7 @@ void Foam::functionObjects::wallHeatFlux::calcHeatFlux
     {
         const volScalarField::Boundary& radHeatFluxBf = qrPtr->boundaryField();
 
-        for (const label patchi : patchSet_)
+        for (const label patchi : patchIDs_)
         {
             wallHeatFluxBf[patchi] -= radHeatFluxBf[patchi];
         }
@@ -106,7 +106,6 @@ Foam::functionObjects::wallHeatFlux::wallHeatFlux
 :
     fvMeshFunctionObject(name, runTime, dict),
     writeFile(obr_, name, typeName, dict),
-    patchSet_(),
     qrName_("qr")
 {
     volScalarField* wallHeatFluxPtr
@@ -118,9 +117,9 @@ Foam::functionObjects::wallHeatFlux::wallHeatFlux
                 scopedName(typeName),
                 mesh_.time().timeName(),
                 mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                IOobject::REGISTER
+                IOobjectOption::NO_READ,
+                IOobjectOption::NO_WRITE,
+                IOobjectOption::REGISTER
             ),
             mesh_,
             dimensionedScalar(dimMass/pow3(dimTime), Zero)
@@ -139,55 +138,58 @@ Foam::functionObjects::wallHeatFlux::wallHeatFlux
 
 bool Foam::functionObjects::wallHeatFlux::read(const dictionary& dict)
 {
+    const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
+
     fvMeshFunctionObject::read(dict);
     writeFile::read(dict);
 
-    const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
-
-    patchSet_ =
-        mesh_.boundaryMesh().patchSet
-        (
-            dict.getOrDefault<wordRes>("patches", wordRes())
-        );
-
     dict.readIfPresent("qr", qrName_);
 
-    Info<< type() << " " << name() << ":" << nl;
-
-    if (patchSet_.empty())
+    wordRes patchNames;
+    labelHashSet patchSet(0);
+    if (dict.readIfPresent("patches", patchNames) && !patchNames.empty())
     {
-        forAll(pbm, patchi)
-        {
-            if (isA<wallPolyPatch>(pbm[patchi]))
-            {
-                patchSet_.insert(patchi);
-            }
-        }
+        patchSet = pbm.patchSet(patchNames);
+    }
 
-        Info<< "    processing all wall patches" << nl << endl;
+    labelHashSet allWalls(pbm.findPatchIDs<wallPolyPatch>());
+
+    Info<< type() << ' ' << name() << ':' << nl;
+
+    if (patchSet.empty())
+    {
+        patchIDs_ = allWalls.sortedToc();
+
+        Info<< "    processing all (" << patchIDs_.size()
+            << ") wall patches" << nl << endl;
     }
     else
     {
-        Info<< "    processing wall patches: " << nl;
-        labelHashSet filteredPatchSet;
-        for (const label patchi : patchSet_)
+        allWalls &= patchSet;
+        patchSet -= allWalls;
+        patchIDs_ = allWalls.sortedToc();
+
+        if (!patchSet.empty())
         {
-            if (isA<wallPolyPatch>(pbm[patchi]))
+            WarningInFunction
+                << "Requested wall heat-flux on ("
+                << patchSet.size() << ") non-wall patches:" << nl;
+
+            for (const label patchi : patchSet.sortedToc())
             {
-                filteredPatchSet.insert(patchi);
-                Info<< "        " << pbm[patchi].name() << endl;
-            }
-            else
-            {
-                WarningInFunction
-                    << "Requested wall heat-flux on non-wall boundary "
-                    << "type patch: " << pbm[patchi].name() << endl;
+                Info<< "        " << pbm[patchi].name() << nl;
             }
+            Info<< nl;
         }
 
-        Info<< endl;
+        Info<< "    processing (" << patchIDs_.size()
+            << ") wall patches:" << nl;
 
-        patchSet_ = filteredPatchSet;
+        for (const label patchi : patchIDs_)
+        {
+            Info<< "        " << pbm[patchi].name() << nl;
+        }
+        Info<< endl;
     }
 
     return true;
@@ -258,18 +260,18 @@ bool Foam::functionObjects::wallHeatFlux::execute()
             << "database" << exit(FatalError);
     }
 
+
     const fvPatchList& patches = mesh_.boundary();
 
     const surfaceScalarField::Boundary& magSf = mesh_.magSf().boundaryField();
 
-    for (const label patchi : patchSet_)
+    for (const label patchi : patchIDs_)
     {
         const fvPatch& pp = patches[patchi];
 
         const scalarField& hfp = wallHeatFlux.boundaryField()[patchi];
 
-        const scalar minHfp = gMin(hfp);
-        const scalar maxHfp = gMax(hfp);
+        const MinMax<scalar> limits = gMinMax(hfp);
         const scalar integralHfp = gSum(magSf[patchi]*hfp);
 
         if (Pstream::master())
@@ -278,21 +280,21 @@ bool Foam::functionObjects::wallHeatFlux::execute()
 
             file()
                 << token::TAB << pp.name()
-                << token::TAB << minHfp
-                << token::TAB << maxHfp
+                << token::TAB << limits.min()
+                << token::TAB << limits.max()
                 << token::TAB << integralHfp
                 << endl;
         }
 
         Log << "    min/max/integ(" << pp.name() << ") = "
-            << minHfp << ", " << maxHfp << ", " << integralHfp << endl;
+            << limits.min() << ", " << limits.max()
+            << ", " << integralHfp << endl;
 
-        this->setResult("min(" + pp.name() + ")", minHfp);
-        this->setResult("max(" + pp.name() + ")", maxHfp);
+        this->setResult("min(" + pp.name() + ")", limits.min());
+        this->setResult("max(" + pp.name() + ")", limits.max());
         this->setResult("int(" + pp.name() + ")", integralHfp);
     }
 
-
     return true;
 }
 
@@ -302,7 +304,7 @@ bool Foam::functionObjects::wallHeatFlux::write()
     const auto& wallHeatFlux =
         lookupObject<volScalarField>(scopedName(typeName));
 
-    Log << type() << " " << name() << " write:" << nl
+    Log << type() << ' ' << name() << " write:" << nl
         << "    writing field " << wallHeatFlux.name() << endl;
 
     wallHeatFlux.write();
diff --git a/src/functionObjects/field/wallHeatFlux/wallHeatFlux.H b/src/functionObjects/field/wallHeatFlux/wallHeatFlux.H
index b4f8e12a37..7938d13d94 100644
--- a/src/functionObjects/field/wallHeatFlux/wallHeatFlux.H
+++ b/src/functionObjects/field/wallHeatFlux/wallHeatFlux.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 OpenFOAM Foundation
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2016-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -89,8 +89,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef functionObjects_wallHeatFlux_H
-#define functionObjects_wallHeatFlux_H
+#ifndef Foam_functionObjects_wallHeatFlux_H
+#define Foam_functionObjects_wallHeatFlux_H
 
 #include "fvMeshFunctionObject.H"
 #include "writeFile.H"
@@ -117,8 +117,8 @@ protected:
 
     // Protected Data
 
-        //- Optional list of wall patches to process
-        labelHashSet patchSet_;
+        //- Wall patches to process (optionally filtered by name)
+        labelList patchIDs_;
 
         //- Name of radiative heat flux name
         word qrName_;
diff --git a/src/functionObjects/field/wallShearStress/wallShearStress.C b/src/functionObjects/field/wallShearStress/wallShearStress.C
index 0f27187334..83ce52c202 100644
--- a/src/functionObjects/field/wallShearStress/wallShearStress.C
+++ b/src/functionObjects/field/wallShearStress/wallShearStress.C
@@ -67,7 +67,7 @@ void Foam::functionObjects::wallShearStress::calcShearStress
 {
     shearStress.dimensions().reset(Reff.dimensions());
 
-    for (const label patchi : patchSet_)
+    for (const label patchi : patchIDs_)
     {
         vectorField& ssp = shearStress.boundaryFieldRef()[patchi];
         const vectorField& Sfp = mesh_.Sf().boundaryField()[patchi];
@@ -105,9 +105,9 @@ Foam::functionObjects::wallShearStress::wallShearStress
                 scopedName(typeName),
                 mesh_.time().timeName(),
                 mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                IOobject::REGISTER
+                IOobjectOption::NO_READ,
+                IOobjectOption::NO_WRITE,
+                IOobjectOption::REGISTER
             ),
             mesh_,
             dimensionedVector(sqr(dimLength)/sqr(dimTime), Zero)
@@ -130,48 +130,51 @@ bool Foam::functionObjects::wallShearStress::read(const dictionary& dict)
 
     const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
 
-    patchSet_ =
-        mesh_.boundaryMesh().patchSet
-        (
-            dict.getOrDefault<wordRes>("patches", wordRes())
-        );
-
-    Info<< type() << " " << name() << ":" << nl;
-
-    if (patchSet_.empty())
+    wordRes patchNames;
+    labelHashSet patchSet(0);
+    if (dict.readIfPresent("patches", patchNames) && !patchNames.empty())
     {
-        forAll(pbm, patchi)
-        {
-            if (isA<wallPolyPatch>(pbm[patchi]))
-            {
-                patchSet_.insert(patchi);
-            }
-        }
+        patchSet = pbm.patchSet(patchNames);
+    }
 
-        Info<< "    processing all wall patches" << nl << endl;
+    labelHashSet allWalls(pbm.findPatchIDs<wallPolyPatch>());
+
+    Info<< type() << ' ' << name() << ':' << nl;
+
+    if (patchSet.empty())
+    {
+        patchIDs_ = allWalls.sortedToc();
+
+        Info<< "    processing all (" << patchIDs_.size()
+            << ") wall patches" << nl << endl;
     }
     else
     {
-        Info<< "    processing wall patches: " << nl;
-        labelHashSet filteredPatchSet;
-        for (const label patchi : patchSet_)
+        allWalls &= patchSet;
+        patchSet -= allWalls;
+        patchIDs_ = allWalls.sortedToc();
+
+        if (!patchSet.empty())
         {
-            if (isA<wallPolyPatch>(pbm[patchi]))
+            WarningInFunction
+                << "Requested wall shear stress on ("
+                << patchSet.size() << ") non-wall patches:" << nl;
+
+            for (const label patchi : patchSet.sortedToc())
             {
-                filteredPatchSet.insert(patchi);
-                Info<< "        " << pbm[patchi].name() << endl;
-            }
-            else
-            {
-                WarningInFunction
-                    << "Requested wall shear stress on non-wall boundary "
-                    << "type patch: " << pbm[patchi].name() << endl;
+                Info<< "        " << pbm[patchi].name() << nl;
             }
+            Info<< nl;
         }
 
-        Info<< endl;
+        Info<< "    processing (" << patchIDs_.size()
+            << ") wall patches:" << nl;
 
-        patchSet_ = filteredPatchSet;
+        for (const label patchi : patchIDs_)
+        {
+            Info<< "        " << pbm[patchi].name() << nl;
+        }
+        Info<< endl;
     }
 
     return true;
@@ -234,14 +237,13 @@ bool Foam::functionObjects::wallShearStress::write()
 
     const fvPatchList& patches = mesh_.boundary();
 
-    for (const label patchi : patchSet_)
+    for (const label patchi : patchIDs_)
     {
         const fvPatch& pp = patches[patchi];
 
         const vectorField& ssp = wallShearStress.boundaryField()[patchi];
 
-        const vector minSsp = gMin(ssp);
-        const vector maxSsp = gMax(ssp);
+        const MinMax<vector> limits = gMinMax(ssp);
 
         if (UPstream::master())
         {
@@ -249,13 +251,13 @@ bool Foam::functionObjects::wallShearStress::write()
 
             file()
                 << token::TAB << pp.name()
-                << token::TAB << minSsp
-                << token::TAB << maxSsp
+                << token::TAB << limits.min()
+                << token::TAB << limits.max()
                 << endl;
         }
 
         Log << "    min/max(" << pp.name() << ") = "
-            << minSsp << ", " << maxSsp << endl;
+            << limits.min() << ", " << limits.max() << endl;
     }
 
     return true;
diff --git a/src/functionObjects/field/wallShearStress/wallShearStress.H b/src/functionObjects/field/wallShearStress/wallShearStress.H
index c0908f1276..4d1c95d614 100644
--- a/src/functionObjects/field/wallShearStress/wallShearStress.H
+++ b/src/functionObjects/field/wallShearStress/wallShearStress.H
@@ -75,7 +75,7 @@ Usage
       Property     | Description                           | Type | Req'd | Dflt
       type         | Type name: wallShearStress            | word |  yes  | -
       libs         | Library name: fieldFunctionObjects    | word |  yes  | -
-      patches    | Names of operand patches   | wordList | no | all wall patches
+      patches    | Names of operand patches   | wordRes | no | all wall patches
       writeFields | Flag to write shear stress field       | bool |  no   | true
     \endtable
 
@@ -111,7 +111,6 @@ SourceFiles
 #include "fvMeshFunctionObject.H"
 #include "writeFile.H"
 #include "volFieldsFwd.H"
-#include "HashSet.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -139,8 +138,8 @@ protected:
 
     // Protected Data
 
-        //- Optional list of patches to process
-        labelHashSet patchSet_;
+        //- Wall patches to process (optionally filtered by name)
+        labelList patchIDs_;
 
 
     // Protected Member Functions
diff --git a/src/functionObjects/forces/forces/forces.C b/src/functionObjects/forces/forces/forces.C
index 5f8c48479c..e1ff6bba1d 100644
--- a/src/functionObjects/forces/forces/forces.C
+++ b/src/functionObjects/forces/forces/forces.C
@@ -206,7 +206,7 @@ void Foam::functionObjects::forces::reset()
     else
     {
         constexpr bool updateAccessTime = false;
-        for (const label patchi : patchSet_)
+        for (const label patchi : patchIDs_)
         {
             force.boundaryFieldRef(updateAccessTime)[patchi] = Zero;
             moment.boundaryFieldRef(updateAccessTime)[patchi] = Zero;
@@ -545,7 +545,6 @@ Foam::functionObjects::forces::forces
     forceFilePtr_(),
     momentFilePtr_(),
     coordSysPtr_(nullptr),
-    patchSet_(),
     rhoRef_(VGREAT),
     pRef_(0),
     pName_("p"),
@@ -585,7 +584,6 @@ Foam::functionObjects::forces::forces
     forceFilePtr_(),
     momentFilePtr_(),
     coordSysPtr_(nullptr),
-    patchSet_(),
     rhoRef_(VGREAT),
     pRef_(0),
     pName_("p"),
@@ -610,6 +608,8 @@ Foam::functionObjects::forces::forces
 
 bool Foam::functionObjects::forces::read(const dictionary& dict)
 {
+    const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
+
     if (!fvMeshFunctionObject::read(dict) || !writeFile::read(dict))
     {
         return false;
@@ -617,13 +617,10 @@ bool Foam::functionObjects::forces::read(const dictionary& dict)
 
     initialised_ = false;
 
-    Info<< type() << " " << name() << ":" << endl;
+    Info<< type() << ' ' << name() << ':' << endl;
 
-    patchSet_ =
-        mesh_.boundaryMesh().patchSet
-        (
-            dict.get<wordRes>("patches")
-        );
+    // Can also use pbm.indices(), but no warnings...
+    patchIDs_ = pbm.patchSet(dict.get<wordRes>("patches")).sortedToc();
 
     dict.readIfPresent("directForceDensity", directForceDensity_);
     if (directForceDensity_)
@@ -702,7 +699,7 @@ void Foam::functionObjects::forces::calcForcesMoments()
         const auto& magSfb = mesh_.magSf().boundaryField();
         const auto& Cb = mesh_.C().boundaryField();
 
-        for (const label patchi : patchSet_)
+        for (const label patchi : patchIDs_)
         {
             const vectorField Md(Cb[patchi] - origin);
 
@@ -738,7 +735,7 @@ void Foam::functionObjects::forces::calcForcesMoments()
         const scalar rhoRef = rho(p);
         const scalar pRef = pRef_/rhoRef;
 
-        for (const label patchi : patchSet_)
+        for (const label patchi : patchIDs_)
         {
             const vectorField Md(Cb[patchi] - origin);
 
diff --git a/src/functionObjects/forces/forces/forces.H b/src/functionObjects/forces/forces/forces.H
index 0611bf33d6..cd96a754e0 100644
--- a/src/functionObjects/forces/forces/forces.H
+++ b/src/functionObjects/forces/forces/forces.H
@@ -231,8 +231,8 @@ protected:
             //- Coordinate system used when evaluating forces and moments
             autoPtr<coordinateSystem> coordSysPtr_;
 
-            //- Names of operand patches
-            labelHashSet patchSet_;
+            //- Selected operand patches
+            labelList patchIDs_;
 
             //- Reference density needed for incompressible calculations
             scalar rhoRef_;
diff --git a/src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/patchInjection/patchInjection.C b/src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/patchInjection/patchInjection.C
index ea0fbb347d..24acdaf03b 100644
--- a/src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/patchInjection/patchInjection.C
+++ b/src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/patchInjection/patchInjection.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015-2017 OpenFOAM Foundation
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -55,40 +55,36 @@ patchInjection::patchInjection
     deltaStable_(coeffDict_.getOrDefault<scalar>("deltaStable", 0))
 {
     const polyBoundaryMesh& pbm = film.regionMesh().boundaryMesh();
-    patchIDs_.setSize
+
+    const label nPatches
     (
-        pbm.size() - film.regionMesh().globalData().processorPatches().size()
+        pbm.size()
+      - film.regionMesh().globalData().processorPatches().size()
     );
 
     wordRes patchNames;
     if (coeffDict_.readIfPresent("patches", patchNames))
     {
-        const labelHashSet patchSet = pbm.patchSet(patchNames);
+        // Can also use pbm.indices(), but no warnings...
+        patchIDs_ = pbm.patchSet(patchNames).sortedToc();
 
         Info<< "        applying to patches:" << nl;
 
-        label pidi = 0;
-        for (const label patchi : patchSet)
+        for (const label patchi : patchIDs_)
         {
-            patchIDs_[pidi++] = patchi;
             Info<< "            " << pbm[patchi].name() << endl;
         }
-        patchIDs_.setSize(pidi);
-        patchInjectedMasses_.setSize(pidi, 0);
     }
     else
     {
         Info<< "            applying to all patches" << endl;
 
-        forAll(patchIDs_, patchi)
-        {
-            patchIDs_[patchi] = patchi;
-        }
-
-        patchInjectedMasses_.setSize(patchIDs_.size(), 0);
+        patchIDs_ = identity(nPatches);
     }
 
-    if (!patchIDs_.size())
+    patchInjectedMasses_.resize(patchIDs_.size(), Zero);
+
+    if (patchIDs_.empty())
     {
         FatalErrorInFunction
             << "No patches selected"
@@ -97,12 +93,6 @@ patchInjection::patchInjection
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-patchInjection::~patchInjection()
-{}
-
-
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 void patchInjection::correct
@@ -113,7 +103,7 @@ void patchInjection::correct
 )
 {
     // Do not correct if no patches selected
-    if (!patchIDs_.size()) return;
+    if (patchIDs_.empty()) return;
 
     const scalarField& delta = film().delta();
     const scalarField& rho = film().rho();
@@ -130,10 +120,8 @@ void patchInjection::correct
         // Accumulate the total mass removed from patch
         scalar dMassPatch = 0;
 
-        forAll(faceCells, fci)
+        for (const label celli : faceCells)
         {
-            label celli = faceCells[fci];
-
             scalar ddelta = max(0.0, delta[celli] - deltaStable_);
             scalar dMass = ddelta*rho[celli]*magSf[celli];
             massToInject[celli] += dMass;
@@ -154,7 +142,7 @@ void patchInjection::correct
             getModelProperty<scalarField>
             (
                 "patchInjectedMasses",
-                scalarField(patchInjectedMasses_.size(), 0)
+                scalarField(patchInjectedMasses_.size(), Zero)
             )
         );
 
diff --git a/src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/patchInjection/patchInjection.H b/src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/patchInjection/patchInjection.H
index c89c9de48a..3eea258e45 100644
--- a/src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/patchInjection/patchInjection.H
+++ b/src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/patchInjection/patchInjection.H
@@ -35,8 +35,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef patchInjection_H
-#define patchInjection_H
+#ifndef Foam_patchInjection_H
+#define Foam_patchInjection_H
 
 #include "injectionModel.H"
 
@@ -57,15 +57,6 @@ class patchInjection
 :
     public injectionModel
 {
-    // Private member functions
-
-        //- No copy construct
-        patchInjection(const patchInjection&) = delete;
-
-        //- No copy assignment
-        void operator=(const patchInjection&) = delete;
-
-
 protected:
 
         //- Stable film thickness - mass only removed if thickness exceeds
@@ -81,6 +72,13 @@ protected:
 
 public:
 
+    //- No copy construct
+    patchInjection(const patchInjection&) = delete;
+
+    //- No copy assignment
+    void operator=(const patchInjection&) = delete;
+
+
     //- Runtime type information
     TypeName("patchInjection");
 
@@ -92,7 +90,7 @@ public:
 
 
     //- Destructor
-    virtual ~patchInjection();
+    virtual ~patchInjection() = default;
 
 
     // Member Functions