From d38de84d212ae30fde2bd9c74d88bac65eb71ba8 Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Mon, 28 Mar 2022 16:30:51 +0200
Subject: [PATCH] ENH: bundle Pstream:: AllGather methods

- bundles frequently used 'gather/scatter' patterns more consistently.

  - combineAllGather     -> combineGather + broadcast
  - listCombineAllGather -> listCombineGather + broadcast
  - mapCombineAllGather  -> mapCombineGather + broadcast
  - allGatherList        -> gatherList + scatterList
  - reduce               -> gather + broadcast (ie, allreduce)

- The allGatherList currently wraps gatherList/scatterList, but may be
  replaced with a different algorithm in the future.

STYLE: PstreamCombineReduceOps.H is mostly unneeded now
---
 applications/test/Tuple2/Test-Tuple2.C        |  10 +-
 applications/test/treeComms/Make/files        |   3 +
 applications/test/treeComms/Make/options      |   2 +
 applications/test/treeComms/Test-treeComms.C  | 241 ++++++++++++++
 .../extrudeToRegionMesh/extrudeToRegionMesh.C |  18 +-
 .../DelaunayMesh/DistributedDelaunayMesh.C    |   5 +-
 .../backgroundMeshDecomposition.C             |   5 +-
 .../conformalVoronoiMeshCalcDualMesh.C        |   5 +-
 .../conformalVoronoiMeshConformToSurface.C    |  18 +-
 .../initialPointsMethod/pointFile/pointFile.C |   7 +-
 .../manipulation/checkMesh/checkTopology.C    |   7 +-
 .../manipulation/createPatch/createPatch.C    |   5 +-
 .../splitMeshRegions/splitMeshRegions.C       |   3 +-
 .../foamFormatConvert/foamFormatConvert.C     |   8 +-
 .../parLagrangianRedistributor.C              |   2 +-
 .../redistributePar/redistributePar.C         |   5 +-
 .../foamToEnsight/findCloudFields.H           |   3 +-
 .../foamToEnsight/foamToEnsight.C             |   2 +-
 .../particleTracks/particleTracks.C           |   3 +-
 .../postChannel/channelIndexTemplates.C       |   3 +-
 .../preProcessing/setFields/setFields.C       |   3 +-
 .../viewFactorsGen/viewFactorsGen.C           |  11 +-
 .../surfaceRedistributePar.C                  |   3 +-
 src/OpenFOAM/db/IOstreams/Pstreams/Pstream.H  | 314 ++++++++++++------
 .../IOstreams/Pstreams/PstreamCombineGather.C | 230 +++++++++----
 .../Pstreams/PstreamCombineReduceOps.H        |  12 +-
 .../db/IOstreams/Pstreams/PstreamGather.C     |  26 +-
 .../db/IOstreams/Pstreams/PstreamGatherList.C |  53 ++-
 src/OpenFOAM/db/IOstreams/Pstreams/UPstream.H |  12 +-
 src/OpenFOAM/fields/Fields/Field/FieldOps.C   |   9 +-
 .../hostCollatedFileOperation.C               |   5 +-
 .../procFacesGAMGProcAgglomeration.C          |   3 +-
 .../ProcessorTopology/ProcessorTopology.C     |   4 +-
 .../polyMesh/globalMeshData/globalMeshData.C  |   5 +-
 .../mapDistributeBaseTemplates.C              |   1 -
 .../mapDistribute/mapDistributeTemplates.C    |   1 -
 .../polyBoundaryMesh/polyBoundaryMesh.C       |   6 +-
 .../meshes/polyMesh/syncTools/syncTools.H     |   4 +-
 .../meshes/polyMesh/zones/ZoneMesh/ZoneMesh.C |   6 +-
 .../powerLawLopesdaCosta.C                    |   8 +-
 .../freeSurfacePointDisplacement.C            |   4 +-
 .../pointHistory/pointHistory.C               |   3 +-
 .../extrudePatchMesh/extrudePatchMesh.C       |   3 +-
 .../fvMeshDistribute/fvMeshDistribute.C       |   1 -
 .../fvMeshDistributeTemplates.C               |   3 +-
 src/dynamicMesh/fvMeshSubset/fvMeshSubset.C   |   8 +-
 .../polyTopoChange/hexRef8/hexRef8.C          |   6 +-
 .../faGlobalMeshData/faGlobalMeshData.C       |   3 +-
 .../faMesh/faMeshDemandDrivenData.C           |   8 +-
 .../turbulentDFSEMInletFvPatchVectorField.C   |  17 +-
 src/functionObjects/field/Curle/Curle.C       |   3 +-
 .../columnAverage/columnAverageTemplates.C    |   6 +-
 .../extractEulerianParticles.C                |   7 +-
 .../field/fieldMinMax/fieldMinMaxTemplates.C  |  20 +-
 .../fieldValue/fieldValueTemplates.C          |   5 +-
 .../regionSizeDistribution.C                  |   6 +-
 src/functionObjects/forces/forces/forces.C    |   6 +-
 .../forces/propellerInfo/propellerInfo.C      |   6 +-
 .../utilities/areaWrite/areaWrite.C           |   7 +-
 src/lagrangian/basic/Cloud/Cloud.C            |   2 +-
 src/lagrangian/basic/Cloud/CloudIO.C          |   3 +-
 .../basic/InteractionLists/InteractionLists.C |   6 +-
 .../CellZoneInjection/CellZoneInjection.C     |   3 +-
 .../PatchInjection/patchInjectionBase.C       |   9 +-
 .../movement/lumpedPointMovement.C            |   9 +-
 .../meshRefinement/meshRefinement.C           |  32 +-
 .../meshRefinement/meshRefinementBaffles.C    |  15 +-
 .../refinementSurfaces/surfaceZonesInfo.C     |   8 +-
 .../snappyHexMeshDriver/snappyRefineDriver.C  |   9 +-
 .../advancingFrontAMIParallelOps.C            |   6 +-
 .../nearestFaceAMI/nearestFaceAMI.C           |   5 +-
 .../cyclicPeriodicAMIPolyPatch.C              |  12 +-
 .../algorithms/PointEdgeWave/PointEdgeWave.C  |   1 -
 .../mappedPolyPatch/mappedPatchBase.C         |  21 +-
 src/meshTools/processorLOD/box/box.C          |   3 +-
 src/meshTools/regionSplit2D/regionSplit2D.C   |   3 +-
 .../cellSources/nearestToCell/nearestToCell.C |   3 +-
 .../faceSources/regionToFace/regionToFace.C   |   4 +-
 .../planeToFaceZone/planeToFaceZone.C         |  12 +-
 .../nearestToPoint/nearestToPoint.C           |   3 +-
 .../sensitivitySurfacePointsIncompressible.C  |   5 +-
 .../NURBS3DVolume/NURBS3DVolume.C             |   6 +-
 .../cellVolumeWeightCellCellStencil.C         |   6 +-
 .../inverseDistanceCellCellStencil.C          |  13 +-
 .../trackingInverseDistanceCellCellStencil.C  |   8 +-
 src/overset/regionsToCell/findRefCells.C      |   6 +-
 .../multiLevelDecomp/multiLevelDecomp.C       |  11 +-
 .../distributedTriSurfaceMesh.C               |  22 +-
 .../faReconstruct/faMeshReconstructor.C       |   3 +-
 .../singleLayerRegion/singleLayerRegion.C     |   3 +-
 .../meshToMesh/meshToMeshParallelOps.C        |   7 +-
 src/sampling/probes/patchProbes.C             |   5 +-
 src/sampling/probes/patchProbesTemplates.C    |   6 +-
 src/sampling/probes/probesTemplates.C         |   6 +-
 src/sampling/sampledSet/cloud/cloudSet.C      |   7 +-
 .../sampledSet/patchCloud/patchCloudSet.C     |   5 +-
 .../sampledSet/patchSeed/patchSeedSet.C       |  10 +-
 .../sampledSet/sampledSets/sampledSets.C      |   7 +-
 .../sampledSet/shortestPath/shortestPathSet.C |   2 +-
 .../sampledMeshedSurface.C                    |   3 +-
 .../sampledSurfaces/sampledSurfaces.C         |   7 +-
 .../surface/isoSurface/isoSurfacePoint.C      |   5 +-
 .../TDACChemistryModel/TDACChemistryModel.C   |   3 +-
 .../solarLoad/faceReflecting/faceReflecting.C |  10 +-
 .../radiationModels/viewFactor/viewFactor.C   |  12 +-
 105 files changed, 935 insertions(+), 605 deletions(-)
 create mode 100644 applications/test/treeComms/Make/files
 create mode 100644 applications/test/treeComms/Make/options
 create mode 100644 applications/test/treeComms/Test-treeComms.C

diff --git a/applications/test/Tuple2/Test-Tuple2.C b/applications/test/Tuple2/Test-Tuple2.C
index 13375a0f71..95df0e5a42 100644
--- a/applications/test/Tuple2/Test-Tuple2.C
+++ b/applications/test/Tuple2/Test-Tuple2.C
@@ -39,7 +39,7 @@ Description
 #include "List.H"
 #include "ListOps.H"
 #include "ops.H"
-#include "PstreamCombineReduceOps.H"
+#include "Pstream.H"
 #include <functional>
 
 using namespace Foam;
@@ -138,8 +138,8 @@ int main()
             maxFirstEqOp<label>()(maxIndexed, item);
         }
 
-        Foam::combineReduce(minIndexed, minFirstEqOp<label>());
-        Foam::combineReduce(maxIndexed, maxFirstEqOp<label>());
+        Pstream::combineAllGather(minIndexed, minFirstEqOp<label>());
+        Pstream::combineAllGather(maxIndexed, maxFirstEqOp<label>());
 
         Info<< "Min indexed: " << minIndexed << nl
             << "Max indexed: " << maxIndexed << nl;
@@ -156,8 +156,8 @@ int main()
             maxIndexed = maxFirstOp<label>()(maxIndexed, item);
         }
 
-        Foam::combineReduce(minIndexed, minFirstEqOp<label>());
-        Foam::combineReduce(maxIndexed, maxFirstEqOp<label>());
+        Pstream::combineAllGather(minIndexed, minFirstEqOp<label>());
+        Pstream::combineAllGather(maxIndexed, maxFirstEqOp<label>());
 
         Info<< "Min indexed: " << minIndexed << nl
             << "Max indexed: " << maxIndexed << nl;
diff --git a/applications/test/treeComms/Make/files b/applications/test/treeComms/Make/files
new file mode 100644
index 0000000000..4cf2c50a9b
--- /dev/null
+++ b/applications/test/treeComms/Make/files
@@ -0,0 +1,3 @@
+Test-treeComms.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-treeComms
diff --git a/applications/test/treeComms/Make/options b/applications/test/treeComms/Make/options
new file mode 100644
index 0000000000..18e6fe47af
--- /dev/null
+++ b/applications/test/treeComms/Make/options
@@ -0,0 +1,2 @@
+/* EXE_INC = */
+/* EXE_LIBS = */
diff --git a/applications/test/treeComms/Test-treeComms.C b/applications/test/treeComms/Test-treeComms.C
new file mode 100644
index 0000000000..fe226b828a
--- /dev/null
+++ b/applications/test/treeComms/Test-treeComms.C
@@ -0,0 +1,241 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Description
+    Print/test tree communication patterns
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "List.H"
+#include "Time.H"
+#include "vector.H"
+#include "Fstream.H"
+#include "IPstream.H"
+#include "OPstream.H"
+#include "IOstreams.H"
+
+using namespace Foam;
+
+void printConnection(Ostream& os, const label proci, const labelUList& below)
+{
+    for (const label connectProci : below)
+    {
+        os << indent << proci << " -- " << connectProci << nl;
+    }
+}
+
+
+// The number of receives - as per gatherList (v2112)
+void printRecvCount_gatherList
+(
+    const List<UPstream::commsStruct>& comms,
+    const label comm = UPstream::worldComm
+)
+{
+    const label np = UPstream::nProcs(comm);
+    labelList nMesg;
+    labelList nRecv;
+
+    if (UPstream::parRun() && np > 1 && Pstream::master(comm))
+    {
+        nMesg.resize(np, Zero);
+        nRecv.resize(np, Zero);
+
+        forAll(comms, proci)
+        {
+            nMesg[proci] += comms[proci].below().size();
+
+            // Receive from my downstairs neighbours
+            for (const label belowID : comms[proci].below())
+            {
+                const labelList& belowLeaves = comms[belowID].allBelow();
+                nRecv[proci] += (belowLeaves.size() + 1);
+            }
+        }
+    }
+
+    label nTotalMesg = sum(nMesg);
+    label nTotalRecv = sum(nRecv);
+
+    Info<< "gatherList communication:" << nl
+        << "messages = " << nTotalMesg << " -> "
+        << flatOutput(nMesg) << nl
+        << "gathers = " << nTotalRecv << " -> "
+        << flatOutput(nRecv) << nl;
+}
+
+
+// The number of sends - as per scatterList (v2112)
+void printSendCount_scatterList
+(
+    const List<UPstream::commsStruct>& comms,
+    const label comm = UPstream::worldComm
+)
+{
+    const label np = UPstream::nProcs(comm);
+    labelList nMesg;
+    labelList nSend;
+
+    if (UPstream::parRun() && np > 1 && Pstream::master(comm))
+    {
+        nMesg.resize(np, Zero);
+        nSend.resize(np, Zero);
+
+        forAll(comms, proci)
+        {
+            nMesg[proci] += comms[proci].below().size();
+
+            // Send to my downstairs neighbours
+            for (const label belowID : comms[proci].below())
+            {
+                const labelList& notBelowLeaves = comms[belowID].allNotBelow();
+                nSend[proci] += notBelowLeaves.size();
+            }
+        }
+    }
+
+    label nTotalMesg = sum(nMesg);
+    label nTotalSend = sum(nSend);
+
+    Info<< "scatterList communication:" << nl
+        << "messages = " << nTotalMesg << " -> "
+        << flatOutput(nMesg) << nl
+        << "scatters = " << nTotalSend << " -> "
+        << flatOutput(nSend) << nl;
+}
+
+
+// Transmission widths (contiguous data)
+void printWidths
+(
+    const List<UPstream::commsStruct>& comms,
+    const label comm = UPstream::worldComm
+)
+{
+    const label np = UPstream::nProcs(comm);
+    labelList maxBelow;
+    labelList maxNotBelow;
+
+    if (UPstream::parRun() && np > 1 && Pstream::master(comm))
+    {
+        maxBelow.resize(np, Zero);
+        maxNotBelow.resize(np, Zero);
+
+        forAll(comms, proci)
+        {
+            label& max0 = maxBelow[proci];
+            label& max1 = maxNotBelow[proci];
+
+            for (const label belowID : comms[proci].below())
+            {
+                // Receive from my downstairs neighbours
+                const labelList& belowLeaves = comms[belowID].allBelow();
+
+                // Send to my downstairs neighbours
+                const labelList& notBelowLeaves = comms[belowID].allNotBelow();
+
+                max0 = max(max0, belowLeaves.size());
+                max1 = max(max1, notBelowLeaves.size());
+            }
+        }
+    }
+
+    Info<< "buffer width:" << nl
+        << "gathers  -> " << flatOutput(maxBelow) << nl
+        << "scatters -> " << flatOutput(maxNotBelow) << nl;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::noCheckProcessorDirectories();
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+
+    if (!Pstream::parRun())
+    {
+        FatalErrorInFunction
+            << "Please run in parallel" << exit(FatalError);
+    }
+    const auto& comms = UPstream::treeCommunication();
+
+    printRecvCount_gatherList(comms);
+    printSendCount_scatterList(comms);
+    printWidths(comms);
+
+    // My communication order
+    const UPstream::commsStruct& myComm = comms[UPstream::myProcNo()];
+
+    // Info<< "allComms: " << comms << nl;
+
+    if (Pstream::master())
+    {
+        OFstream os("treeComm.dot");
+
+        os << "// tree communication" << nl << nl;
+        os.beginBlock("graph");
+
+        printConnection(os, 0, myComm.below());
+        // Pout<< flatOutput(myComm.allBelow()) << nl;
+
+        for (const int proci : Pstream::subProcs())
+        {
+            IPstream fromProc(Pstream::commsTypes::scheduled, proci);
+            labelList below(fromProc);
+
+            printConnection(os, proci, below);
+        }
+
+        os.endBlock();
+
+        Info<< "Wrote processorTopology graph: "
+            << runTime.relativePath(os.name()) << nl;
+
+        Info<< nl
+            << "Use dot (or other graphviz tools)" << nl;
+    }
+    else
+    {
+        OPstream toMaster
+        (
+             Pstream::commsTypes::scheduled,
+             Pstream::masterNo()
+        );
+
+        toMaster << myComm.below();
+        // Pout<< flatOutput(myComm.allBelow()) << nl;
+    }
+
+    Info<< "\nEnd\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMesh.C b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMesh.C
index 0eee01dc5a..b2acecd41e 100644
--- a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMesh.C
+++ b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMesh.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2020 OpenCFD Ltd.
+    Copyright (C) 2015-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -661,11 +661,8 @@ void countExtrudePatches
     }
     // Synchronise decision. Actual numbers are not important, just make
     // sure that they're > 0 on all processors.
-    Pstream::listCombineGather(zoneSidePatch, plusEqOp<label>());
-    Pstream::listCombineGather(zoneZonePatch, plusEqOp<label>());
-
-    Pstream::broadcast(zoneSidePatch);
-    Pstream::broadcast(zoneZonePatch);
+    Pstream::listCombineAllGather(zoneSidePatch, plusEqOp<label>());
+    Pstream::listCombineAllGather(zoneZonePatch, plusEqOp<label>());
 }
 
 
@@ -1465,8 +1462,7 @@ int main(int argc, char *argv[])
     {
         List<wordList> allNames(Pstream::nProcs());
         allNames[Pstream::myProcNo()] = mesh.boundaryMesh().names();
-        Pstream::gatherList(allNames);
-        Pstream::scatterList(allNames);
+        Pstream::allGatherList(allNames);
 
         FatalErrorInFunction
             << "Patches are not synchronised on all processors."
@@ -1858,8 +1854,7 @@ int main(int argc, char *argv[])
     const primitiveFacePatch extrudePatch(std::move(zoneFaces), mesh.points());
 
 
-    Pstream::listCombineGather(isInternal, orEqOp<bool>());
-    Pstream::broadcast(isInternal);
+    Pstream::listCombineAllGather(isInternal, orEqOp<bool>());
 
     // Check zone either all internal or all external faces
     checkZoneInside(mesh, zoneNames, zoneID, extrudeMeshFaces, isInternal);
@@ -2320,8 +2315,7 @@ int main(int argc, char *argv[])
         }
 
         // Reduce
-        Pstream::mapCombineGather(globalSum, plusEqOp<point>());
-        Pstream::broadcast(globalSum);
+        Pstream::mapCombineAllGather(globalSum, plusEqOp<point>());
 
         forAll(localToGlobalRegion, localRegionI)
         {
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/DelaunayMesh/DistributedDelaunayMesh.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/DelaunayMesh/DistributedDelaunayMesh.C
index 3a93a47583..b51546a528 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/DelaunayMesh/DistributedDelaunayMesh.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/DelaunayMesh/DistributedDelaunayMesh.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2016 OpenFOAM Foundation
-    Copyright (C) 2020-2021 OpenCFD Ltd.
+    Copyright (C) 2020-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -157,8 +157,7 @@ bool Foam::DistributedDelaunayMesh<Triangulation>::distributeBoundBoxes
     // Give the bounds of every processor to every other processor
     allBackgroundMeshBounds_()[Pstream::myProcNo()] = bb;
 
-    Pstream::gatherList(allBackgroundMeshBounds_());
-    Pstream::scatterList(allBackgroundMeshBounds_());
+    Pstream::allGatherList(allBackgroundMeshBounds_());
 
     return true;
 }
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
index bc44ceab9c..afb0baa2cf 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2017-2020 OpenCFD Ltd.
+    Copyright (C) 2017-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -713,8 +713,7 @@ void Foam::backgroundMeshDecomposition::buildPatchAndTree()
     // Give the bounds of every processor to every other processor
     allBackgroundMeshBounds_[Pstream::myProcNo()] = overallBb;
 
-    Pstream::gatherList(allBackgroundMeshBounds_);
-    Pstream::scatterList(allBackgroundMeshBounds_);
+    Pstream::allGatherList(allBackgroundMeshBounds_);
 
     // find global bounding box
     globalBackgroundBounds_ = treeBoundBox(boundBox::invertedBox);
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
index 14a2f32fd5..ad2ce6d098 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2016 OpenFOAM Foundation
-    Copyright (C) 2018-2021 OpenCFD Ltd.
+    Copyright (C) 2018-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -1564,8 +1564,7 @@ Foam::label Foam::conformalVoronoiMesh::createPatchInfo
         }
 
         // Because the previous test was insufficient, combine the lists.
-        Pstream::gatherList(procUsedList);
-        Pstream::scatterList(procUsedList);
+        Pstream::allGatherList(procUsedList);
 
         forAll(procUsedList, proci)
         {
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
index 46166f4f1a..0fe385e0bd 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2016 OpenFOAM Foundation
-    Copyright (C) 2020-2021 OpenCFD Ltd.
+    Copyright (C) 2020-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -692,11 +692,8 @@ Foam::label Foam::conformalVoronoiMesh::synchroniseSurfaceTrees
     );
 
     List<pointIndexHitAndFeatureDynList> procSurfLocations(Pstream::nProcs());
-
     procSurfLocations[Pstream::myProcNo()] = surfaceHits;
-
-    Pstream::gatherList(procSurfLocations);
-    Pstream::scatterList(procSurfLocations);
+    Pstream::allGatherList(procSurfLocations);
 
     List<labelHashSet> hits(Pstream::nProcs());
 
@@ -732,8 +729,7 @@ Foam::label Foam::conformalVoronoiMesh::synchroniseSurfaceTrees
         }
     }
 
-    Pstream::listCombineGather(hits, plusEqOp<labelHashSet>());
-    Pstream::broadcast(hits);
+    Pstream::listCombineAllGather(hits, plusEqOp<labelHashSet>());
 
     forAll(surfaceHits, eI)
     {
@@ -780,11 +776,8 @@ Foam::label Foam::conformalVoronoiMesh::synchroniseEdgeTrees
     );
 
     List<pointIndexHitAndFeatureDynList> procEdgeLocations(Pstream::nProcs());
-
     procEdgeLocations[Pstream::myProcNo()] = featureEdgeHits;
-
-    Pstream::gatherList(procEdgeLocations);
-    Pstream::scatterList(procEdgeLocations);
+    Pstream::allGatherList(procEdgeLocations);
 
     List<labelHashSet> hits(Pstream::nProcs());
 
@@ -823,8 +816,7 @@ Foam::label Foam::conformalVoronoiMesh::synchroniseEdgeTrees
         }
     }
 
-    Pstream::listCombineGather(hits, plusEqOp<labelHashSet>());
-    Pstream::broadcast(hits);
+    Pstream::listCombineAllGather(hits, plusEqOp<labelHashSet>());
 
     forAll(featureEdgeHits, eI)
     {
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/initialPointsMethod/pointFile/pointFile.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/initialPointsMethod/pointFile/pointFile.C
index 683154f913..da33aea938 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/initialPointsMethod/pointFile/pointFile.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/initialPointsMethod/pointFile/pointFile.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2016 OpenFOAM Foundation
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -142,10 +142,7 @@ Foam::List<Vb::Point> Foam::pointFile::initialPoints() const
             List<boolList> allProcPt(Pstream::nProcs());
 
             allProcPt[Pstream::myProcNo()] = procPt;
-
-            Pstream::gatherList(allProcPt);
-
-            Pstream::scatterList(allProcPt);
+            Pstream::allGatherList(allProcPt);
 
             forAll(procPt, ptI)
             {
diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
index 6edf38ea6f..d1b8240b5e 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
@@ -474,8 +474,11 @@ Foam::label Foam::checkTopology
                     }
                 }
 
-                Pstream::listCombineGather(regionDisconnected, andEqOp<bool>());
-                Pstream::broadcast(regionDisconnected);
+                Pstream::listCombineAllGather
+                (
+                    regionDisconnected,
+                    andEqOp<bool>()
+                );
             }
 
 
diff --git a/applications/utilities/mesh/manipulation/createPatch/createPatch.C b/applications/utilities/mesh/manipulation/createPatch/createPatch.C
index 5c112c81f8..9a8dcd47db 100644
--- a/applications/utilities/mesh/manipulation/createPatch/createPatch.C
+++ b/applications/utilities/mesh/manipulation/createPatch/createPatch.C
@@ -713,9 +713,8 @@ void syncPoints
             sharedPts[pd.sharedPointAddr()[i]] = points[meshPointi];
         }
 
-        // Combine on master.
-        Pstream::listCombineGather(sharedPts, cop);
-        Pstream::broadcast(sharedPts);
+        // Combine - globally consistent
+        Pstream::listCombineAllGather(sharedPts, cop);
 
         // Now we will all have the same information. Merge it back with
         // my local information.
diff --git a/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C b/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
index b1d4691206..e81301ab51 100644
--- a/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
+++ b/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
@@ -1111,8 +1111,7 @@ label findCorrespondingRegion
         }
     }
 
-    Pstream::listCombineGather(cellsInZone, plusEqOp<label>());
-    Pstream::broadcast(cellsInZone);
+    Pstream::listCombineAllGather(cellsInZone, plusEqOp<label>());
 
     // Pick region with largest overlap of zoneI
     label regionI = findMax(cellsInZone);
diff --git a/applications/utilities/miscellaneous/foamFormatConvert/foamFormatConvert.C b/applications/utilities/miscellaneous/foamFormatConvert/foamFormatConvert.C
index f21eac25a2..518e2510c6 100644
--- a/applications/utilities/miscellaneous/foamFormatConvert/foamFormatConvert.C
+++ b/applications/utilities/miscellaneous/foamFormatConvert/foamFormatConvert.C
@@ -223,7 +223,7 @@ bool writeOptionalMeshObject
 
     // Make sure all know if there is a valid class name
     wordList classNames(1, io.headerClassName());
-    combineReduce(classNames, uniqueEqOp<word>());
+    Pstream::combineAllGather(classNames, uniqueEqOp<word>());
 
     // Check for correct type
     if (classNames[0] == T::typeName)
@@ -422,7 +422,7 @@ int main(int argc, char *argv[])
             )
         );
 
-        combineReduce(lagrangianDirs, uniqueEqOp<fileName>());
+        Pstream::combineAllGather(lagrangianDirs, uniqueEqOp<fileName>());
 
         if (!lagrangianDirs.empty())
         {
@@ -459,7 +459,7 @@ int main(int argc, char *argv[])
                 )
             );
 
-            combineReduce(cloudDirs, uniqueEqOp<fileName>());
+            Pstream::combineAllGather(cloudDirs, uniqueEqOp<fileName>());
 
             forAll(cloudDirs, i)
             {
@@ -485,7 +485,7 @@ int main(int argc, char *argv[])
                 );
 
                 // Combine with all other cloud objects
-                combineReduce(cloudFields, uniqueEqOp<word>());
+                Pstream::combineAllGather(cloudFields, uniqueEqOp<word>());
 
                 for (const word& name : cloudFields)
                 {
diff --git a/applications/utilities/parallelProcessing/redistributePar/parLagrangianRedistributor.C b/applications/utilities/parallelProcessing/redistributePar/parLagrangianRedistributor.C
index ee52b18bbb..b5dbf55215 100644
--- a/applications/utilities/parallelProcessing/redistributePar/parLagrangianRedistributor.C
+++ b/applications/utilities/parallelProcessing/redistributePar/parLagrangianRedistributor.C
@@ -286,7 +286,7 @@ Foam::parLagrangianRedistributor::redistributeLagrangianPositions
         nsTransPs[sendProcI] = subMap[sendProcI].size();
     }
     // Send sizes across. Note: blocks.
-    combineReduce(sizes, Pstream::listEq());
+    Pstream::combineAllGather(sizes, Pstream::listEq());
 
     labelListList constructMap(Pstream::nProcs());
     label constructSize = 0;
diff --git a/applications/utilities/parallelProcessing/redistributePar/redistributePar.C b/applications/utilities/parallelProcessing/redistributePar/redistributePar.C
index 8c49938be0..dfc2572961 100644
--- a/applications/utilities/parallelProcessing/redistributePar/redistributePar.C
+++ b/applications/utilities/parallelProcessing/redistributePar/redistributePar.C
@@ -2649,7 +2649,7 @@ int main(int argc, char *argv[])
     bool nfs = true;
     {
         List<fileName> roots(1, args.rootPath());
-        combineReduce(roots, ListOps::uniqueEqOp<fileName>());
+        Pstream::combineAllGather(roots, ListOps::uniqueEqOp<fileName>());
         nfs = (roots.size() == 1);
     }
 
@@ -3280,8 +3280,7 @@ int main(int argc, char *argv[])
                     )
                 );
                 meshDir[Pstream::myProcNo()] = fName.path();
-                Pstream::gatherList(meshDir);
-                Pstream::scatterList(meshDir);
+                Pstream::allGatherList(meshDir);
                 //Info<< "Per processor faces dirs:" << nl
                 //    << "    " << meshDir << nl << endl;
             }
diff --git a/applications/utilities/postProcessing/dataConversion/foamToEnsight/findCloudFields.H b/applications/utilities/postProcessing/dataConversion/foamToEnsight/findCloudFields.H
index c50e6e2a74..debff3288a 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToEnsight/findCloudFields.H
+++ b/applications/utilities/postProcessing/dataConversion/foamToEnsight/findCloudFields.H
@@ -104,12 +104,11 @@ if (timeDirs.size() && doLagrangian)
     {
         for (auto& cloudFields : regionCloudFields)
         {
-            Pstream::mapCombineGather
+            Pstream::mapCombineAllGather
             (
                 cloudFields,
                 HashTableOps::plusEqOp<word>()
             );
-            Pstream::broadcast(cloudFields);
         }
     }
 }
diff --git a/applications/utilities/postProcessing/dataConversion/foamToEnsight/foamToEnsight.C b/applications/utilities/postProcessing/dataConversion/foamToEnsight/foamToEnsight.C
index 45b816fe24..12b4891c3d 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToEnsight/foamToEnsight.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToEnsight/foamToEnsight.C
@@ -122,7 +122,7 @@ Usage
 #include "IOobjectList.H"
 #include "IOmanip.H"
 #include "OFstream.H"
-#include "PstreamCombineReduceOps.H"
+#include "Pstream.H"
 #include "HashOps.H"
 #include "regionProperties.H"
 
diff --git a/applications/utilities/postProcessing/lagrangian/particleTracks/particleTracks.C b/applications/utilities/postProcessing/lagrangian/particleTracks/particleTracks.C
index a21db2813a..a86ab41748 100644
--- a/applications/utilities/postProcessing/lagrangian/particleTracks/particleTracks.C
+++ b/applications/utilities/postProcessing/lagrangian/particleTracks/particleTracks.C
@@ -187,8 +187,7 @@ int main(int argc, char *argv[])
         const label maxNProcs = returnReduce(maxIds.size(), maxOp<label>());
         maxIds.resize(maxNProcs, -1);
 
-        Pstream::listCombineGather(maxIds, maxEqOp<label>());
-        Pstream::broadcast(maxIds);
+        Pstream::listCombineAllGather(maxIds, maxEqOp<label>());
 
         // From ids to count
         const labelList numIds = maxIds + 1;
diff --git a/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndexTemplates.C b/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndexTemplates.C
index 046df992ae..61c70dd170 100644
--- a/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndexTemplates.C
+++ b/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndexTemplates.C
@@ -40,8 +40,7 @@ Foam::Field<T> Foam::channelIndex::regionSum(const Field<T>& cellField) const
     }
 
     // Global sum
-    Pstream::listCombineGather(regionField, plusEqOp<T>());
-    Pstream::broadcast(regionField);
+    Pstream::listCombineAllGather(regionField, plusEqOp<T>());
 
     return regionField;
 }
diff --git a/applications/utilities/preProcessing/setFields/setFields.C b/applications/utilities/preProcessing/setFields/setFields.C
index e9487faacf..ea4ab60d77 100644
--- a/applications/utilities/preProcessing/setFields/setFields.C
+++ b/applications/utilities/preProcessing/setFields/setFields.C
@@ -287,8 +287,7 @@ bool setFaceFieldType
             }
         }
 
-        Pstream::listCombineGather(nChanged, plusEqOp<label>());
-        Pstream::broadcast(nChanged);
+        Pstream::listCombineAllGather(nChanged, plusEqOp<label>());
 
         auto& fieldBf = field.boundaryFieldRef();
 
diff --git a/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C b/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
index 231f58424d..14d4534b1b 100644
--- a/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
+++ b/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2016-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -456,12 +456,9 @@ int main(int argc, char *argv[])
     remoteCoarseSf[Pstream::myProcNo()] = localCoarseSf;
     remoteCoarseAgg[Pstream::myProcNo()] = localAgg;
 
-    Pstream::gatherList(remoteCoarseCf);
-    Pstream::scatterList(remoteCoarseCf);
-    Pstream::gatherList(remoteCoarseSf);
-    Pstream::scatterList(remoteCoarseSf);
-    Pstream::gatherList(remoteCoarseAgg);
-    Pstream::scatterList(remoteCoarseAgg);
+    Pstream::allGatherList(remoteCoarseCf);
+    Pstream::allGatherList(remoteCoarseSf);
+    Pstream::allGatherList(remoteCoarseAgg);
 
 
     globalIndex globalNumbering(nCoarseFaces);
diff --git a/applications/utilities/surface/surfaceRedistributePar/surfaceRedistributePar.C b/applications/utilities/surface/surfaceRedistributePar/surfaceRedistributePar.C
index 8728809ffd..fbf65a4ffd 100644
--- a/applications/utilities/surface/surfaceRedistributePar/surfaceRedistributePar.C
+++ b/applications/utilities/surface/surfaceRedistributePar/surfaceRedistributePar.C
@@ -198,8 +198,7 @@ int main(int argc, char *argv[])
                 boundBox(mesh.points(), false)
             ).extend(rndGen, 1e-3)
         );
-        Pstream::gatherList(meshBb);
-        Pstream::scatterList(meshBb);
+        Pstream::allGatherList(meshBb);
     }
 
     IOobject io
diff --git a/src/OpenFOAM/db/IOstreams/Pstreams/Pstream.H b/src/OpenFOAM/db/IOstreams/Pstreams/Pstream.H
index 268b9221c6..75cd715a3b 100644
--- a/src/OpenFOAM/db/IOstreams/Pstreams/Pstream.H
+++ b/src/OpenFOAM/db/IOstreams/Pstreams/Pstream.H
@@ -129,12 +129,13 @@ public:
 
     // Static Functions
 
-        // Broadcast
+    // Broadcast
 
         //- Broadcast buffer or string content to all processes in communicator
         using UPstream::broadcast;
 
-        //- Generic broadcast using streams to serialize/de-serialize
+        //- Generic broadcast using streams to serialize/de-serialize.
+        //  Not normally used directly.
         template<class T>
         static void genericBroadcast
         (
@@ -178,125 +179,232 @@ public:
         );
 
 
-        // Gather
+    // Gather
 
-            //- Gather data.
-            //- Apply \c bop to combine \c value from different processors
-            template<class T, class BinaryOp>
-            static void gather
-            (
-                const List<commsStruct>& comms,
-                T& value,
-                const BinaryOp& bop,
-                const int tag,
-                const label comm
-            );
+        //- Gather (reduce) data, appyling \c bop to combine \c value
+        //- from different processors. The basis for Foam::reduce().
+        //  Uses the specified communication schedule.
+        template<class T, class BinaryOp>
+        static void gather
+        (
+            const List<commsStruct>& comms,
+            T& value,
+            const BinaryOp& bop,
+            const int tag,
+            const label comm
+        );
 
-            //- Like above but switches between linear/tree communication
-            template<class T, class BinaryOp>
-            static void gather
-            (
-                T& value,
-                const BinaryOp& bop,
-                const int tag = UPstream::msgType(),
-                const label comm = UPstream::worldComm
-            );
+        //- Gather (reduce) data, applying \c bop to combine \c value
+        //- from different processors. The basis for Foam::reduce().
+        //  Uses linear/tree communication.
+        template<class T, class BinaryOp>
+        static void gather
+        (
+            T& value,
+            const BinaryOp& bop,
+            const int tag = UPstream::msgType(),
+            const label comm = UPstream::worldComm
+        );
 
 
-        // Gather/combine data
-        // Inplace combine values from processors.
-        // (Uses construct from Istream instead of <<)
+    // Gather/combine data
+    // Inplace combine values from processors.
+    // (Uses construct from Istream instead of <<)
 
-            template<class T, class CombineOp>
-            static void combineGather
-            (
-                const List<commsStruct>& comms,
-                T& value,
-                const CombineOp& cop,
-                const int tag,
-                const label comm
-            );
+        //- Gather data, applying \c cop to inplace combine \c value
+        //- from different processors.
+        //  Uses the specified communication schedule.
+        template<class T, class CombineOp>
+        static void combineGather
+        (
+            const List<commsStruct>& comms,
+            T& value,
+            const CombineOp& cop,
+            const int tag,
+            const label comm
+        );
 
-            //- Like above but switches between linear/tree communication
-            template<class T, class CombineOp>
-            static void combineGather
-            (
-                T& value,
-                const CombineOp& cop,
-                const int tag = UPstream::msgType(),
-                const label comm = UPstream::worldComm
-            );
+        //- Gather data, applying \c cop to inplace combine \c value
+        //- from different processors.
+        //  Uses linear/tree communication.
+        template<class T, class CombineOp>
+        static void combineGather
+        (
+            T& value,
+            const CombineOp& cop,
+            const int tag = UPstream::msgType(),
+            const label comm = UPstream::worldComm
+        );
+
+        //- Gather data, applying \c cop to inplace combine \c value
+        //- from different processors.
+        //- After completion all processors have the same data.
+        //  Uses the specified communication schedule.
+        //  Wraps combineGather/broadcast (may change in the future).
+        template<class T, class CombineOp>
+        static void combineAllGather
+        (
+            const List<commsStruct>& comms,
+            T& value,
+            const CombineOp& cop,
+            const int tag = UPstream::msgType(),
+            const label comm = UPstream::worldComm
+        );
+
+        //- Gather data, applying \c cop to inplace combine \c value
+        //- from different processors.
+        //- After completion all processors have the same data.
+        //  Uses linear/tree communication.
+        //  Wraps combineGather/broadcast (may change in the future).
+        template<class T, class CombineOp>
+        static void combineAllGather
+        (
+            T& value,
+            const CombineOp& cop,
+            const int tag = UPstream::msgType(),
+            const label comm = UPstream::worldComm
+        );
 
 
-        // Combine variants working on whole List at a time.
+    // Combine variants working on whole List at a time.
 
-            template<class T, class CombineOp>
-            static void listCombineGather
-            (
-                const List<commsStruct>& comms,
-                List<T>& value,
-                const CombineOp& cop,
-                const int tag,
-                const label comm
-            );
+        template<class T, class CombineOp>
+        static void listCombineGather
+        (
+            const List<commsStruct>& comms,
+            List<T>& values,
+            const CombineOp& cop,
+            const int tag,
+            const label comm
+        );
 
-            //- Like above but switches between linear/tree communication
-            template<class T, class CombineOp>
-            static void listCombineGather
-            (
-                List<T>& values,
-                const CombineOp& cop,
-                const int tag = UPstream::msgType(),
-                const label comm = UPstream::worldComm
-            );
+        //- Like above but switches between linear/tree communication
+        template<class T, class CombineOp>
+        static void listCombineGather
+        (
+            List<T>& values,
+            const CombineOp& cop,
+            const int tag = UPstream::msgType(),
+            const label comm = UPstream::worldComm
+        );
+
+        //- After completion all processors have the same data.
+        template<class T, class CombineOp>
+        static void listCombineAllGather
+        (
+            const List<commsStruct>& comms,
+            List<T>& values,
+            const CombineOp& cop,
+            const int tag,
+            const label comm
+        );
+
+        //- After completion all processors have the same data.
+        template<class T, class CombineOp>
+        static void listCombineAllGather
+        (
+            List<T>& values,
+            const CombineOp& cop,
+            const int tag = UPstream::msgType(),
+            const label comm = UPstream::worldComm
+        );
 
 
-        // Combine variants working on whole map at a time.
-        // Container needs iterators, find() and insert methods defined.
+    // Combine variants working on whole map at a time.
+    // Container needs iterators, find() and insert methods defined.
 
-            template<class Container, class CombineOp>
-            static void mapCombineGather
-            (
-                const List<commsStruct>& comms,
-                Container& values,
-                const CombineOp& cop,
-                const int tag,
-                const label comm
-            );
+        template<class Container, class CombineOp>
+        static void mapCombineGather
+        (
+            const List<commsStruct>& comms,
+            Container& values,
+            const CombineOp& cop,
+            const int tag,
+            const label comm
+        );
 
-            //- Like above but switches between linear/tree communication
-            template<class Container, class CombineOp>
-            static void mapCombineGather
-            (
-                Container& values,
-                const CombineOp& cop,
-                const int tag = UPstream::msgType(),
-                const label comm = UPstream::worldComm
-            );
+        //- Like above but switches between linear/tree communication
+        template<class Container, class CombineOp>
+        static void mapCombineGather
+        (
+            Container& values,
+            const CombineOp& cop,
+            const int tag = UPstream::msgType(),
+            const label comm = UPstream::worldComm
+        );
 
 
-        // Gather/scatter keeping the individual processor data separate.
-        // The values is a List of size UPstream::nProcs() where
-        // values[UPstream::myProcNo()] is the data for the current processor.
+        //- After completion all processors have the same data.
+        template<class Container, class CombineOp>
+        static void mapCombineAllGather
+        (
+            const List<commsStruct>& comms,
+            Container& values,
+            const CombineOp& cop,
+            const int tag,
+            const label comm
+        );
 
-            //- Gather data but keep individual values separate
-            template<class T>
-            static void gatherList
-            (
-                const List<commsStruct>& comms,
-                List<T>& values,
-                const int tag,
-                const label comm
-            );
+        //- After completion all processors have the same data.
+        template<class Container, class CombineOp>
+        static void mapCombineAllGather
+        (
+            Container& values,
+            const CombineOp& cop,
+            const int tag = UPstream::msgType(),
+            const label comm = UPstream::worldComm
+        );
 
-            //- Like above but switches between linear/tree communication
-            template<class T>
-            static void gatherList
-            (
-                List<T>& values,
-                const int tag = UPstream::msgType(),
-                const label comm = UPstream::worldComm
-            );
+
+    // Gather/scatter keeping the individual processor data separate.
+    // The values is a List of size UPstream::nProcs() where
+    // values[UPstream::myProcNo()] is the data for the current processor.
+
+        //- Gather data, but keep individual values separate.
+        //- Uses the specified communication schedule.
+        template<class T>
+        static void gatherList
+        (
+            const List<commsStruct>& comms,
+            List<T>& values,
+            const int tag,
+            const label comm
+        );
+
+        //- Gather data, but keep individual values separate.
+        //- Uses linear/tree communication.
+        template<class T>
+        static void gatherList
+        (
+            List<T>& values,
+            const int tag = UPstream::msgType(),
+            const label comm = UPstream::worldComm
+        );
+
+        //- Gather data, but keep individual values separate.
+        //- Uses the specified communication schedule.
+        //  After completion all processors have the same data.
+        //  Wraps gatherList/scatterList (may change in the future).
+        template<class T>
+        static void allGatherList
+        (
+            const List<commsStruct>& comms,
+            List<T>& values,
+            const int tag,
+            const label comm
+        );
+
+        //- Gather data, but keep individual values separate.
+        //- Uses linear/tree communication.
+        //  After completion all processors have the same data.
+        //  Wraps gatherList/scatterList (may change in the future).
+        template<class T>
+        static void allGatherList
+        (
+            List<T>& values,
+            const int tag = UPstream::msgType(),
+            const label comm = UPstream::worldComm
+        );
 
 
         // Scatter
diff --git a/src/OpenFOAM/db/IOstreams/Pstreams/PstreamCombineGather.C b/src/OpenFOAM/db/IOstreams/Pstreams/PstreamCombineGather.C
index 8429e3e27a..946838de21 100644
--- a/src/OpenFOAM/db/IOstreams/Pstreams/PstreamCombineGather.C
+++ b/src/OpenFOAM/db/IOstreams/Pstreams/PstreamCombineGather.C
@@ -143,26 +143,6 @@ void Foam::Pstream::combineGather
 }
 
 
-template<class T, class CombineOp>
-void Foam::Pstream::combineGather
-(
-    T& value,
-    const CombineOp& cop,
-    const int tag,
-    const label comm
-)
-{
-    combineGather
-    (
-        UPstream::whichCommunication(comm),
-        value,
-        cop,
-        tag,
-        comm
-    );
-}
-
-
 template<class T>
 void Foam::Pstream::combineScatter
 (
@@ -244,6 +224,26 @@ void Foam::Pstream::combineScatter
 }
 
 
+template<class T, class CombineOp>
+void Foam::Pstream::combineGather
+(
+    T& value,
+    const CombineOp& cop,
+    const int tag,
+    const label comm
+)
+{
+    Pstream::combineGather
+    (
+        UPstream::whichCommunication(comm),
+        value,
+        cop,
+        tag,
+        comm
+    );
+}
+
+
 template<class T>
 void Foam::Pstream::combineScatter
 (
@@ -260,6 +260,42 @@ void Foam::Pstream::combineScatter
 }
 
 
+template<class T, class CombineOp>
+void Foam::Pstream::combineAllGather
+(
+    const List<UPstream::commsStruct>& comms,
+    T& value,
+    const CombineOp& cop,
+    const int tag,
+    const label comm
+)
+{
+    Pstream::combineGather(comms, value, cop, tag, comm);
+    Pstream::broadcast(value, comm);
+}
+
+
+template<class T, class CombineOp>
+void Foam::Pstream::combineAllGather
+(
+    T& value,
+    const CombineOp& cop,
+    const int tag,
+    const label comm
+)
+{
+    if (UPstream::parRun() && UPstream::nProcs(comm) > 1)
+    {
+        const auto& comms = UPstream::whichCommunication(comm);
+
+        Pstream::combineGather(comms, value, cop, tag, comm);
+        Pstream::broadcast(value, comm);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 template<class T, class CombineOp>
 void Foam::Pstream::listCombineGather
 (
@@ -366,26 +402,6 @@ void Foam::Pstream::listCombineGather
 }
 
 
-template<class T, class CombineOp>
-void Foam::Pstream::listCombineGather
-(
-    List<T>& values,
-    const CombineOp& cop,
-    const int tag,
-    const label comm
-)
-{
-    listCombineGather
-    (
-        UPstream::whichCommunication(comm),
-        values,
-        cop,
-        tag,
-        comm
-    );
-}
-
-
 template<class T>
 void Foam::Pstream::listCombineScatter
 (
@@ -467,6 +483,26 @@ void Foam::Pstream::listCombineScatter
 }
 
 
+template<class T, class CombineOp>
+void Foam::Pstream::listCombineGather
+(
+    List<T>& values,
+    const CombineOp& cop,
+    const int tag,
+    const label comm
+)
+{
+    Pstream::listCombineGather
+    (
+        UPstream::whichCommunication(comm),
+        values,
+        cop,
+        tag,
+        comm
+    );
+}
+
+
 template<class T>
 void Foam::Pstream::listCombineScatter
 (
@@ -478,7 +514,7 @@ void Foam::Pstream::listCombineScatter
     #ifndef Foam_Pstream_scatter_nobroadcast
     Pstream::broadcast(values, comm);
     #else
-    listCombineScatter
+    Pstream::listCombineScatter
     (
         UPstream::whichCommunication(comm),
         values,
@@ -489,6 +525,42 @@ void Foam::Pstream::listCombineScatter
 }
 
 
+template<class T, class CombineOp>
+void Foam::Pstream::listCombineAllGather
+(
+    const List<UPstream::commsStruct>& comms,
+    List<T>& values,
+    const CombineOp& cop,
+    const int tag,
+    const label comm
+)
+{
+    Pstream::listCombineGather(comms, values, cop, tag, comm);
+    Pstream::broadcast(values, comm);
+}
+
+
+template<class T, class CombineOp>
+void Foam::Pstream::listCombineAllGather
+(
+    List<T>& values,
+    const CombineOp& cop,
+    const int tag,
+    const label comm
+)
+{
+    if (UPstream::parRun() && UPstream::nProcs(comm) > 1)
+    {
+        const auto& comms = UPstream::whichCommunication(comm);
+
+        Pstream::listCombineGather(comms, values, cop, tag, comm);
+        Pstream::broadcast(values, comm);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 template<class Container, class CombineOp>
 void Foam::Pstream::mapCombineGather
 (
@@ -570,26 +642,6 @@ void Foam::Pstream::mapCombineGather
 }
 
 
-template<class Container, class CombineOp>
-void Foam::Pstream::mapCombineGather
-(
-    Container& values,
-    const CombineOp& cop,
-    const int tag,
-    const label comm
-)
-{
-    mapCombineGather
-    (
-        UPstream::whichCommunication(comm),
-        values,
-        cop,
-        tag,
-        comm
-    );
-}
-
-
 template<class Container>
 void Foam::Pstream::mapCombineScatter
 (
@@ -652,6 +704,26 @@ void Foam::Pstream::mapCombineScatter
 }
 
 
+template<class Container, class CombineOp>
+void Foam::Pstream::mapCombineGather
+(
+    Container& values,
+    const CombineOp& cop,
+    const int tag,
+    const label comm
+)
+{
+    Pstream::mapCombineGather
+    (
+        UPstream::whichCommunication(comm),
+        values,
+        cop,
+        tag,
+        comm
+    );
+}
+
+
 template<class Container>
 void Foam::Pstream::mapCombineScatter
 (
@@ -663,7 +735,7 @@ void Foam::Pstream::mapCombineScatter
     #ifndef Foam_Pstream_scatter_nobroadcast
     Pstream::broadcast(values, comm);
     #else
-    mapCombineScatter
+    Pstream::mapCombineScatter
     (
         UPstream::whichCommunication(comm),
         values,
@@ -674,4 +746,38 @@ void Foam::Pstream::mapCombineScatter
 }
 
 
+template<class Container, class CombineOp>
+void Foam::Pstream::mapCombineAllGather
+(
+    const List<UPstream::commsStruct>& comms,
+    Container& values,
+    const CombineOp& cop,
+    const int tag,
+    const label comm
+)
+{
+    Pstream::mapCombineGather(comms, values, cop, tag, comm);
+    Pstream::broadcast(values, comm);
+}
+
+
+template<class Container, class CombineOp>
+void Foam::Pstream::mapCombineAllGather
+(
+    Container& values,
+    const CombineOp& cop,
+    const int tag,
+    const label comm
+)
+{
+    if (UPstream::parRun() && UPstream::nProcs(comm) > 1)
+    {
+        const auto& comms = UPstream::whichCommunication(comm);
+
+        Pstream::mapCombineGather(comms, values, cop, tag, comm);
+        Pstream::broadcast(values, comm);
+    }
+}
+
+
 // ************************************************************************* //
diff --git a/src/OpenFOAM/db/IOstreams/Pstreams/PstreamCombineReduceOps.H b/src/OpenFOAM/db/IOstreams/Pstreams/PstreamCombineReduceOps.H
index 0c200c9420..209f28eeef 100644
--- a/src/OpenFOAM/db/IOstreams/Pstreams/PstreamCombineReduceOps.H
+++ b/src/OpenFOAM/db/IOstreams/Pstreams/PstreamCombineReduceOps.H
@@ -38,7 +38,6 @@ Description
 #ifndef Foam_PstreamCombineReduceOps_H
 #define Foam_PstreamCombineReduceOps_H
 
-#include "UPstream.H"
 #include "Pstream.H"
 #include "ops.H"
 
@@ -59,8 +58,7 @@ void combineReduce
     const label comm
 )
 {
-    Pstream::combineGather(comms, value, cop, tag, comm);
-    Pstream::broadcast(value, comm);
+    Pstream::combineAllGather(comms, value, cop, tag, comm);
 }
 
 
@@ -73,13 +71,7 @@ void combineReduce
     const label comm = UPstream::worldComm
 )
 {
-    if (UPstream::parRun() && UPstream::nProcs(comm) > 1)
-    {
-        const auto& comms = UPstream::whichCommunication(comm);
-
-        Pstream::combineGather(comms, value, cop, tag, comm);
-        Pstream::broadcast(value, comm);
-    }
+    Pstream::combineAllGather(value, cop, tag, comm);
 }
 
 
diff --git a/src/OpenFOAM/db/IOstreams/Pstreams/PstreamGather.C b/src/OpenFOAM/db/IOstreams/Pstreams/PstreamGather.C
index 8bb4208b6d..1ea68fd28b 100644
--- a/src/OpenFOAM/db/IOstreams/Pstreams/PstreamGather.C
+++ b/src/OpenFOAM/db/IOstreams/Pstreams/PstreamGather.C
@@ -118,19 +118,6 @@ void Foam::Pstream::gather
 }
 
 
-template<class T, class BinaryOp>
-void Foam::Pstream::gather
-(
-    T& value,
-    const BinaryOp& bop,
-    const int tag,
-    const label comm
-)
-{
-    gather(UPstream::whichCommunication(comm), value, bop, tag, comm);
-}
-
-
 template<class T>
 void Foam::Pstream::scatter
 (
@@ -214,6 +201,19 @@ void Foam::Pstream::scatter
 }
 
 
+template<class T, class BinaryOp>
+void Foam::Pstream::gather
+(
+    T& value,
+    const BinaryOp& bop,
+    const int tag,
+    const label comm
+)
+{
+    Pstream::gather(UPstream::whichCommunication(comm), value, bop, tag, comm);
+}
+
+
 template<class T>
 void Foam::Pstream::scatter(T& value, const int tag, const label comm)
 {
diff --git a/src/OpenFOAM/db/IOstreams/Pstreams/PstreamGatherList.C b/src/OpenFOAM/db/IOstreams/Pstreams/PstreamGatherList.C
index db7bee31c9..cbecc6e6a2 100644
--- a/src/OpenFOAM/db/IOstreams/Pstreams/PstreamGatherList.C
+++ b/src/OpenFOAM/db/IOstreams/Pstreams/PstreamGatherList.C
@@ -187,13 +187,6 @@ void Foam::Pstream::gatherList
 }
 
 
-template<class T>
-void Foam::Pstream::gatherList(List<T>& values, const int tag, const label comm)
-{
-    gatherList(UPstream::whichCommunication(comm), values, tag, comm);
-}
-
-
 template<class T>
 void Foam::Pstream::scatterList
 (
@@ -324,6 +317,18 @@ void Foam::Pstream::scatterList
 }
 
 
+template<class T>
+void Foam::Pstream::gatherList
+(
+    List<T>& values,
+    const int tag,
+    const label comm
+)
+{
+    Pstream::gatherList(UPstream::whichCommunication(comm), values, tag, comm);
+}
+
+
 template<class T>
 void Foam::Pstream::scatterList
 (
@@ -332,7 +337,39 @@ void Foam::Pstream::scatterList
     const label comm
 )
 {
-    scatterList(UPstream::whichCommunication(comm), values, tag, comm);
+    Pstream::scatterList(UPstream::whichCommunication(comm), values, tag, comm);
+}
+
+
+template<class T>
+void Foam::Pstream::allGatherList
+(
+    const List<UPstream::commsStruct>& comms,
+    List<T>& values,
+    const int tag,
+    const label comm
+)
+{
+    Pstream::gatherList(comms, values, tag, comm);
+    Pstream::scatterList(comms, values, tag, comm);
+}
+
+
+template<class T>
+void Foam::Pstream::allGatherList
+(
+    List<T>& values,
+    const int tag,
+    const label comm
+)
+{
+    if (UPstream::parRun() && UPstream::nProcs(comm) > 1)
+    {
+        const auto& comms = UPstream::whichCommunication(comm);
+
+        Pstream::gatherList(comms, values, tag, comm);
+        Pstream::scatterList(comms, values, tag, comm);
+    }
 }
 
 
diff --git a/src/OpenFOAM/db/IOstreams/Pstreams/UPstream.H b/src/OpenFOAM/db/IOstreams/Pstreams/UPstream.H
index e82ca1b3c7..177365747f 100644
--- a/src/OpenFOAM/db/IOstreams/Pstreams/UPstream.H
+++ b/src/OpenFOAM/db/IOstreams/Pstreams/UPstream.H
@@ -82,17 +82,17 @@ public:
         {
             // Private Data
 
-                //- procID of above processor
+                //- The procID of the processor directly above
                 label above_;
 
-                //- procIDs of processors directly below me
+                //- The procIDs of all processors directly below
                 labelList below_;
 
                 //- procIDs of all processors below (so not just directly below)
                 labelList allBelow_;
 
                 //- procIDs of all processors not below.
-                //  (inverse set of allBelow_ and minus myProcNo)
+                //  Inverse set of allBelow_ without myProcNo.
                 labelList allNotBelow_;
 
 
@@ -125,21 +125,27 @@ public:
 
             // Member Functions
 
+                //- The procID of the processor directly above
                 label above() const noexcept
                 {
                     return above_;
                 }
 
+                //- The procIDs of all processors directly below
                 const labelList& below() const noexcept
                 {
                     return below_;
                 }
 
+                //- The procIDs of all processors below
+                //- (so not just directly below)
                 const labelList& allBelow() const noexcept
                 {
                     return allBelow_;
                 }
 
+                //- The procIDs of all processors that are above.
+                //- The inverse set of allBelow without myProcNo.
                 const labelList& allNotBelow() const noexcept
                 {
                     return allNotBelow_;
diff --git a/src/OpenFOAM/fields/Fields/Field/FieldOps.C b/src/OpenFOAM/fields/Fields/Field/FieldOps.C
index 061dff27a6..4d010e9e0b 100644
--- a/src/OpenFOAM/fields/Fields/Field/FieldOps.C
+++ b/src/OpenFOAM/fields/Fields/Field/FieldOps.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -25,7 +25,8 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "PstreamCombineReduceOps.H"
+#include "Pstream.H"
+#include "ops.H"
 #include <algorithm>
 
 // * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
@@ -171,7 +172,7 @@ Foam::Tuple2<T1,T2> Foam::FieldOps::findMinData
         result.second() = data[i];
     }
 
-    Foam::combineReduce(result, minFirstEqOp<T1>());
+    Pstream::combineAllGather(result, minFirstEqOp<T1>());
     return result;
 }
 
@@ -192,7 +193,7 @@ Foam::Tuple2<T1,T2> Foam::FieldOps::findMaxData
         result.second() = data[i];
     }
 
-    Foam::combineReduce(result, maxFirstEqOp<T1>());
+    Pstream::combineAllGather(result, maxFirstEqOp<T1>());
     return result;
 }
 
diff --git a/src/OpenFOAM/global/fileOperations/collatedFileOperation/hostCollatedFileOperation.C b/src/OpenFOAM/global/fileOperations/collatedFileOperation/hostCollatedFileOperation.C
index 55a02996dc..aab4e831bc 100644
--- a/src/OpenFOAM/global/fileOperations/collatedFileOperation/hostCollatedFileOperation.C
+++ b/src/OpenFOAM/global/fileOperations/collatedFileOperation/hostCollatedFileOperation.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2017-2018 OpenFOAM Foundation
-    Copyright (C) 2021 OpenCFD Ltd.
+    Copyright (C) 2021-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -106,8 +106,7 @@ Foam::labelList Foam::fileOperations::hostCollatedFileOperation::subRanks
 
         stringList hosts(Pstream::nProcs());
         hosts[Pstream::myProcNo()] = myHostName;
-        Pstream::gatherList(hosts);
-        Pstream::scatterList(hosts);
+        Pstream::allGatherList(hosts);
 
         // Collect procs with same hostname
         forAll(hosts, proci)
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGProcAgglomerations/procFacesGAMGProcAgglomeration/procFacesGAMGProcAgglomeration.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGProcAgglomerations/procFacesGAMGProcAgglomeration/procFacesGAMGProcAgglomeration.C
index 1678a9a5ed..bb2fc9cd5f 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGProcAgglomerations/procFacesGAMGProcAgglomeration/procFacesGAMGProcAgglomeration.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGProcAgglomerations/procFacesGAMGProcAgglomeration/procFacesGAMGProcAgglomeration.C
@@ -81,8 +81,7 @@ Foam::procFacesGAMGProcAgglomeration::singleCellMesh
         }
     }
 
-    Pstream::gatherList(procFaces, Pstream::msgType(), mesh.comm());
-    Pstream::scatterList(procFaces, Pstream::msgType(), mesh.comm());
+    Pstream::allGatherList(procFaces, Pstream::msgType(), mesh.comm());
 
     autoPtr<lduPrimitiveMesh> singleCellMeshPtr;
 
diff --git a/src/OpenFOAM/meshes/ProcessorTopology/ProcessorTopology.C b/src/OpenFOAM/meshes/ProcessorTopology/ProcessorTopology.C
index 1343bd05d4..359f7ae0b9 100644
--- a/src/OpenFOAM/meshes/ProcessorTopology/ProcessorTopology.C
+++ b/src/OpenFOAM/meshes/ProcessorTopology/ProcessorTopology.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -113,8 +114,7 @@ Foam::ProcessorTopology<Container, ProcPatch>::ProcessorTopology
             procNeighbours(this->size(), patches);
 
         // Distribute to all processors
-        Pstream::gatherList(*this, Pstream::msgType(), comm);
-        Pstream::scatterList(*this, Pstream::msgType(), comm);
+        Pstream::allGatherList(*this, Pstream::msgType(), comm);
     }
 
     if
diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C
index c4ea5b3e3d..0741da6531 100644
--- a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C
+++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C
@@ -28,7 +28,6 @@ License
 
 #include "globalMeshData.H"
 #include "Pstream.H"
-#include "PstreamCombineReduceOps.H"
 #include "processorPolyPatch.H"
 #include "globalPoints.H"
 #include "polyMesh.H"
@@ -1926,8 +1925,8 @@ Foam::pointField Foam::globalMeshData::geometricSharedPoints() const
     // Get coords of my shared points
     pointField sharedPoints(mesh_.points(), sharedPointLabels());
 
-    // Append from all processors
-    combineReduce(sharedPoints, ListOps::appendEqOp<point>());
+    // Append from all processors, globally consistent
+    Pstream::combineAllGather(sharedPoints, ListOps::appendEqOp<point>());
 
     // Merge tolerance
     scalar tolDim = matchTol_ * mesh_.bounds().mag();
diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeBaseTemplates.C b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeBaseTemplates.C
index 03e6fcf2ea..6d191752d0 100644
--- a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeBaseTemplates.C
+++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeBaseTemplates.C
@@ -28,7 +28,6 @@ License
 
 #include "Pstream.H"
 #include "PstreamBuffers.H"
-#include "PstreamCombineReduceOps.H"
 #include "flipOp.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C
index c22f8cd01b..5e643d1012 100644
--- a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C
+++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C
@@ -27,7 +27,6 @@ License
 
 #include "Pstream.H"
 #include "PstreamBuffers.H"
-#include "PstreamCombineReduceOps.H"
 #include "globalIndexAndTransform.H"
 #include "transformField.H"
 #include "flipOp.H"
diff --git a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C
index 6c52daaa22..d71112b0da 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C
@@ -1018,13 +1018,11 @@ bool Foam::polyBoundaryMesh::checkParallelSync(const bool report) const
 
     List<wordList> allNames(Pstream::nProcs());
     allNames[Pstream::myProcNo()] = names;
-    Pstream::gatherList(allNames);
-    Pstream::scatterList(allNames);
+    Pstream::allGatherList(allNames);
 
     List<wordList> allTypes(Pstream::nProcs());
     allTypes[Pstream::myProcNo()] = types;
-    Pstream::gatherList(allTypes);
-    Pstream::scatterList(allTypes);
+    Pstream::allGatherList(allTypes);
 
     // Have every processor check but print error on master
     // (in processor sequence).
diff --git a/src/OpenFOAM/meshes/polyMesh/syncTools/syncTools.H b/src/OpenFOAM/meshes/polyMesh/syncTools/syncTools.H
index 08323f6f3a..be8f76ef98 100644
--- a/src/OpenFOAM/meshes/polyMesh/syncTools/syncTools.H
+++ b/src/OpenFOAM/meshes/polyMesh/syncTools/syncTools.H
@@ -32,7 +32,7 @@ Description
 
     Require
     - combineOperator (e.g. sumEqOp - not sumOp!) that is defined for the
-      type and combineReduce(UList\<T\>, combineOperator) should be defined.
+      type be defined.
     - null value which gets overridden by any valid value.
     - transform function
 
@@ -57,7 +57,7 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declarations
+// Forward Declarations
 class polyBoundaryMesh;
 
 /*---------------------------------------------------------------------------*\
diff --git a/src/OpenFOAM/meshes/polyMesh/zones/ZoneMesh/ZoneMesh.C b/src/OpenFOAM/meshes/polyMesh/zones/ZoneMesh/ZoneMesh.C
index 22103ed846..96c959d3b2 100644
--- a/src/OpenFOAM/meshes/polyMesh/zones/ZoneMesh/ZoneMesh.C
+++ b/src/OpenFOAM/meshes/polyMesh/zones/ZoneMesh/ZoneMesh.C
@@ -765,13 +765,11 @@ bool Foam::ZoneMesh<ZoneType, MeshType>::checkParallelSync
     // Collect all names
     List<wordList> allNames(Pstream::nProcs());
     allNames[Pstream::myProcNo()] = this->names();
-    Pstream::gatherList(allNames);
-    Pstream::scatterList(allNames);
+    Pstream::allGatherList(allNames);
 
     List<wordList> allTypes(Pstream::nProcs());
     allTypes[Pstream::myProcNo()] = this->types();
-    Pstream::gatherList(allTypes);
-    Pstream::scatterList(allTypes);
+    Pstream::allGatherList(allTypes);
 
     // Have every processor check but only master print error.
 
diff --git a/src/atmosphericModels/porosityModels/powerLawLopesdaCosta/powerLawLopesdaCosta.C b/src/atmosphericModels/porosityModels/powerLawLopesdaCosta/powerLawLopesdaCosta.C
index 15d68357d1..e5c903e2d5 100644
--- a/src/atmosphericModels/porosityModels/powerLawLopesdaCosta/powerLawLopesdaCosta.C
+++ b/src/atmosphericModels/porosityModels/powerLawLopesdaCosta/powerLawLopesdaCosta.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2018 OpenFOAM Foundation
-    Copyright (C) 2020,2022 OpenCFD Ltd.
+    Copyright (C) 2020-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -175,10 +175,8 @@ Foam::porosityModels::powerLawLopesdaCostaZone::powerLawLopesdaCostaZone
         groundSurfaceProcTris[Pstream::myProcNo()] = groundSurface;
         groundSurfaceProcPoints[Pstream::myProcNo()] = groundSurface.points();
 
-        Pstream::gatherList(groundSurfaceProcTris);
-        Pstream::scatterList(groundSurfaceProcTris);
-        Pstream::gatherList(groundSurfaceProcPoints);
-        Pstream::scatterList(groundSurfaceProcPoints);
+        Pstream::allGatherList(groundSurfaceProcTris);
+        Pstream::allGatherList(groundSurfaceProcPoints);
 
         label nTris = 0;
         forAll(groundSurfaceProcTris, i)
diff --git a/src/dynamicFaMesh/interfaceTrackingFvMesh/freeSurfacePointDisplacement.C b/src/dynamicFaMesh/interfaceTrackingFvMesh/freeSurfacePointDisplacement.C
index 66407a3371..46e1deffda 100644
--- a/src/dynamicFaMesh/interfaceTrackingFvMesh/freeSurfacePointDisplacement.C
+++ b/src/dynamicFaMesh/interfaceTrackingFvMesh/freeSurfacePointDisplacement.C
@@ -30,7 +30,6 @@ License
 #include "emptyFaPatch.H"
 #include "wedgeFaPatch.H"
 #include "wallFvPatch.H"
-#include "PstreamCombineReduceOps.H"
 #include "coordinateSystem.H"
 #include "unitConversion.H"
 #include "scalarMatrices.H"
@@ -466,8 +465,7 @@ Foam::interfaceTrackingFvMesh::pointDisplacement()
                 }
             }
 
-            Pstream::gatherList(procLsPoints);
-            Pstream::scatterList(procLsPoints);
+            Pstream::allGatherList(procLsPoints);
 
             if (curSharedPointIndex != -1)
             {
diff --git a/src/dynamicFaMesh/interfaceTrackingFvMesh/functionObjects/pointHistory/pointHistory.C b/src/dynamicFaMesh/interfaceTrackingFvMesh/functionObjects/pointHistory/pointHistory.C
index fe61625241..d1060b286f 100644
--- a/src/dynamicFaMesh/interfaceTrackingFvMesh/functionObjects/pointHistory/pointHistory.C
+++ b/src/dynamicFaMesh/interfaceTrackingFvMesh/functionObjects/pointHistory/pointHistory.C
@@ -119,8 +119,7 @@ Foam::pointHistory::pointHistory
         }
     }
 
-    Pstream::gatherList(minDist);
-    Pstream::scatterList(minDist);
+    Pstream::allGatherList(minDist);
 
     processor_ = -1;
     scalar min = GREAT;
diff --git a/src/dynamicMesh/extrudePatchMesh/extrudePatchMesh.C b/src/dynamicMesh/extrudePatchMesh/extrudePatchMesh.C
index 02e1c285ff..628acf544f 100644
--- a/src/dynamicMesh/extrudePatchMesh/extrudePatchMesh.C
+++ b/src/dynamicMesh/extrudePatchMesh/extrudePatchMesh.C
@@ -214,8 +214,7 @@ void Foam::extrudePatchMesh::extrudeMesh(const List<polyPatch*>& regionPatches)
             }
 
             // Reduce
-            Pstream::mapCombineGather(globalSum, plusEqOp<point>());
-            Pstream::broadcast(globalSum);
+            Pstream::mapCombineAllGather(globalSum, plusEqOp<point>());
 
             forAll(localToGlobalRegion, localRegionI)
             {
diff --git a/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C b/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C
index 044a935d90..69c6c08313 100644
--- a/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C
+++ b/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C
@@ -27,7 +27,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "fvMeshDistribute.H"
-#include "PstreamCombineReduceOps.H"
 #include "fvMeshAdder.H"
 #include "faceCoupleInfo.H"
 #include "processorFvPatchField.H"
diff --git a/src/dynamicMesh/fvMeshDistribute/fvMeshDistributeTemplates.C b/src/dynamicMesh/fvMeshDistribute/fvMeshDistributeTemplates.C
index d55137a27a..38958bda5d 100644
--- a/src/dynamicMesh/fvMeshDistribute/fvMeshDistributeTemplates.C
+++ b/src/dynamicMesh/fvMeshDistribute/fvMeshDistributeTemplates.C
@@ -399,8 +399,7 @@ void Foam::fvMeshDistribute::getFieldNames
     {
         List<wordList> allNames(Pstream::nProcs());
         allNames[Pstream::myProcNo()] = list;
-        Pstream::gatherList(allNames);
-        Pstream::scatterList(allNames);
+        Pstream::allGatherList(allNames);
 
         for (const int proci : Pstream::subProcs())
         {
diff --git a/src/dynamicMesh/fvMeshSubset/fvMeshSubset.C b/src/dynamicMesh/fvMeshSubset/fvMeshSubset.C
index 59419c2a20..a60e205521 100644
--- a/src/dynamicMesh/fvMeshSubset/fvMeshSubset.C
+++ b/src/dynamicMesh/fvMeshSubset/fvMeshSubset.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2015-2021 OpenCFD Ltd.
+    Copyright (C) 2015-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -1011,14 +1011,12 @@ void Foam::fvMeshSubset::setCellSubset
         List<wordList> patchNames(Pstream::nProcs());
         patchNames[Pstream::myProcNo()] = oldPatches.names();
         patchNames[Pstream::myProcNo()].setSize(nextPatchID);
-        Pstream::gatherList(patchNames);
-        Pstream::scatterList(patchNames);
+        Pstream::allGatherList(patchNames);
 
         // Get patch sizes (up to nextPatchID).
         // Note that up to nextPatchID the globalPatchMap is an identity so
         // no need to index through that.
-        Pstream::listCombineGather(globalPatchSizes, plusEqOp<label>());
-        Pstream::broadcast(globalPatchSizes);
+        Pstream::listCombineAllGather(globalPatchSizes, plusEqOp<label>());
 
         // Now all processors have all the patchnames.
         // Decide: if all processors have the same patch names and size is zero
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C
index bf65c55859..1073465ab1 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C
@@ -436,8 +436,7 @@ Foam::scalar Foam::hexRef8::getLevel0EdgeLength() const
 
     // Get the minimum per level over all processors. Note minimum so if
     // cells are not cubic we use the smallest edge side.
-    Pstream::listCombineGather(typEdgeLenSqr, minEqOp<scalar>());
-    Pstream::broadcast(typEdgeLenSqr);
+    Pstream::listCombineAllGather(typEdgeLenSqr, minEqOp<scalar>());
 
     if (debug)
     {
@@ -471,8 +470,7 @@ Foam::scalar Foam::hexRef8::getLevel0EdgeLength() const
         }
     }
 
-    Pstream::listCombineGather(maxEdgeLenSqr, maxEqOp<scalar>());
-    Pstream::broadcast(maxEdgeLenSqr);
+    Pstream::listCombineAllGather(maxEdgeLenSqr, maxEqOp<scalar>());
 
     if (debug)
     {
diff --git a/src/finiteArea/faMesh/faGlobalMeshData/faGlobalMeshData.C b/src/finiteArea/faMesh/faGlobalMeshData/faGlobalMeshData.C
index 8d341b65e1..2df0739f11 100644
--- a/src/finiteArea/faMesh/faGlobalMeshData/faGlobalMeshData.C
+++ b/src/finiteArea/faMesh/faGlobalMeshData/faGlobalMeshData.C
@@ -33,7 +33,6 @@ Author
 #include "faGlobalMeshData.H"
 #include "faMesh.H"
 #include "globalMeshData.H"
-#include "PstreamCombineReduceOps.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -110,7 +109,7 @@ void Foam::faGlobalMeshData::updateMesh()
 
     sharedPointLabels_ = sharedPointLabels.toc();
 
-    combineReduce(globalList, plusEqOp<labelField>());
+    Pstream::combineAllGather(globalList, plusEqOp<labelField>());
 
     nGlobalPoints_ = 0;
     for (label i=0; i<globalList.size(); ++i)
diff --git a/src/finiteArea/faMesh/faMeshDemandDrivenData.C b/src/finiteArea/faMesh/faMeshDemandDrivenData.C
index e56b1c13c2..6b523434ea 100644
--- a/src/finiteArea/faMesh/faMeshDemandDrivenData.C
+++ b/src/finiteArea/faMesh/faMeshDemandDrivenData.C
@@ -34,7 +34,6 @@ License
 #include "fac.H"
 #include "processorFaPatch.H"
 #include "wedgeFaPatch.H"
-#include "PstreamCombineReduceOps.H"
 #include "cartesianCS.H"
 #include "scalarMatrices.H"
 #include "processorFaPatchFields.H"
@@ -1146,7 +1145,7 @@ void Foam::faMesh::calcPointAreaNormals_orig(vectorField& result) const
             gpNormals[addr[i]] += spNormals[i];
         }
 
-        combineReduce(gpNormals, plusEqOp<vectorField>());
+        Pstream::combineAllGather(gpNormals, plusEqOp<vectorField>());
 
         // Extract local data
         forAll(addr, i)
@@ -1398,7 +1397,7 @@ void Foam::faMesh::calcPointAreaNormals(vectorField& result) const
             gpNormals[addr[i]] += spNormals[i];
         }
 
-        combineReduce(gpNormals, plusEqOp<vectorField>());
+        Pstream::combineAllGather(gpNormals, plusEqOp<vectorField>());
 
         // Extract local data
         forAll(addr, i)
@@ -1918,8 +1917,7 @@ void Foam::faMesh::calcPointAreaNormalsByQuadricsFit(vectorField& result) const
                 tol = 0.001*mag(bb.max() - bb.min());
             }
 
-            Pstream::gatherList(procLsPoints);
-            Pstream::scatterList(procLsPoints);
+            Pstream::allGatherList(procLsPoints);
 
             if (curSharedPointIndex != -1)
             {
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C
index bd80ea9766..ed9d48e091 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015 OpenFOAM Foundation
-    Copyright (C) 2016-2021 OpenCFD Ltd.
+    Copyright (C) 2016-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -161,15 +161,10 @@ void Foam::turbulentDFSEMInletFvPatchVectorField::initialisePatch()
         }
     }
 
-    for (auto& s : sumTriMagSf_)
-    {
-        s = 0.0;
-    }
-
+    sumTriMagSf_ = Zero;
     sumTriMagSf_[Pstream::myProcNo() + 1] = sum(triMagSf);
 
-    Pstream::listCombineGather(sumTriMagSf_, maxEqOp<scalar>());
-    Pstream::broadcast(sumTriMagSf_);
+    Pstream::listCombineAllGather(sumTriMagSf_, maxEqOp<scalar>());
 
     for (label i = 1; i < triMagSf.size(); ++i)
     {
@@ -494,8 +489,7 @@ void Foam::turbulentDFSEMInletFvPatchVectorField::calcOverlappingProcEddies
 
     List<boundBox> patchBBs(Pstream::nProcs());
     patchBBs[Pstream::myProcNo()] = patchBounds_;
-    Pstream::gatherList(patchBBs);
-    Pstream::scatterList(patchBBs);
+    Pstream::allGatherList(patchBBs);
 
     // Per processor indices into all segments to send
     List<DynamicList<label>> dynSendMap(Pstream::nProcs());
@@ -537,8 +531,7 @@ void Foam::turbulentDFSEMInletFvPatchVectorField::calcOverlappingProcEddies
     {
         sendSizes[Pstream::myProcNo()][procI] = sendMap[procI].size();
     }
-    Pstream::gatherList(sendSizes);
-    Pstream::scatterList(sendSizes);
+    Pstream::allGatherList(sendSizes);
 
     // Determine order of receiving
     labelListList constructMap(Pstream::nProcs());
diff --git a/src/functionObjects/field/Curle/Curle.C b/src/functionObjects/field/Curle/Curle.C
index eb5cef166b..9442cfa7d5 100644
--- a/src/functionObjects/field/Curle/Curle.C
+++ b/src/functionObjects/field/Curle/Curle.C
@@ -223,8 +223,7 @@ bool Foam::functionObjects::Curle::execute()
 
     pDash /= 4*mathematical::pi;
 
-    Pstream::listCombineGather(pDash, plusEqOp<scalar>());
-    Pstream::broadcast(pDash);
+    Pstream::listCombineAllGather(pDash, plusEqOp<scalar>());
 
     if (surfaceWriterPtr_)
     {
diff --git a/src/functionObjects/field/columnAverage/columnAverageTemplates.C b/src/functionObjects/field/columnAverage/columnAverageTemplates.C
index 33c94412b3..50d40abe8f 100644
--- a/src/functionObjects/field/columnAverage/columnAverageTemplates.C
+++ b/src/functionObjects/field/columnAverage/columnAverageTemplates.C
@@ -87,10 +87,8 @@ bool Foam::functionObjects::columnAverage::columnAverageField
         }
 
         // Global sum
-        Pstream::listCombineGather(regionField, plusEqOp<Type>());
-        Pstream::listCombineGather(regionCount, plusEqOp<label>());
-        Pstream::broadcast(regionField);
-        Pstream::broadcast(regionCount);
+        Pstream::listCombineAllGather(regionField, plusEqOp<Type>());
+        Pstream::listCombineAllGather(regionCount, plusEqOp<label>());
 
         forAll(regionField, regioni)
         {
diff --git a/src/functionObjects/field/extractEulerianParticles/extractEulerianParticles/extractEulerianParticles.C b/src/functionObjects/field/extractEulerianParticles/extractEulerianParticles/extractEulerianParticles.C
index 4fbcabdbaa..b97268833d 100644
--- a/src/functionObjects/field/extractEulerianParticles/extractEulerianParticles/extractEulerianParticles.C
+++ b/src/functionObjects/field/extractEulerianParticles/extractEulerianParticles/extractEulerianParticles.C
@@ -333,8 +333,7 @@ void Foam::functionObjects::extractEulerianParticles::calculateAddressing
 
     // Create map from new regions to slots in particles list
     // - filter through new-to-new addressing to identify new particles
-    Pstream::listCombineGather(newToNewRegion, maxEqOp<label>());
-    Pstream::broadcast(newToNewRegion);
+    Pstream::listCombineAllGather(newToNewRegion, maxEqOp<label>());
 
     label nParticle = -1;
     labelHashSet newRegions;
@@ -353,8 +352,8 @@ void Foam::functionObjects::extractEulerianParticles::calculateAddressing
 
     // Accumulate old region data or create a new particle if there is no
     // mapping from the old-to-new region
-    Pstream::listCombineGather(oldToNewRegion, maxEqOp<label>());
-    Pstream::broadcast(oldToNewRegion);
+    Pstream::listCombineAllGather(oldToNewRegion, maxEqOp<label>());
+
     List<eulerianParticle> newParticles(newRegionToParticleMap.size());
     forAll(oldToNewRegion, oldRegioni)
     {
diff --git a/src/functionObjects/field/fieldMinMax/fieldMinMaxTemplates.C b/src/functionObjects/field/fieldMinMax/fieldMinMaxTemplates.C
index cdc9d0a5a9..092969850f 100644
--- a/src/functionObjects/field/fieldMinMax/fieldMinMaxTemplates.C
+++ b/src/functionObjects/field/fieldMinMax/fieldMinMaxTemplates.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2015-2019 OpenCFD Ltd.
+    Copyright (C) 2015-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -186,19 +186,13 @@ void Foam::functionObjects::fieldMinMax::calcMinMaxFieldType
     }
 
     // Collect info from all processors and output
-    Pstream::gatherList(minVs);
-    Pstream::scatterList(minVs);
-    Pstream::gatherList(minCells);
-    Pstream::scatterList(minCells);
-    Pstream::gatherList(minCs);
-    Pstream::scatterList(minCs);
+    Pstream::allGatherList(minVs);
+    Pstream::allGatherList(minCells);
+    Pstream::allGatherList(minCs);
 
-    Pstream::gatherList(maxVs);
-    Pstream::scatterList(maxVs);
-    Pstream::gatherList(maxCells);
-    Pstream::scatterList(maxCells);
-    Pstream::gatherList(maxCs);
-    Pstream::scatterList(maxCs);
+    Pstream::allGatherList(maxVs);
+    Pstream::allGatherList(maxCells);
+    Pstream::allGatherList(maxCs);
 
     minId = findMin(minVs);
     const Type& minValue = minVs[minId];
diff --git a/src/functionObjects/field/fieldValues/fieldValue/fieldValueTemplates.C b/src/functionObjects/field/fieldValues/fieldValue/fieldValueTemplates.C
index 2f8bb0fc38..69d2614ae1 100644
--- a/src/functionObjects/field/fieldValues/fieldValue/fieldValueTemplates.C
+++ b/src/functionObjects/field/fieldValues/fieldValue/fieldValueTemplates.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2017 OpenCFD Ltd.
+    Copyright (C) 2015-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -40,8 +40,7 @@ void Foam::functionObjects::fieldValue::combineFields(Field<Type>& field)
         List<Field<Type>> allValues(Pstream::nProcs());
         allValues[Pstream::myProcNo()] = field;
 
-        Pstream::gatherList(allValues);
-        Pstream::scatterList(allValues);
+        Pstream::allGatherList(allValues);
 
         field =
             ListListOps::combine<Field<Type>>
diff --git a/src/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C b/src/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C
index 250939c015..96626fda2d 100644
--- a/src/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C
+++ b/src/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C
@@ -66,8 +66,7 @@ static Map<Type> regionSum(const regionSplit& regions, const Field<Type>& fld)
         regionToSum(regioni, Type(Zero)) += fld[celli];
     }
 
-    Pstream::mapCombineGather(regionToSum, plusEqOp<Type>());
-    Pstream::broadcast(regionToSum);
+    Pstream::mapCombineAllGather(regionToSum, plusEqOp<Type>());
 
     return regionToSum;
 }
@@ -215,8 +214,7 @@ Foam::functionObjects::regionSizeDistribution::findPatchRegions
 
 
     // Make sure all the processors have the same set of regions
-    Pstream::mapCombineGather(patchRegions, minEqOp<label>());
-    Pstream::broadcast(patchRegions);
+    Pstream::mapCombineAllGather(patchRegions, minEqOp<label>());
 
     return patchRegions;
 }
diff --git a/src/functionObjects/forces/forces/forces.C b/src/functionObjects/forces/forces/forces.C
index 8040ee589a..2f3648012f 100644
--- a/src/functionObjects/forces/forces/forces.C
+++ b/src/functionObjects/forces/forces/forces.C
@@ -1063,10 +1063,8 @@ void Foam::functionObjects::forces::calcForcesMoment()
         }
     }
 
-    Pstream::listCombineGather(force_, plusEqOp<vectorField>());
-    Pstream::listCombineGather(moment_, plusEqOp<vectorField>());
-    Pstream::broadcast(force_);
-    Pstream::broadcast(moment_);
+    Pstream::listCombineAllGather(force_, plusEqOp<vectorField>());
+    Pstream::listCombineAllGather(moment_, plusEqOp<vectorField>());
 }
 
 
diff --git a/src/functionObjects/forces/propellerInfo/propellerInfo.C b/src/functionObjects/forces/propellerInfo/propellerInfo.C
index 11da8c638e..e74a2472db 100644
--- a/src/functionObjects/forces/propellerInfo/propellerInfo.C
+++ b/src/functionObjects/forces/propellerInfo/propellerInfo.C
@@ -468,8 +468,7 @@ void Foam::functionObjects::propellerInfo::updateSampleDiskCells()
         }
     }
 
-    Pstream::listCombineGather(pointMask_, orEqOp<bool>());
-    Pstream::broadcast(pointMask_);
+    Pstream::listCombineAllGather(pointMask_, orEqOp<bool>());
 }
 
 
@@ -774,8 +773,7 @@ Foam::tmp<Foam::Field<Type>> Foam::functionObjects::propellerInfo::interpolate
         }
     }
 
-    Pstream::listCombineGather(field, maxEqOp<Type>());
-    Pstream::broadcast(field);
+    Pstream::listCombineAllGather(field, maxEqOp<Type>());
 
     return tfield;
 }
diff --git a/src/functionObjects/utilities/areaWrite/areaWrite.C b/src/functionObjects/utilities/areaWrite/areaWrite.C
index 6c4a7bdaf3..d0bc87ca40 100644
--- a/src/functionObjects/utilities/areaWrite/areaWrite.C
+++ b/src/functionObjects/utilities/areaWrite/areaWrite.C
@@ -266,11 +266,8 @@ bool Foam::areaWrite::write()
             selected = areaMesh.thisDb().classes(fieldSelection_);
         }
 
-        if (Pstream::parRun())
-        {
-            Pstream::mapCombineGather(selected, HashSetOps::plusEqOp<word>());
-            Pstream::broadcast(selected);
-        }
+        // Parallel consistency (no-op in serial)
+        Pstream::mapCombineAllGather(selected, HashSetOps::plusEqOp<word>());
 
         missed.clear();
 
diff --git a/src/lagrangian/basic/Cloud/Cloud.C b/src/lagrangian/basic/Cloud/Cloud.C
index 3bac8a8f8d..c5f062124b 100644
--- a/src/lagrangian/basic/Cloud/Cloud.C
+++ b/src/lagrangian/basic/Cloud/Cloud.C
@@ -29,7 +29,7 @@ License
 #include "Cloud.H"
 #include "processorPolyPatch.H"
 #include "globalMeshData.H"
-#include "PstreamCombineReduceOps.H"
+#include "PstreamBuffers.H"
 #include "mapPolyMesh.H"
 #include "Time.H"
 #include "OFstream.H"
diff --git a/src/lagrangian/basic/Cloud/CloudIO.C b/src/lagrangian/basic/Cloud/CloudIO.C
index 4e571c68e9..ba78da6857 100644
--- a/src/lagrangian/basic/Cloud/CloudIO.C
+++ b/src/lagrangian/basic/Cloud/CloudIO.C
@@ -104,8 +104,7 @@ void Foam::Cloud<ParticleType>::writeCloudUniformProperties() const
     labelList np(Pstream::nProcs(), Zero);
     np[Pstream::myProcNo()] = ParticleType::particleCount_;
 
-    Pstream::listCombineGather(np, maxEqOp<label>());
-    Pstream::broadcast(np);
+    Pstream::listCombineAllGather(np, maxEqOp<label>());
 
     uniformPropsDict.add
     (
diff --git a/src/lagrangian/basic/InteractionLists/InteractionLists.C b/src/lagrangian/basic/InteractionLists/InteractionLists.C
index 85653a21d0..eb63c6feca 100644
--- a/src/lagrangian/basic/InteractionLists/InteractionLists.C
+++ b/src/lagrangian/basic/InteractionLists/InteractionLists.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -58,9 +58,7 @@ void Foam::InteractionLists<ParticleType>::buildInteractionLists()
 
     allExtendedProcBbs[Pstream::myProcNo()] = extendedProcBb;
 
-    Pstream::gatherList(allExtendedProcBbs);
-
-    Pstream::scatterList(allExtendedProcBbs);
+    Pstream::allGatherList(allExtendedProcBbs);
 
     List<treeBoundBox> extendedProcBbsInRange;
     List<label> extendedProcBbsTransformIndex;
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/CellZoneInjection/CellZoneInjection.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/CellZoneInjection/CellZoneInjection.C
index 250f3e6691..cade4faea9 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/CellZoneInjection/CellZoneInjection.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/CellZoneInjection/CellZoneInjection.C
@@ -120,8 +120,7 @@ void Foam::CellZoneInjection<CloudType>::setPositions
         globalPositions.localStart(Pstream::myProcNo())
     ) = positions;
 
-    Pstream::listCombineGather(allPositions, minEqOp<point>());
-    Pstream::broadcast(allPositions);
+    Pstream::listCombineAllGather(allPositions, minEqOp<point>());
 
     // Gather local cell tet and tet-point Ids, but leave non-local ids set -1
     SubList<label>
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/patchInjectionBase.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/patchInjectionBase.C
index c1c216f033..7f8e012bf6 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/patchInjectionBase.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/patchInjectionBase.C
@@ -111,15 +111,10 @@ void Foam::patchInjectionBase::updateMesh(const polyMesh& mesh)
         }
     }
 
-    forAll(sumTriMagSf_, i)
-    {
-        sumTriMagSf_[i] = 0.0;
-    }
-
+    sumTriMagSf_ = Zero;
     sumTriMagSf_[Pstream::myProcNo() + 1] = sum(triMagSf);
 
-    Pstream::listCombineGather(sumTriMagSf_, maxEqOp<scalar>());
-    Pstream::broadcast(sumTriMagSf_);
+    Pstream::listCombineAllGather(sumTriMagSf_, maxEqOp<scalar>());
 
     for (label i = 1; i < triMagSf.size(); i++)
     {
diff --git a/src/lumpedPointMotion/movement/lumpedPointMovement.C b/src/lumpedPointMotion/movement/lumpedPointMovement.C
index 3c6f2b4759..5b2d8b6427 100644
--- a/src/lumpedPointMotion/movement/lumpedPointMovement.C
+++ b/src/lumpedPointMotion/movement/lumpedPointMovement.C
@@ -880,8 +880,7 @@ Foam::List<Foam::scalar> Foam::lumpedPointMovement::areas
         }
     }
 
-    Pstream::listCombineGather(zoneAreas, plusEqOp<scalar>());
-    Pstream::broadcast(zoneAreas);
+    Pstream::listCombineAllGather(zoneAreas, plusEqOp<scalar>());
 
     return zoneAreas;
 }
@@ -1009,10 +1008,8 @@ bool Foam::lumpedPointMovement::forcesAndMoments
         Info<<"No pressure field" << endl;
     }
 
-    Pstream::listCombineGather(forces, plusEqOp<vector>());
-    Pstream::listCombineGather(moments, plusEqOp<vector>());
-    Pstream::broadcast(forces);
-    Pstream::broadcast(moments);
+    Pstream::listCombineAllGather(forces, plusEqOp<vector>());
+    Pstream::listCombineAllGather(moments, plusEqOp<vector>());
 
     return true;
 }
diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C
index c9eec5a15b..799c2045a0 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C
@@ -1953,8 +1953,7 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
             labelList nProcCells(distributor.countCells(distribution));
             Pout<< "Wanted distribution:" << nProcCells << endl;
 
-            Pstream::listCombineGather(nProcCells, plusEqOp<label>());
-            Pstream::broadcast(nProcCells);
+            Pstream::listCombineAllGather(nProcCells, plusEqOp<label>());
 
             Pout<< "Wanted resulting decomposition:" << endl;
             forAll(nProcCells, proci)
@@ -2210,24 +2209,25 @@ void Foam::meshRefinement::checkCoupledFaceZones(const polyMesh& mesh)
     {
         List<wordList> zoneNames(Pstream::nProcs());
         zoneNames[Pstream::myProcNo()] = fZones.names();
-        Pstream::gatherList(zoneNames);
-        Pstream::scatterList(zoneNames);
+        Pstream::allGatherList(zoneNames);
+
         // All have same data now. Check.
         forAll(zoneNames, proci)
         {
-            if (proci != Pstream::myProcNo())
+            if
+            (
+                proci != Pstream::myProcNo()
+             && (zoneNames[proci] != zoneNames[Pstream::myProcNo()])
+            )
             {
-                if (zoneNames[proci] != zoneNames[Pstream::myProcNo()])
-                {
-                    FatalErrorInFunction
-                        << "faceZones are not synchronised on processors." << nl
-                        << "Processor " << proci << " has faceZones "
-                        << zoneNames[proci] << nl
-                        << "Processor " << Pstream::myProcNo()
-                        << " has faceZones "
-                        << zoneNames[Pstream::myProcNo()] << nl
-                        << exit(FatalError);
-                }
+                FatalErrorInFunction
+                    << "faceZones are not synchronised on processors." << nl
+                    << "Processor " << proci << " has faceZones "
+                    << zoneNames[proci] << nl
+                    << "Processor " << Pstream::myProcNo()
+                    << " has faceZones "
+                    << zoneNames[Pstream::myProcNo()] << nl
+                    << exit(FatalError);
             }
         }
     }
diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBaffles.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBaffles.C
index 6ab46e3ca4..a37213e7b8 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBaffles.C
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBaffles.C
@@ -820,8 +820,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createZoneBaffles
 
         if (nTotalBaffles > 0)
         {
-            Pstream::listCombineGather(nBaffles, plusEqOp<label>());
-            Pstream::broadcast(nBaffles);
+            Pstream::listCombineAllGather(nBaffles, plusEqOp<label>());
 
             Info<< nl
                 << setf(ios_base::left)
@@ -1947,8 +1946,7 @@ void Foam::meshRefinement::findCellZoneTopo
     // - region numbers are identical on all processors
     // - keepRegion is identical ,,
     // - cellZones are identical ,,
-    Pstream::listCombineGather(regionToCellZone, maxEqOp<label>());
-    Pstream::broadcast(regionToCellZone);
+    Pstream::listCombineAllGather(regionToCellZone, maxEqOp<label>());
 
 
     // Find the region containing the keepPoint
@@ -1998,8 +1996,7 @@ void Foam::meshRefinement::findCellZoneTopo
         // - cellZones are identical ,,
         // This done at top of loop to account for geometric matching
         // not being synchronised.
-        Pstream::listCombineGather(regionToCellZone, maxEqOp<label>());
-        Pstream::broadcast(regionToCellZone);
+        Pstream::listCombineAllGather(regionToCellZone, maxEqOp<label>());
 
 
         bool changed = false;
@@ -5732,8 +5729,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
 
         // 2.Combine faceZoneNames allocated on different processors
 
-        Pstream::mapCombineGather(zonesToFaceZone, eqOp<word>());
-        Pstream::broadcast(zonesToFaceZone);
+        Pstream::mapCombineAllGather(zonesToFaceZone, eqOp<word>());
 
 
         // 3. Allocate faceZones from (now synchronised) faceZoneNames
@@ -5956,8 +5952,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
                     nPosOrientation.find(faceToConnectedZone[faceI])() += n;
                 }
             }
-            Pstream::mapCombineGather(nPosOrientation, plusEqOp<label>());
-            Pstream::broadcast(nPosOrientation);
+            Pstream::mapCombineAllGather(nPosOrientation, plusEqOp<label>());
 
 
             Info<< "Split " << nFreeStanding << " free-standing zone faces"
diff --git a/src/mesh/snappyHexMesh/refinementSurfaces/surfaceZonesInfo.C b/src/mesh/snappyHexMesh/refinementSurfaces/surfaceZonesInfo.C
index 8903cfaa2d..c70880e793 100644
--- a/src/mesh/snappyHexMesh/refinementSurfaces/surfaceZonesInfo.C
+++ b/src/mesh/snappyHexMesh/refinementSurfaces/surfaceZonesInfo.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2015 OpenFOAM Foundation
-    Copyright (C) 2015-2020 OpenCFD Ltd.
+    Copyright (C) 2015-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -485,8 +485,7 @@ Foam::labelList Foam::surfaceZonesInfo::addCellZonesToMesh
     // Check they are synced
     List<wordList> allCellZones(Pstream::nProcs());
     allCellZones[Pstream::myProcNo()] = mesh.cellZones().names();
-    Pstream::gatherList(allCellZones);
-    Pstream::scatterList(allCellZones);
+    Pstream::allGatherList(allCellZones);
 
     for (label proci = 1; proci < allCellZones.size(); proci++)
     {
@@ -646,8 +645,7 @@ Foam::labelListList Foam::surfaceZonesInfo::addFaceZonesToMesh
     // Check they are synced
     List<wordList> allFaceZones(Pstream::nProcs());
     allFaceZones[Pstream::myProcNo()] = faceZones.names();
-    Pstream::gatherList(allFaceZones);
-    Pstream::scatterList(allFaceZones);
+    Pstream::allGatherList(allFaceZones);
 
     for (label proci = 1; proci < allFaceZones.size(); proci++)
     {
diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C
index 7a16523f90..fe5f8772b9 100644
--- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C
@@ -3247,12 +3247,9 @@ void Foam::snappyRefineDriver::deleteSmallRegions
         nCellsPerRegion[regioni]++;
         nCellsPerZone[zonei]++;
     }
-    Pstream::listCombineGather(nCellsPerRegion, plusEqOp<label>());
-    Pstream::listCombineGather(regionToZone, maxEqOp<label>());
-    Pstream::listCombineGather(nCellsPerZone, plusEqOp<label>());
-    Pstream::broadcast(nCellsPerRegion);
-    Pstream::broadcast(regionToZone);
-    Pstream::broadcast(nCellsPerZone);
+    Pstream::listCombineAllGather(nCellsPerRegion, plusEqOp<label>());
+    Pstream::listCombineAllGather(regionToZone, maxEqOp<label>());
+    Pstream::listCombineAllGather(nCellsPerZone, plusEqOp<label>());
 
 
     // Mark small regions. Note that all processors have the same information
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/advancingFrontAMI/advancingFrontAMIParallelOps.C b/src/meshTools/AMIInterpolation/AMIInterpolation/advancingFrontAMI/advancingFrontAMIParallelOps.C
index 565be1efea..f199d57314 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/advancingFrontAMI/advancingFrontAMIParallelOps.C
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/advancingFrontAMI/advancingFrontAMIParallelOps.C
@@ -281,8 +281,7 @@ Foam::autoPtr<Foam::mapDistribute> Foam::advancingFrontAMI::calcProcMap
         procBb[Pstream::myProcNo()] = treeBoundBoxList();
     }
 
-    Pstream::gatherList(procBb);
-    Pstream::scatterList(procBb);
+    Pstream::allGatherList(procBb);
 
     if (debug)
     {
@@ -353,8 +352,7 @@ Foam::autoPtr<Foam::mapDistribute> Foam::advancingFrontAMI::calcProcMap
     {
         sendSizes[Pstream::myProcNo()][proci] = sendMap[proci].size();
     }
-    Pstream::gatherList(sendSizes);
-    Pstream::scatterList(sendSizes);
+    Pstream::allGatherList(sendSizes);
 
 
     // Determine order of receiving
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/nearestFaceAMI/nearestFaceAMI.C b/src/meshTools/AMIInterpolation/AMIInterpolation/nearestFaceAMI/nearestFaceAMI.C
index 9a9c52fd60..64a3243ed0 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/nearestFaceAMI/nearestFaceAMI.C
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/nearestFaceAMI/nearestFaceAMI.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020,2022 OpenCFD Ltd.
+    Copyright (C) 2020-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -50,8 +50,7 @@ Foam::autoPtr<Foam::mapDistribute> Foam::nearestFaceAMI::calcFaceMap
     List<boundBox> procBbs(Pstream::nProcs());
     procBbs[Pstream::myProcNo()] =
         boundBox(tgtPatch.points(), tgtPatch.meshPoints(), true);
-    Pstream::gatherList(procBbs);
-    Pstream::scatterList(procBbs);
+    Pstream::allGatherList(procBbs);
 
     // Identify which of my local src faces intersect with each processor
     // tgtPatch bb within the current match's search distance
diff --git a/src/meshTools/AMIInterpolation/patches/cyclicPeriodicAMI/cyclicPeriodicAMIPolyPatch/cyclicPeriodicAMIPolyPatch.C b/src/meshTools/AMIInterpolation/patches/cyclicPeriodicAMI/cyclicPeriodicAMIPolyPatch/cyclicPeriodicAMIPolyPatch.C
index 21a23c264b..fb7e1af8fb 100644
--- a/src/meshTools/AMIInterpolation/patches/cyclicPeriodicAMI/cyclicPeriodicAMIPolyPatch/cyclicPeriodicAMIPolyPatch.C
+++ b/src/meshTools/AMIInterpolation/patches/cyclicPeriodicAMI/cyclicPeriodicAMIPolyPatch/cyclicPeriodicAMIPolyPatch.C
@@ -131,13 +131,11 @@ void Foam::cyclicPeriodicAMIPolyPatch::syncTransforms() const
                 // Sync a list of separation vectors
                 List<vectorField> sep(Pstream::nProcs());
                 sep[Pstream::myProcNo()] = periodicPatch.separation();
-                Pstream::gatherList(sep);
-                Pstream::scatterList(sep);
+                Pstream::allGatherList(sep);
 
                 List<boolList> coll(Pstream::nProcs());
                 coll[Pstream::myProcNo()] = periodicPatch.collocated();
-                Pstream::gatherList(coll);
-                Pstream::scatterList(coll);
+                Pstream::allGatherList(coll);
 
                 // If locally we have zero faces pick the first one that has a
                 // separation vector
@@ -166,13 +164,11 @@ void Foam::cyclicPeriodicAMIPolyPatch::syncTransforms() const
                 // Sync a list of forward and reverse transforms
                 List<tensorField> forwardT(Pstream::nProcs());
                 forwardT[Pstream::myProcNo()] = periodicPatch.forwardT();
-                Pstream::gatherList(forwardT);
-                Pstream::scatterList(forwardT);
+                Pstream::allGatherList(forwardT);
 
                 List<tensorField> reverseT(Pstream::nProcs());
                 reverseT[Pstream::myProcNo()] = periodicPatch.reverseT();
-                Pstream::gatherList(reverseT);
-                Pstream::scatterList(reverseT);
+                Pstream::allGatherList(reverseT);
 
                 // If locally we have zero faces pick the first one that has a
                 // transformation vector
diff --git a/src/meshTools/algorithms/PointEdgeWave/PointEdgeWave.C b/src/meshTools/algorithms/PointEdgeWave/PointEdgeWave.C
index 2a6357835b..5791929cce 100644
--- a/src/meshTools/algorithms/PointEdgeWave/PointEdgeWave.C
+++ b/src/meshTools/algorithms/PointEdgeWave/PointEdgeWave.C
@@ -32,7 +32,6 @@ License
 #include "cyclicPolyPatch.H"
 #include "UIPstream.H"
 #include "UOPstream.H"
-#include "PstreamCombineReduceOps.H"
 #include "debug.H"
 #include "typeInfo.H"
 #include "globalMeshData.H"
diff --git a/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C b/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C
index cd685968c5..c3dab6608c 100644
--- a/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C
+++ b/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C
@@ -202,8 +202,7 @@ void Foam::mappedPatchBase::collectSamples
         List<pointField> globalFc(nProcs);
         globalFc[myRank] = facePoints;
 
-        Pstream::gatherList(globalFc, Pstream::msgType(), myComm);
-        Pstream::scatterList(globalFc, Pstream::msgType(), myComm);
+        Pstream::allGatherList(globalFc, Pstream::msgType(), myComm);
 
         // Rework into straight list
         patchFc = ListListOps::combine<pointField>
@@ -216,8 +215,8 @@ void Foam::mappedPatchBase::collectSamples
     {
         List<pointField> globalSamples(nProcs);
         globalSamples[myRank] = samplePoints(facePoints);
-        Pstream::gatherList(globalSamples, Pstream::msgType(), myComm);
-        Pstream::scatterList(globalSamples, Pstream::msgType(), myComm);
+        Pstream::allGatherList(globalSamples, Pstream::msgType(), myComm);
+
         // Rework into straight list
         samples = ListListOps::combine<pointField>
         (
@@ -230,8 +229,7 @@ void Foam::mappedPatchBase::collectSamples
         labelListList globalFaces(nProcs);
         globalFaces[myRank] = identity(patch_.size());
         // Distribute to all processors
-        Pstream::gatherList(globalFaces, Pstream::msgType(), myComm);
-        Pstream::scatterList(globalFaces, Pstream::msgType(), myComm);
+        Pstream::allGatherList(globalFaces, Pstream::msgType(), myComm);
 
         patchFaces = ListListOps::combine<labelList>
         (
@@ -607,14 +605,12 @@ void Foam::mappedPatchBase::findSamples
     wordList samplePatches(nProcs);
     {
         samplePatches[myRank] = samplePatch_;
-        Pstream::gatherList(samplePatches, Pstream::msgType(), myComm);
-        Pstream::scatterList(samplePatches, Pstream::msgType(), myComm);
+        Pstream::allGatherList(samplePatches, Pstream::msgType(), myComm);
     }
     wordList sampleRegions(nProcs);
     {
         sampleRegions[myRank] = sampleRegion_;
-        Pstream::gatherList(sampleRegions, Pstream::msgType(), myComm);
-        Pstream::scatterList(sampleRegions, Pstream::msgType(), myComm);
+        Pstream::allGatherList(sampleRegions, Pstream::msgType(), myComm);
     }
 
 
@@ -672,15 +668,14 @@ void Foam::mappedPatchBase::findSamples
     }
 
 
-    // Find nearest. Combine on master.
-    Pstream::listCombineGather
+    // Find nearest - globally consistent
+    Pstream::listCombineAllGather
     (
         nearest,
         nearestWorldEqOp(),
         Pstream::msgType(),
         myComm
     );
-    Pstream::broadcast(nearest, myComm);
 
     //if (debug)
     //{
diff --git a/src/meshTools/processorLOD/box/box.C b/src/meshTools/processorLOD/box/box.C
index 190edc3373..50156e4c71 100644
--- a/src/meshTools/processorLOD/box/box.C
+++ b/src/meshTools/processorLOD/box/box.C
@@ -550,8 +550,7 @@ Foam::autoPtr<Foam::mapDistribute> Foam::processorLODs::box::createLODMap
     {
         sendSizes[localProci][proci] = sendElems[proci].size();
     }
-    Pstream::gatherList(sendSizes);
-    Pstream::scatterList(sendSizes);
+    Pstream::allGatherList(sendSizes);
 
 
     // Determine order of receiving
diff --git a/src/meshTools/regionSplit2D/regionSplit2D.C b/src/meshTools/regionSplit2D/regionSplit2D.C
index 45d584860a..d11518a062 100644
--- a/src/meshTools/regionSplit2D/regionSplit2D.C
+++ b/src/meshTools/regionSplit2D/regionSplit2D.C
@@ -139,8 +139,7 @@ Foam::regionSplit2D::regionSplit2D
     // Ensure regionToCompactAddr consistent across all processors
     // - not concerned about the op (keys are unique)
     // - map size will be the number of regions in the set of faces
-    Pstream::mapCombineGather(regionToCompactAddr, minEqOp<label>());
-    Pstream::broadcast(regionToCompactAddr);
+    Pstream::mapCombineAllGather(regionToCompactAddr, minEqOp<label>());
 
     nRegions_ = regionToCompactAddr.size();
 
diff --git a/src/meshTools/topoSet/cellSources/nearestToCell/nearestToCell.C b/src/meshTools/topoSet/cellSources/nearestToCell/nearestToCell.C
index 11a97ee049..d8560bc828 100644
--- a/src/meshTools/topoSet/cellSources/nearestToCell/nearestToCell.C
+++ b/src/meshTools/topoSet/cellSources/nearestToCell/nearestToCell.C
@@ -85,8 +85,7 @@ void Foam::nearestToCell::combine(topoSet& set, const bool add) const
         );
     }
 
-    Pstream::listCombineGather(nearest, mappedPatchBase::nearestEqOp());
-    Pstream::broadcast(nearest);
+    Pstream::listCombineAllGather(nearest, mappedPatchBase::nearestEqOp());
 
     for (const auto& near : nearest)
     {
diff --git a/src/meshTools/topoSet/faceSources/regionToFace/regionToFace.C b/src/meshTools/topoSet/faceSources/regionToFace/regionToFace.C
index 960e08a2ac..8d3e28b7c8 100644
--- a/src/meshTools/topoSet/faceSources/regionToFace/regionToFace.C
+++ b/src/meshTools/topoSet/faceSources/regionToFace/regionToFace.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2017 OpenFOAM Foundation
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -177,7 +177,7 @@ void Foam::regionToFace::combine(topoSet& set, const bool add) const
     }
 
     // Globally reduce
-    combineReduce(ni, mappedPatchBase::nearestEqOp());
+    Pstream::combineAllGather(ni, mappedPatchBase::nearestEqOp());
 
     if (verbose_)
     {
diff --git a/src/meshTools/topoSet/faceZoneSources/planeToFaceZone/planeToFaceZone.C b/src/meshTools/topoSet/faceZoneSources/planeToFaceZone/planeToFaceZone.C
index ca50d80545..dfa9df9afa 100644
--- a/src/meshTools/topoSet/faceZoneSources/planeToFaceZone/planeToFaceZone.C
+++ b/src/meshTools/topoSet/faceZoneSources/planeToFaceZone/planeToFaceZone.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2020 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -153,8 +153,7 @@ void Foam::planeToFaceZone::combine(faceZoneSet& fzSet, const bool add) const
                     bitSet(),  // No border edges
                     newSetFaceRegions
                 );
-            Pstream::gatherList(procNRegions);
-            Pstream::scatterList(procNRegions);
+            Pstream::allGatherList(procNRegions);
 
             // Cumulative sum the number of regions on each processor to get an
             // offset which makes the local region ID-s globally unique
@@ -218,8 +217,7 @@ void Foam::planeToFaceZone::combine(faceZoneSet& fzSet, const bool add) const
                     regionRegions[regioni].unset(regioni);
                 }
             }
-            Pstream::listCombineGather(regionRegions, bitOrEqOp<bitSet>());
-            Pstream::broadcast(regionRegions);
+            Pstream::listCombineAllGather(regionRegions, bitOrEqOp<bitSet>());
 
             // Collapse the region connections into a map between each region
             // and the lowest numbered region that it connects to
@@ -260,8 +258,8 @@ void Foam::planeToFaceZone::combine(faceZoneSet& fzSet, const bool add) const
             {
                 ++ regionNFaces[regioni];
             }
-            Pstream::listCombineGather(regionNFaces, plusEqOp<label>());
-            Pstream::broadcast(regionNFaces);
+            Pstream::listCombineAllGather(regionNFaces, plusEqOp<label>());
+
             Info<< "    Found " << nRegions << " contiguous regions with "
                 << regionNFaces << " faces" << endl;
         }
diff --git a/src/meshTools/topoSet/pointSources/nearestToPoint/nearestToPoint.C b/src/meshTools/topoSet/pointSources/nearestToPoint/nearestToPoint.C
index 195321b7bd..d2deace74a 100644
--- a/src/meshTools/topoSet/pointSources/nearestToPoint/nearestToPoint.C
+++ b/src/meshTools/topoSet/pointSources/nearestToPoint/nearestToPoint.C
@@ -104,8 +104,7 @@ void Foam::nearestToPoint::combine(topoSet& set, const bool add) const
     }
 
 
-    Pstream::listCombineGather(nearest, mappedPatchBase::nearestEqOp());
-    Pstream::broadcast(nearest);
+    Pstream::listCombineAllGather(nearest, mappedPatchBase::nearestEqOp());
 
     for (const auto& near : nearest)
     {
diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurfacePoints/sensitivitySurfacePointsIncompressible.C b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurfacePoints/sensitivitySurfacePointsIncompressible.C
index 586ff819f2..358da1ef3c 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurfacePoints/sensitivitySurfacePointsIncompressible.C
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurfacePoints/sensitivitySurfacePointsIncompressible.C
@@ -7,7 +7,7 @@
 -------------------------------------------------------------------------------
     Copyright (C) 2007-2020 PCOpt/NTUA
     Copyright (C) 2013-2020 FOSS GP
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -715,8 +715,7 @@ void sensitivitySurfacePoints::assembleSensitivities()
                 patchScalarSens[3*ptI + 2] =
                     wallPointSensNormalVecPtr_()[patchI][ptI].z();
             }
-            Pstream::gatherList(procPatchSens);
-            Pstream::scatterList(procPatchSens);
+            Pstream::allGatherList(procPatchSens);
 
             forAll(procPatchSens, procI)
             {
diff --git a/src/optimisation/adjointOptimisation/adjoint/parameterization/NURBS/NURBS3DVolume/NURBS3DVolume/NURBS3DVolume.C b/src/optimisation/adjointOptimisation/adjoint/parameterization/NURBS/NURBS3DVolume/NURBS3DVolume/NURBS3DVolume.C
index 411d3c8ecf..54cd1a2ac6 100644
--- a/src/optimisation/adjointOptimisation/adjoint/parameterization/NURBS/NURBS3DVolume/NURBS3DVolume/NURBS3DVolume.C
+++ b/src/optimisation/adjointOptimisation/adjoint/parameterization/NURBS/NURBS3DVolume/NURBS3DVolume/NURBS3DVolume.C
@@ -995,8 +995,7 @@ Foam::vectorField Foam::NURBS3DVolume::computeControlPointSensitivities
     }
 
     // Sum contributions from all processors
-    Pstream::listCombineGather(controlPointDerivs, plusEqOp<vector>());
-    Pstream::broadcast(controlPointDerivs);
+    Pstream::listCombineAllGather(controlPointDerivs, plusEqOp<vector>());
 
     return controlPointDerivs;
 }
@@ -1081,8 +1080,7 @@ Foam::vectorField Foam::NURBS3DVolume::computeControlPointSensitivities
         }
     }
     // Sum contributions from all processors
-    Pstream::listCombineGather(controlPointDerivs, plusEqOp<vector>());
-    Pstream::broadcast(controlPointDerivs);
+    Pstream::listCombineAllGather(controlPointDerivs, plusEqOp<vector>());
 
     return controlPointDerivs;
 }
diff --git a/src/overset/cellCellStencil/cellVolumeWeight/cellVolumeWeightCellCellStencil.C b/src/overset/cellCellStencil/cellVolumeWeight/cellVolumeWeightCellCellStencil.C
index 4620c44e77..89233a7bf3 100644
--- a/src/overset/cellCellStencil/cellVolumeWeight/cellVolumeWeightCellCellStencil.C
+++ b/src/overset/cellCellStencil/cellVolumeWeight/cellVolumeWeightCellCellStencil.C
@@ -377,8 +377,7 @@ void Foam::cellCellStencils::cellVolumeWeight::findHoles
     {
         // Synchronise region status on processors
         // (could instead swap status through processor patches)
-        Pstream::listCombineGather(regionType, maxEqOp<label>());
-        Pstream::broadcast(regionType);
+        Pstream::listCombineAllGather(regionType, maxEqOp<label>());
 
         // Communicate region status through interpolative cells
         labelList cellRegionType(labelUIndList(regionType, cellRegion));
@@ -755,8 +754,7 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
     {
         nCellsPerZone[zoneID[cellI]]++;
     }
-    Pstream::listCombineGather(nCellsPerZone, plusEqOp<label>());
-    Pstream::broadcast(nCellsPerZone);
+    Pstream::listCombineAllGather(nCellsPerZone, plusEqOp<label>());
 
 
     Info<< typeName << " : detected " << nZones
diff --git a/src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.C b/src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.C
index 9d0c0fd557..5ab6e8673a 100644
--- a/src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.C
+++ b/src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2020 OpenCFD Ltd.
+    Copyright (C) 2017-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -1151,8 +1151,7 @@ void Foam::cellCellStencils::inverseDistance::findHoles
     {
         // Synchronise region status on processors
         // (could instead swap status through processor patches)
-        Pstream::listCombineGather(regionType, maxEqOp<label>());
-        Pstream::broadcast(regionType);
+        Pstream::listCombineAllGather(regionType, maxEqOp<label>());
 
         DebugInfo<< FUNCTION_NAME << " : Gathered region type" << endl;
 
@@ -1756,13 +1755,10 @@ bool Foam::cellCellStencils::inverseDistance::update()
     {
         nCellsPerZone[zoneID[cellI]]++;
     }
-    Pstream::listCombineGather(nCellsPerZone, plusEqOp<label>());
-    Pstream::broadcast(nCellsPerZone);
-
+    Pstream::listCombineAllGather(nCellsPerZone, plusEqOp<label>());
 
     const boundBox& allBb(mesh_.bounds());
 
-
     PtrList<fvMeshSubset> meshParts(nZones);
     List<treeBoundBoxList> meshBb(nZones);
 
@@ -1803,8 +1799,7 @@ bool Foam::cellCellStencils::inverseDistance::update()
             }
         }
 
-        Pstream::gatherList(procBb);
-        Pstream::scatterList(procBb);
+        Pstream::allGatherList(procBb);
 
         // Move local bounding boxes to per-mesh indexing
         forAll(meshBb, zoneI)
diff --git a/src/overset/cellCellStencil/trackingInverseDistance/trackingInverseDistanceCellCellStencil.C b/src/overset/cellCellStencil/trackingInverseDistance/trackingInverseDistanceCellCellStencil.C
index 0f0316aea6..2e163d3e37 100644
--- a/src/overset/cellCellStencil/trackingInverseDistance/trackingInverseDistanceCellCellStencil.C
+++ b/src/overset/cellCellStencil/trackingInverseDistance/trackingInverseDistanceCellCellStencil.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2020 OpenCFD Ltd.
+    Copyright (C) 2017-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -517,8 +517,7 @@ Foam::cellCellStencils::trackingInverseDistance::trackingInverseDistance
     {
         nCellsPerZone[zoneID[celli]]++;
     }
-    Pstream::listCombineGather(nCellsPerZone, plusEqOp<label>());
-    Pstream::broadcast(nCellsPerZone);
+    Pstream::listCombineAllGather(nCellsPerZone, plusEqOp<label>());
 
     meshParts_.setSize(nZones);
     forAll(meshParts_, zonei)
@@ -627,8 +626,7 @@ bool Foam::cellCellStencils::trackingInverseDistance::update()
             }
         }
 
-        Pstream::gatherList(procBb);
-        Pstream::scatterList(procBb);
+        Pstream::allGatherList(procBb);
 
         // Move local bounding boxes to per-mesh indexing
         forAll(meshBb, zonei)
diff --git a/src/overset/regionsToCell/findRefCells.C b/src/overset/regionsToCell/findRefCells.C
index 922aa3fbb8..b949b80ec1 100644
--- a/src/overset/regionsToCell/findRefCells.C
+++ b/src/overset/regionsToCell/findRefCells.C
@@ -68,8 +68,7 @@ void Foam::setRefCells
             }
         }
 
-        Pstream::listCombineGather(regionNeedReference, orEqOp<bool>());
-        Pstream::broadcast(regionNeedReference);
+        Pstream::listCombineAllGather(regionNeedReference, orEqOp<bool>());
     }
 
 
@@ -189,8 +188,7 @@ void Foam::setRefCells
             }
         }
 
-        Pstream::listCombineGather(hasRef, plusEqOp<label>());
-        Pstream::broadcast(hasRef);
+        Pstream::listCombineAllGather(hasRef, plusEqOp<label>());
 
         forAll(hasRef, regionI)
         {
diff --git a/src/parallel/decompose/decompositionMethods/multiLevelDecomp/multiLevelDecomp.C b/src/parallel/decompose/decompositionMethods/multiLevelDecomp/multiLevelDecomp.C
index 5ac913fc80..e541e49a7c 100644
--- a/src/parallel/decompose/decompositionMethods/multiLevelDecomp/multiLevelDecomp.C
+++ b/src/parallel/decompose/decompositionMethods/multiLevelDecomp/multiLevelDecomp.C
@@ -418,8 +418,12 @@ void Foam::multiLevelDecomp::decompose
             );
 
             label nPoints = returnReduce(domainPoints.size(), sumOp<label>());
-            Pstream::listCombineGather(nOutsideConnections, plusEqOp<label>());
-            Pstream::broadcast(nOutsideConnections);
+
+            Pstream::listCombineAllGather
+            (
+                nOutsideConnections,
+                plusEqOp<label>()
+            );
             label nPatches = 0;
             label nFaces = 0;
             for (const label nConnect : nOutsideConnections)
@@ -526,12 +530,11 @@ void Foam::multiLevelDecomp::decompose
                 }
 
                 reduce(nPoints, sumOp<label>());
-                Pstream::listCombineGather
+                Pstream::listCombineAllGather
                 (
                     nOutsideConnections,
                     plusEqOp<label>()
                 );
-                Pstream::broadcast(nOutsideConnections);
 
                 label nPatches = 0;
                 label nFaces = 0;
diff --git a/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C b/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C
index 3bc4dd90ff..8ebb39b778 100644
--- a/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C
+++ b/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C
@@ -323,7 +323,7 @@ Foam::word Foam::distributedTriSurfaceMesh::findLocalInstance
 bool Foam::distributedTriSurfaceMesh::read()
 {
     // Get bb of all domains.
-    procBb_.setSize(Pstream::nProcs());
+    procBb_.resize_nocopy(Pstream::nProcs());
 
     if (dict_.empty())
     {
@@ -334,8 +334,6 @@ bool Foam::distributedTriSurfaceMesh::read()
 
         procBb_[Pstream::myProcNo()] =
             treeBoundBoxList(1, treeBoundBox(localBb));
-        Pstream::gatherList(procBb_);
-        Pstream::scatterList(procBb_);
 
         dict_.add("bounds", procBb_[Pstream::myProcNo()]);
 
@@ -372,8 +370,6 @@ bool Foam::distributedTriSurfaceMesh::read()
     else
     {
         dict_.readEntry("bounds", procBb_[Pstream::myProcNo()]);
-        Pstream::gatherList(procBb_);
-        Pstream::scatterList(procBb_);
 
         // Wanted distribution type
         distType_ = distributionTypeNames_.get("distributionType", dict_);
@@ -392,6 +388,8 @@ bool Foam::distributedTriSurfaceMesh::read()
         );
     }
 
+    Pstream::allGatherList(procBb_);
+
     return true;
 }
 
@@ -4603,8 +4601,7 @@ Foam::distributedTriSurfaceMesh::localQueries
     {
         sendSizes[Pstream::myProcNo()][proci] = sendMap[proci].size();
     }
-    Pstream::gatherList(sendSizes);
-    Pstream::scatterList(sendSizes);
+    Pstream::allGatherList(sendSizes);
 
 
     // Determine receive map
@@ -4690,16 +4687,11 @@ void Foam::distributedTriSurfaceMesh::distribute
     {
         List<List<treeBoundBox>> newProcBb(Pstream::nProcs());
 
-        switch(distType_)
+        switch (distType_)
         {
             case FOLLOW:
-                newProcBb[Pstream::myProcNo()].setSize(bbs.size());
-                forAll(bbs, i)
-                {
-                    newProcBb[Pstream::myProcNo()][i] = bbs[i];
-                }
-                Pstream::gatherList(newProcBb);
-                Pstream::scatterList(newProcBb);
+                newProcBb[Pstream::myProcNo()] = bbs;
+                Pstream::allGatherList(newProcBb);
             break;
 
             case INDEPENDENT:
diff --git a/src/parallel/reconstruct/faReconstruct/faMeshReconstructor.C b/src/parallel/reconstruct/faReconstruct/faMeshReconstructor.C
index b3d64b4293..91d6636e77 100644
--- a/src/parallel/reconstruct/faReconstruct/faMeshReconstructor.C
+++ b/src/parallel/reconstruct/faReconstruct/faMeshReconstructor.C
@@ -30,7 +30,6 @@ License
 #include "globalMeshData.H"
 #include "edgeHashes.H"
 #include "Time.H"
-#include "PstreamCombineReduceOps.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -386,7 +385,7 @@ void Foam::faMeshReconstructor::calcAddressing
 
 
         // Collect from all processors
-        combineReduce
+        Pstream::combineAllGather
         (
             patchEdgeLabels,
             ListOps::appendEqOp<label>()
diff --git a/src/regionModels/regionModel/singleLayerRegion/singleLayerRegion.C b/src/regionModels/regionModel/singleLayerRegion/singleLayerRegion.C
index c3fd5afaed..fb2f837509 100644
--- a/src/regionModels/regionModel/singleLayerRegion/singleLayerRegion.C
+++ b/src/regionModels/regionModel/singleLayerRegion/singleLayerRegion.C
@@ -142,8 +142,7 @@ void Foam::regionModels::singleLayerRegion::initialise()
         }
     }
 
-    Pstream::listCombineGather(passivePatchIDs_, maxEqOp<label>());
-    Pstream::broadcast(passivePatchIDs_);
+    Pstream::listCombineAllGather(passivePatchIDs_, maxEqOp<label>());
 
     magSf.field() = 0.5*(magSf + passiveMagSf);
     magSf.correctBoundaryConditions();
diff --git a/src/sampling/meshToMesh/meshToMeshParallelOps.C b/src/sampling/meshToMesh/meshToMeshParallelOps.C
index da39fe2da7..dff4eab292 100644
--- a/src/sampling/meshToMesh/meshToMeshParallelOps.C
+++ b/src/sampling/meshToMesh/meshToMeshParallelOps.C
@@ -163,9 +163,7 @@ Foam::autoPtr<Foam::mapDistribute> Foam::meshToMesh::calcProcMap
                 procBb[Pstream::myProcNo()] = treeBoundBoxList();
             }
 
-
-            Pstream::gatherList(procBb);
-            Pstream::scatterList(procBb);
+            Pstream::allGatherList(procBb);
 
 
             if (debug)
@@ -257,8 +255,7 @@ Foam::autoPtr<Foam::mapDistribute> Foam::meshToMesh::calcProcMap
             {
                 sendSizes[Pstream::myProcNo()][proci] = sendMap[proci].size();
             }
-            Pstream::gatherList(sendSizes);
-            Pstream::scatterList(sendSizes);
+            Pstream::allGatherList(sendSizes);
 
 
             // determine order of receiving
diff --git a/src/sampling/probes/patchProbes.C b/src/sampling/probes/patchProbes.C
index 1afa7a111e..7ea65d3c0f 100644
--- a/src/sampling/probes/patchProbes.C
+++ b/src/sampling/probes/patchProbes.C
@@ -165,9 +165,8 @@ void Foam::patchProbes::findElements(const fvMesh& mesh)
     }
 
 
-    // Find nearest.
-    Pstream::listCombineGather(nearest, mappedPatchBase::nearestEqOp());
-    Pstream::broadcast(nearest);
+    // Find nearest - globally consistent
+    Pstream::listCombineAllGather(nearest, mappedPatchBase::nearestEqOp());
 
     oldPoints_.resize(this->size());
 
diff --git a/src/sampling/probes/patchProbesTemplates.C b/src/sampling/probes/patchProbesTemplates.C
index 8cc8089e5b..1666236eb7 100644
--- a/src/sampling/probes/patchProbesTemplates.C
+++ b/src/sampling/probes/patchProbesTemplates.C
@@ -111,8 +111,7 @@ Foam::patchProbes::sample(const VolumeField<Type>& vField) const
         }
     }
 
-    Pstream::listCombineGather(values, isNotEqOp<Type>());
-    Pstream::broadcast(values);
+    Pstream::listCombineAllGather(values, isNotEqOp<Type>());
 
     return tvalues;
 }
@@ -142,8 +141,7 @@ Foam::patchProbes::sample(const SurfaceField<Type>& sField) const
         }
     }
 
-    Pstream::listCombineGather(values, isNotEqOp<Type>());
-    Pstream::broadcast(values);
+    Pstream::listCombineAllGather(values, isNotEqOp<Type>());
 
     return tvalues;
 }
diff --git a/src/sampling/probes/probesTemplates.C b/src/sampling/probes/probesTemplates.C
index 9a23fee796..ea5bef4b5e 100644
--- a/src/sampling/probes/probesTemplates.C
+++ b/src/sampling/probes/probesTemplates.C
@@ -225,8 +225,7 @@ Foam::probes::sample(const VolumeField<Type>& vField) const
         }
     }
 
-    Pstream::listCombineGather(values, isNotEqOp<Type>());
-    Pstream::broadcast(values);
+    Pstream::listCombineAllGather(values, isNotEqOp<Type>());
 
     return tvalues;
 }
@@ -249,8 +248,7 @@ Foam::probes::sample(const SurfaceField<Type>& sField) const
         }
     }
 
-    Pstream::listCombineGather(values, isNotEqOp<Type>());
-    Pstream::broadcast(values);
+    Pstream::listCombineAllGather(values, isNotEqOp<Type>());
 
     return tvalues;
 }
diff --git a/src/sampling/sampledSet/cloud/cloudSet.C b/src/sampling/sampledSet/cloud/cloudSet.C
index e8e9e4f62f..48f61a1eb8 100644
--- a/src/sampling/sampledSet/cloud/cloudSet.C
+++ b/src/sampling/sampledSet/cloud/cloudSet.C
@@ -84,11 +84,8 @@ void Foam::cloudSet::calcSamples
             minFoundProc[i] = foundProc[i];
         }
     }
-    Pstream::listCombineGather(minFoundProc, minEqOp<label>());
-    Pstream::listCombineGather(maxFoundProc, maxEqOp<label>());
-
-    Pstream::broadcast(minFoundProc);
-    Pstream::broadcast(maxFoundProc);
+    Pstream::listCombineAllGather(minFoundProc, minEqOp<label>());
+    Pstream::listCombineAllGather(maxFoundProc, maxEqOp<label>());
 
 
     DynamicField<point> missingPoints(sampleCoords_.size());
diff --git a/src/sampling/sampledSet/patchCloud/patchCloudSet.C b/src/sampling/sampledSet/patchCloud/patchCloudSet.C
index 5be2459bde..f97a1848ff 100644
--- a/src/sampling/sampledSet/patchCloud/patchCloudSet.C
+++ b/src/sampling/sampledSet/patchCloud/patchCloudSet.C
@@ -154,9 +154,8 @@ void Foam::patchCloudSet::calcSamples
     }
 
 
-    // Find nearest.
-    Pstream::listCombineGather(nearest, mappedPatchBase::nearestEqOp());
-    Pstream::broadcast(nearest);
+    // Find nearest - globally consistent
+    Pstream::listCombineAllGather(nearest, mappedPatchBase::nearestEqOp());
 
 
     if (debug && Pstream::master())
diff --git a/src/sampling/sampledSet/patchSeed/patchSeedSet.C b/src/sampling/sampledSet/patchSeed/patchSeedSet.C
index c878b24ad7..0c33688a4f 100644
--- a/src/sampling/sampledSet/patchSeed/patchSeedSet.C
+++ b/src/sampling/sampledSet/patchSeed/patchSeedSet.C
@@ -177,10 +177,12 @@ void Foam::patchSeedSet::calcSamples
 
             // 2. Reduce on master. Select nearest processor.
 
-            // Find nearest. Combine on master.
-            Pstream::listCombineGather(nearest, mappedPatchBase::nearestEqOp());
-            Pstream::broadcast(nearest);
-
+            // Find nearest - globally consistent
+            Pstream::listCombineAllGather
+            (
+                nearest,
+                mappedPatchBase::nearestEqOp()
+            );
 
             // 3. Pick up my local faces that have won
 
diff --git a/src/sampling/sampledSet/sampledSets/sampledSets.C b/src/sampling/sampledSet/sampledSets/sampledSets.C
index c39e3854c1..e62a5145bb 100644
--- a/src/sampling/sampledSet/sampledSets/sampledSets.C
+++ b/src/sampling/sampledSet/sampledSets/sampledSets.C
@@ -199,12 +199,9 @@ Foam::IOobjectList Foam::sampledSets::preCheckFields(unsigned request)
         selected = mesh_.classes(fieldSelection_);
     }
 
+    // Parallel consistency (no-op in serial)
     // Probably not needed...
-    // if (Pstream::parRun())
-    // {
-    //     Pstream::mapCombineGather(selected, HashSetOps::plusEqOp<word>());
-    //     Pstream::broadcast(selected);
-    // }
+    /// Pstream::mapCombineAllGather(selected, HashSetOps::plusEqOp<word>());
 
 
     DynamicList<label> missed(fieldSelection_.size());
diff --git a/src/sampling/sampledSet/shortestPath/shortestPathSet.C b/src/sampling/sampledSet/shortestPath/shortestPathSet.C
index 11ea78d57b..69096f591b 100644
--- a/src/sampling/sampledSet/shortestPath/shortestPathSet.C
+++ b/src/sampling/sampledSet/shortestPath/shortestPathSet.C
@@ -264,7 +264,7 @@ void Foam::shortestPathSet::sync
             celli,
             Tuple2<point, bool>(origin, findMinDistance)
         );
-        Foam::combineReduce
+        Pstream::combineAllGather
         (
             searchData,
             [](ProcData& x, const ProcData& y)
diff --git a/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.C b/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.C
index 673b49cf32..8151e09754 100644
--- a/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.C
+++ b/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.C
@@ -227,8 +227,7 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
     // See which processor has the nearest. Mark and subset
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    Pstream::listCombineGather(nearest, minFirstEqOp<scalar>{});
-    Pstream::broadcast(nearest);
+    Pstream::listCombineAllGather(nearest, minFirstEqOp<scalar>{});
 
     labelList cellOrFaceLabels(fc.size(), -1);
 
diff --git a/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.C b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.C
index bcdb34fcdd..98f115ec28 100644
--- a/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.C
+++ b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.C
@@ -118,11 +118,8 @@ Foam::IOobjectList Foam::sampledSurfaces::preCheckFields()
         selected = obr_.classes(fieldSelection_);
     }
 
-    if (Pstream::parRun())
-    {
-        Pstream::mapCombineGather(selected, HashSetOps::plusEqOp<word>());
-        Pstream::broadcast(selected);
-    }
+    // Parallel consistency (no-op in serial)
+    Pstream::mapCombineAllGather(selected, HashSetOps::plusEqOp<word>());
 
 
     DynamicList<label> missed(fieldSelection_.size());
diff --git a/src/sampling/surface/isoSurface/isoSurfacePoint.C b/src/sampling/surface/isoSurface/isoSurfacePoint.C
index 7d345270bb..514dddf753 100644
--- a/src/sampling/surface/isoSurface/isoSurfacePoint.C
+++ b/src/sampling/surface/isoSurface/isoSurfacePoint.C
@@ -276,9 +276,8 @@ void Foam::isoSurfacePoint::syncUnseparatedPoints
             sharedPts[pd.sharedPointAddr()[i]] = pointValues[meshPointi];
         }
 
-        // Combine on master.
-        Pstream::listCombineGather(sharedPts, minEqOp<point>());
-        Pstream::broadcast(sharedPts);
+        // Globally consistent
+        Pstream::listCombineAllGather(sharedPts, minEqOp<point>());
 
         // Now we will all have the same information. Merge it back with
         // my local information.
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C
index 1fbee4961d..b249e9ac54 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C
@@ -867,8 +867,7 @@ Foam::scalar Foam::TDACChemistryModel<ReactionThermo, ThermoType>::solve
     if (Pstream::parRun())
     {
         List<bool> active(composition.active());
-        Pstream::listCombineGather(active, orEqOp<bool>());
-        Pstream::broadcast(active);
+        Pstream::listCombineAllGather(active, orEqOp<bool>());
 
         forAll(active, i)
         {
diff --git a/src/thermophysicalModels/radiation/radiationModels/solarLoad/faceReflecting/faceReflecting.C b/src/thermophysicalModels/radiation/radiationModels/solarLoad/faceReflecting/faceReflecting.C
index 6b7c840fb3..b7a7e98014 100644
--- a/src/thermophysicalModels/radiation/radiationModels/solarLoad/faceReflecting/faceReflecting.C
+++ b/src/thermophysicalModels/radiation/radiationModels/solarLoad/faceReflecting/faceReflecting.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -327,12 +327,10 @@ void Foam::faceReflecting::calculate()
     }
 
     // Distribute ray indexes to all proc's
-    Pstream::listCombineGather(refDisDirsIndex, maxEqOp<label>());
-    Pstream::mapCombineGather(refFacesDirIndex, minEqOp<label>());
-
     // Make sure all the processors have the same information
-    Pstream::broadcast(refDisDirsIndex);
-    Pstream::broadcast(refFacesDirIndex);
+
+    Pstream::listCombineAllGather(refDisDirsIndex, maxEqOp<label>());
+    Pstream::mapCombineAllGather(refFacesDirIndex, minEqOp<label>());
 
     scalar maxBounding = 5.0*mag(mesh_.bounds().max() - mesh_.bounds().min());
 
diff --git a/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C b/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C
index 4912800c66..7042b98f62 100644
--- a/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C
+++ b/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2016-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -500,13 +500,9 @@ void Foam::radiation::viewFactor::calculate()
             qrExt[compactGlobalIds[i]] = compactCoarseHo[i];
         }
 
-        Pstream::listCombineGather(T4, maxEqOp<scalar>());
-        Pstream::listCombineGather(E, maxEqOp<scalar>());
-        Pstream::listCombineGather(qrExt, maxEqOp<scalar>());
-
-        Pstream::broadcast(T4);
-        Pstream::broadcast(E);
-        Pstream::broadcast(qrExt);
+        Pstream::listCombineAllGather(T4, maxEqOp<scalar>());
+        Pstream::listCombineAllGather(E, maxEqOp<scalar>());
+        Pstream::listCombineAllGather(qrExt, maxEqOp<scalar>());
 
         if (Pstream::master())
         {