ENH: InterfaceCompositionModels: add new thermo types for kineticGasEvaporation

More consistent use of combineReduce, simpler and/or reductions

See merge request Development/openfoam!566

applications/solvers/multiphase/icoReactingMultiphaseInterFoam/laserDTRM/laserDTRM.C

@@ -157,12 +157,7 @@ void Foam::radiation::laserDTRM::initialiseReflection()
             );
         }
 
-        if (reflections_.size())
-        {
-            reflectionSwitch_ = true;
-        }
-
-        reflectionSwitch_ = returnReduce(reflectionSwitch_, orOp<bool>());
+        reflectionSwitch_ = returnReduceOr(reflections_.size());
     }
 }
 
@@ -299,14 +294,12 @@ void Foam::radiation::laserDTRM::initialise()
                     DTRMCloud_.addParticle(pPtr);
                 }
 
-                if (returnReduce(cellI, maxOp<label>()) == -1)
+                if (nMissed < 10 && returnReduceAnd(cellI < 0))
                 {
-                    if (++nMissed <= 10)
-                    {
-                        WarningInFunction
-                            << "Cannot find owner cell for focalPoint at "
-                            << p0 << endl;
-                    }
+                    ++nMissed;
+                    WarningInFunction
+                        << "Cannot find owner cell for focalPoint at "
+                        << p0 << endl;
                 }
             }
         }

applications/test/GAMGAgglomeration/Test-GAMGAgglomeration.C

@@ -113,7 +113,6 @@ int main(int argc, char *argv[])
         label coarseSize = max(addr)+1;
 
         Info<< "Level : " << level << endl
-            << returnReduce(addr.size(), sumOp<label>()) << endl
             << "    current size      : "
             << returnReduce(addr.size(), sumOp<label>()) << endl
             << "    agglomerated size : "

applications/test/IndirectList/Test-IndirectList.C

@@ -182,7 +182,7 @@ int main(int argc, char *argv[])
 
         // MPI barrier
         bool barrier = true;
-        Pstream::scatter(barrier);
+        Pstream::broadcast(barrier);
     }
 
 

applications/test/Tuple2/Test-Tuple2.C

@@ -138,8 +138,8 @@ int main()
             maxFirstEqOp<label>()(maxIndexed, item);
         }
 
-        Pstream::combineAllGather(minIndexed, minFirstEqOp<label>());
-        Pstream::combineAllGather(maxIndexed, maxFirstEqOp<label>());
+        Pstream::combineReduce(minIndexed, minFirstEqOp<label>());
+        Pstream::combineReduce(maxIndexed, maxFirstEqOp<label>());
 
         Info<< "Min indexed: " << minIndexed << nl
             << "Max indexed: " << maxIndexed << nl;
@@ -156,8 +156,8 @@ int main()
             maxIndexed = maxFirstOp<label>()(maxIndexed, item);
         }
 
-        Pstream::combineAllGather(minIndexed, minFirstEqOp<label>());
-        Pstream::combineAllGather(maxIndexed, maxFirstEqOp<label>());
+        Pstream::combineReduce(minIndexed, minFirstEqOp<label>());
+        Pstream::combineReduce(maxIndexed, maxFirstEqOp<label>());
 
         Info<< "Min indexed: " << minIndexed << nl
             << "Max indexed: " << maxIndexed << nl;

applications/test/decomposedBlockData/Test-decomposedBlockData.C

@@ -60,7 +60,7 @@ int main(int argc, char *argv[])
     Info<< "Reading " << file << nl << endl;
     decomposedBlockData data
     (
-        Pstream::worldComm,
+        UPstream::worldComm,
         IOobject
         (
             file,

applications/test/globalIndex/Test-globalIndex.C

@@ -204,7 +204,7 @@ int main(int argc, char *argv[])
     labelPair inOut;
     pointField allCcs(globalNumbering.gather(mesh.cellCentres()));
     inOut[0] = allCcs.size();
-    Pstream::scatter(allCcs);
+    Pstream::broadcast(allCcs);
     inOut[1] = allCcs.size();
     Pout<< "    " << inOut << endl;
 

applications/test/minMax1/Test-minMax1.C

@@ -234,11 +234,7 @@ int main(int argc, char *argv[])
 
         Pout<< "hashed: " << hashed << nl;
 
-        Pstream::mapCombineGather
-        (
-            hashed,
-            plusEqOp<scalarMinMax>()
-        );
+        Pstream::mapCombineReduce(hashed, plusEqOp<scalarMinMax>());
 
         Info<< "reduced: " << hashed << nl;
 

applications/test/parallel-communicators/Test-parallel-communicators.C

@@ -167,7 +167,7 @@ int main(int argc, char *argv[])
         (
             localValue,
             sumOp<scalar>(),
-            Pstream::msgType(),
+            UPstream::msgType(),
             comm
         );
         Pout<< "sum             :" << sum << endl;

applications/test/parallel-nonBlocking/Test-parallel-nonBlocking.C

@@ -125,13 +125,13 @@ int main(int argc, char *argv[])
     scalar data1 = 1.0;
     label request1 = -1;
     {
-        Foam::reduce(data1, sumOp<scalar>(), Pstream::msgType(), request1);
+        Foam::reduce(data1, sumOp<scalar>(), UPstream::msgType(), request1);
     }
 
     scalar data2 = 0.1;
     label request2 = -1;
     {
-        Foam::reduce(data2, sumOp<scalar>(), Pstream::msgType(), request2);
+        Foam::reduce(data2, sumOp<scalar>(), UPstream::msgType(), request2);
     }
 
 

applications/utilities/mesh/advanced/PDRMesh/PDRMesh.C

@@ -585,7 +585,7 @@ void createBaffles
 // Wrapper around find patch. Also makes sure same patch in parallel.
 label findPatch(const polyBoundaryMesh& patches, const word& patchName)
 {
-    label patchi = patches.findPatchID(patchName);
+    const label patchi = patches.findPatchID(patchName);
 
     if (patchi == -1)
     {
@@ -597,16 +597,15 @@ label findPatch(const polyBoundaryMesh& patches, const word& patchName)
 
     // Check same patch for all procs
     {
-        label newPatch = patchi;
-        reduce(newPatch, minOp<label>());
+        const label newPatchi = returnReduce(patchi, minOp<label>());
 
-        if (newPatch != patchi)
+        if (newPatchi != patchi)
         {
             FatalErrorInFunction
                 << "Patch " << patchName
                 << " should have the same patch index on all processors." << nl
                 << "On my processor it has index " << patchi
-                << " ; on some other processor it has index " << newPatch
+                << " ; on some other processor it has index " << newPatchi
                 << exit(FatalError);
         }
     }

applications/utilities/mesh/advanced/selectCells/selectCells.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -396,12 +397,11 @@ int main(int argc, char *argv[])
     meshSearch queryMesh(mesh);
 
     // Check all 'outside' points
-    forAll(outsidePts, outsideI)
+    for (const point& outsidePoint : outsidePts)
     {
-        const point& outsidePoint = outsidePts[outsideI];
+        const label celli = queryMesh.findCell(outsidePoint, -1, false);
 
-        label celli = queryMesh.findCell(outsidePoint, -1, false);
-        if (returnReduce(celli, maxOp<label>()) == -1)
+        if (returnReduceAnd(celli < 0))
         {
             FatalErrorInFunction
                 << "outsidePoint " << outsidePoint

applications/utilities/mesh/conversion/ensightToFoam/Make/files

@@ -0,0 +1,4 @@
+ensightMeshReader.C
+ensightToFoam.C
+
+EXE = $(FOAM_APPBIN)/ensightToFoam

applications/utilities/mesh/conversion/ensightToFoam/Make/options

@@ -0,0 +1,11 @@
+EXE_INC = \
+    -I$(LIB_SRC)/fileFormats/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/conversion/lnInclude
+
+EXE_LIBS = \
+    -lfileFormats \
+    -lsurfMesh \
+    -lmeshTools \
+    -lconversion

applications/utilities/mesh/conversion/ensightToFoam/ensightMeshReader.H

@@ -0,0 +1,194 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Class
+    Foam::fileFormats::ensightMeshReader
+
+Description
+
+Notes
+
+SourceFiles
+    ensightMeshReader.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_ensightMeshReader_H
+#define Foam_ensightMeshReader_H
+
+#include "meshReader.H"
+//#include "ensightReadFile.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward Declarations
+class ensightReadFile;
+
+namespace fileFormats
+{
+
+/*---------------------------------------------------------------------------*\
+                Class fileFormats::ensightMeshReader Declaration
+\*---------------------------------------------------------------------------*/
+
+class ensightMeshReader
+:
+    public meshReader
+{
+    // Private Data
+
+        //- Merge distance
+        const scalar mergeTol_;
+
+        //- Check and correct handedness
+        const bool setHandedness_;
+
+
+protected:
+
+    // Protected Data
+
+        //- mesh point to original node_id
+        labelList nodeIds_;
+
+        //- mesh cell to original element_id
+        labelList elementIds_;
+
+
+    // Protected Member Functions
+
+        //- Rotate face so lowest vertex is first
+        const face& rotateFace
+        (
+            const face& f,
+            face& rotatedFace
+        ) const;
+
+        //- Read set of vertices. Optional mapping
+        void readVerts
+        (
+            ensightReadFile& is,
+            const label nVerts,
+            const Map<label>& nodeIdToPoints,
+            DynamicList<label>& verts
+        ) const;
+
+        //- Read set of element/node IDs
+        void readIDs
+        (
+            ensightReadFile& is,
+            const bool doRead,
+            const label nShapes,
+            labelList& foamToElem,
+            Map<label>& elemToFoam
+        ) const;
+
+        //- Swap handedness of hex if needed
+        void setHandedness
+        (
+            const cellModel& model,
+            DynamicList<label>& verts,
+            const pointField& points
+        ) const;
+
+        //- Read a single part until eof (return true) or until start of next
+        //  part (return false)
+        bool readGoldPart
+        (
+            ensightReadFile& is,
+            const bool read_node_ids,
+            const bool read_elem_ids,
+
+            pointField& points,
+            labelList& pointToNodeIds,
+            Map<label>& nodeIdToPoints,
+
+            // 3D-elems : cells (cell-to-faces)
+            faceListList& cellFaces,
+            labelList& cellToElemIds,
+            Map<label>& elemIdToCells,
+
+            // 2D-elems : faces
+            faceList& faces,
+            labelList& faceToElemIDs,
+            Map<label>& elemIdToFaces
+        ) const;
+
+        //- Read the mesh from the file(s)
+        virtual bool readGeometry(const scalar scaleFactor = 1.0);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("ensightMeshReader");
+
+
+    // Constructors
+
+        //- Construct from case name
+        ensightMeshReader
+        (
+            const fileName& geomFile,
+            const objectRegistry& registry,
+            const scalar mergeTol = SMALL,
+            const scalar scaleFactor = 1.0,
+            const bool setHandedness = true
+        );
+
+
+    //- Destructor
+    virtual ~ensightMeshReader() = default;
+
+
+    // Access
+
+        //- Original node id (if supplied) or -1
+        const labelList& nodeIds() const noexcept
+        {
+            return nodeIds_;
+        }
+
+        //- Original element id (if supplied) or -1
+        const labelList& elementIds() const noexcept
+        {
+            return elementIds_;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fileFormats
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/mesh/conversion/ensightToFoam/ensightToFoam.C

@@ -0,0 +1,118 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Application
+    ensightToFoam
+
+Group
+    grpMeshConversionUtilities
+
+Description
+    Convert an Ensight Gold mesh into OpenFOAM format.
+
+Usage
+    \b ensightToFoam [OPTION] \<ensightGeometryFile\>
+
+    Options:
+      - \par -mergeTol \<factor\>
+        Specify an alternative merging tolerance as a fraction of
+        the bounding box of the points.
+
+      - \par -scale \<factor\>
+        Specify an optional geometry scaling factor.
+
+      - \par -keepHandedness
+        Do not automatically flip negative volume cells
+
+See also
+    Foam::meshReader and Foam::fileFormats::STARCDMeshReader
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "Time.H"
+#include "ensightMeshReader.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "Convert Ensight mesh to OpenFOAM"
+    );
+
+    argList::noParallel();
+    argList::addArgument(".geo file", "The file containing the geometry");
+    argList::addOption
+    (
+        "mergeTol",
+        "factor",
+        "Merge tolerance as a fraction of bounding box - 0 to disable merging"
+    );
+    argList::addOption
+    (
+        "scale",
+        "factor",
+        "Geometry scaling factor - default is 1"
+    );
+    argList::addBoolOption
+    (
+        "keepHandedness",
+        "Do not automatically flip inverted cells"
+        " (default is to do a geometric test)"
+    );
+
+    argList args(argc, argv);
+    Time runTime(args.rootPath(), args.caseName());
+
+    // Increase the precision of the points data
+    IOstream::defaultPrecision(max(10u, IOstream::defaultPrecision()));
+
+    const fileName geomFile(args.get<fileName>(1));
+
+    {
+        fileFormats::ensightMeshReader reader
+        (
+            geomFile,
+            runTime,
+            args.getOrDefault<scalar>("mergeTol", 1e-10),
+            args.getOrDefault<scalar>("scale", 1.0),
+            args.found("keepHandedness")
+        );
+
+        autoPtr<polyMesh> mesh = reader.mesh(runTime);
+        mesh().setInstance(runTime.constant());
+        mesh().write();
+    }
+
+    Info<< "\nEnd\n" << endl;
+
+    return 0;
+}
+
+// ************************************************************************* //

applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C

@@ -538,8 +538,11 @@ int main(int argc, char *argv[])
 
         // Add any patches.
 
-        label nAdded = nPatches - mesh.boundaryMesh().size();
-        reduce(nAdded, sumOp<label>());
+        const label nAdded = returnReduce
+        (
+            nPatches - mesh.boundaryMesh().size(),
+            sumOp<label>()
+        );
 
         Info<< "Adding overall " << nAdded << " processor patches." << endl;
 
@@ -946,9 +949,8 @@ int main(int argc, char *argv[])
 
         // Put all modifications into meshMod
         bool anyChange = collapser.setRefinement(allPointInfo, meshMod);
-        reduce(anyChange, orOp<bool>());
 
-        if (anyChange)
+        if (returnReduceOr(anyChange))
         {
             // Construct new mesh from polyTopoChange.
             autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false);
@@ -1118,8 +1120,7 @@ int main(int argc, char *argv[])
     processorMeshes::removeFiles(mesh);
 
     // Need writing cellSet
-    label nAdded = returnReduce(addedCellsSet.size(), sumOp<label>());
-    if (nAdded > 0)
+    if (returnReduceOr(addedCellsSet.size()))
     {
         cellSet addedCells(mesh, "addedCells", addedCellsSet);
         Info<< "Writing added cells to cellSet " << addedCells.name()

applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMesh.C

@@ -345,7 +345,7 @@ void deleteEmptyPatches(fvMesh& mesh)
             else
             {
                 // Common patch.
-                if (returnReduce(patches[patchi].empty(), andOp<bool>()))
+                if (returnReduceAnd(patches[patchi].empty()))
                 {
                     Pout<< "Deleting patch " << patchi
                         << " name:" << patches[patchi].name()
@@ -661,8 +661,8 @@ void countExtrudePatches
     }
     // Synchronise decision. Actual numbers are not important, just make
     // sure that they're > 0 on all processors.
-    Pstream::listCombineAllGather(zoneSidePatch, plusEqOp<label>());
-    Pstream::listCombineAllGather(zoneZonePatch, plusEqOp<label>());
+    Pstream::listCombineReduce(zoneSidePatch, plusEqOp<label>());
+    Pstream::listCombineReduce(zoneZonePatch, plusEqOp<label>());
 }
 
 
@@ -1848,7 +1848,7 @@ int main(int argc, char *argv[])
     const primitiveFacePatch extrudePatch(std::move(zoneFaces), mesh.points());
 
 
-    Pstream::listCombineAllGather(isInternal, orEqOp<bool>());
+    Pstream::listCombineReduce(isInternal, orEqOp<bool>());
 
     // Check zone either all internal or all external faces
     checkZoneInside(mesh, zoneNames, zoneID, extrudeMeshFaces, isInternal);
@@ -2309,7 +2309,7 @@ int main(int argc, char *argv[])
         }
 
         // Reduce
-        Pstream::mapCombineAllGather(globalSum, plusEqOp<point>());
+        Pstream::mapCombineReduce(globalSum, plusEqOp<point>());
 
         forAll(localToGlobalRegion, localRegionI)
         {

applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/DelaunayMesh/DistributedDelaunayMesh.C

@@ -588,9 +588,8 @@ Foam::label Foam::DistributedDelaunayMesh<Triangulation>::referVertices
     reduce(preInsertionSize, sumOp<label>());
     reduce(postInsertionSize, sumOp<label>());
 
-    label nTotalToInsert = referredVertices.size();
-
-    reduce(nTotalToInsert, sumOp<label>());
+    label nTotalToInsert =
+        returnReduce(referredVertices.size(), sumOp<label>());
 
     if (preInsertionSize + nTotalToInsert != postInsertionSize)
     {

applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C

@@ -217,7 +217,7 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
                     );
                 }
 
-                if (returnReduce(newCellsToRefine.size(), sumOp<label>()) == 0)
+                if (returnReduceAnd(newCellsToRefine.empty()))
                 {
                     break;
                 }
@@ -878,7 +878,7 @@ Foam::backgroundMeshDecomposition::distribute
             }
         }
 
-        if (returnReduce(cellsToRefine.size(), sumOp<label>()) == 0)
+        if (returnReduceAnd(cellsToRefine.empty()))
         {
             break;
         }

applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.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.
@@ -87,23 +87,17 @@ void Foam::conformalVoronoiMesh::cellSizeMeshOverlapsBackground() const
 
     boundBox cellSizeMeshBb = cellSizeMesh.bounds();
 
-    bool fullyContained = true;
+    bool fullyContained = cellSizeMeshBb.contains(bb);
 
-    if (!cellSizeMeshBb.contains(bb))
+    if (!fullyContained)
     {
-        Pout<< "Triangulation not fully contained in cell size mesh."
-            << endl;
-
-        Pout<< "Cell Size Mesh Bounds = " << cellSizeMesh.bounds() << endl;
-        Pout<< "foamyHexMesh Bounds         = " << bb << endl;
-
-        fullyContained = false;
+        Pout<< "Triangulation not fully contained in cell size mesh." << endl
+            << "Cell Size Mesh Bounds = " << cellSizeMeshBb << endl
+            << "foamyHexMesh Bounds         = " << bb << endl;
     }
 
-    reduce(fullyContained, andOp<unsigned int>());
-
     Info<< "Triangulation is "
-        << (fullyContained ? "fully" : "not fully")
+        << (returnReduceAnd(fullyContained) ? "fully" : "not fully")
         << " contained in the cell size mesh"
         << endl;
 }
@@ -115,12 +109,7 @@ void Foam::conformalVoronoiMesh::insertInternalPoints
     bool distribute
 )
 {
-    label nPoints = points.size();
-
-    if (Pstream::parRun())
-    {
-        reduce(nPoints, sumOp<label>());
-    }
+    const label nPoints = returnReduce(points.size(), sumOp<label>());
 
     Info<< "    " << nPoints << " points to insert..." << endl;
 
@@ -145,16 +134,15 @@ void Foam::conformalVoronoiMesh::insertInternalPoints
         map().distribute(points);
     }
 
-    label nVert = number_of_vertices();
+    label preReinsertionSize(number_of_vertices());
 
     insert(points.begin(), points.end());
 
-    label nInserted(number_of_vertices() - nVert);
-
-    if (Pstream::parRun())
-    {
-        reduce(nInserted, sumOp<label>());
-    }
+    const label nInserted = returnReduce
+    (
+        label(number_of_vertices()) - preReinsertionSize,
+        sumOp<label>()
+    );
 
     Info<< "    " << nInserted << " points inserted"
         << ", failed to insert " << nPoints - nInserted

applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C

@@ -753,7 +753,7 @@ Foam::conformalVoronoiMesh::createPolyMeshFromPoints
 
     forAll(patches, p)
     {
-        label totalPatchSize = patchDicts[p].get<label>("nFaces");
+        label nPatchFaces = patchDicts[p].get<label>("nFaces");
 
         if
         (
@@ -762,7 +762,7 @@ Foam::conformalVoronoiMesh::createPolyMeshFromPoints
         )
         {
             // Do not create empty processor patches
-            if (totalPatchSize > 0)
+            if (nPatchFaces)
             {
                 patchDicts[p].set("transform", "coincidentFullMatch");
 
@@ -781,9 +781,8 @@ Foam::conformalVoronoiMesh::createPolyMeshFromPoints
         else
         {
             // Check that the patch is not empty on every processor
-            reduce(totalPatchSize, sumOp<label>());
 
-            if (totalPatchSize > 0)
+            if (returnReduceOr(nPatchFaces))
             {
                 patches[nValidPatches] = polyPatch::New
                 (

applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C

@@ -729,7 +729,7 @@ Foam::label Foam::conformalVoronoiMesh::synchroniseSurfaceTrees
         }
     }
 
-    Pstream::listCombineAllGather(hits, plusEqOp<labelHashSet>());
+    Pstream::listCombineReduce(hits, plusEqOp<labelHashSet>());
 
     forAll(surfaceHits, eI)
     {
@@ -816,7 +816,7 @@ Foam::label Foam::conformalVoronoiMesh::synchroniseEdgeTrees
         }
     }
 
-    Pstream::listCombineAllGather(hits, plusEqOp<labelHashSet>());
+    Pstream::listCombineReduce(hits, plusEqOp<labelHashSet>());
 
     forAll(featureEdgeHits, eI)
     {

applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2017 OpenFOAM Foundation
-    Copyright (C) 2015-2018 OpenCFD Ltd.
+    Copyright (C) 2015-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -707,9 +707,7 @@ void Foam::conformalVoronoiMesh::reorderProcessorPatches
 
     Info<< incrIndent << indent << "Faces matched." << endl;
 
-    reduce(anyChanged, orOp<bool>());
-
-    if (anyChanged)
+    if (returnReduceOr(anyChanged))
     {
         label nReorderedFaces = 0;
 

applications/utilities/mesh/generation/foamyMesh/foamyHexMeshBackgroundMesh/foamyHexMeshBackgroundMesh.C

@@ -727,14 +727,12 @@ int main(int argc, char *argv[])
 
     pointField mergedPoints;
     faceList mergedFaces;
-    labelList pointMergeMap;
     PatchTools::gatherAndMerge
     (
         tolDim,
         primitivePatch(SubList<face>(isoFaces), isoPoints),
         mergedPoints,
-        mergedFaces,
-        pointMergeMap
+        mergedFaces
     );
 
     if (Pstream::master())

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C

@@ -416,10 +416,12 @@ void extractSurface
 
     // Allocate zone/patch for all patches
     HashTable<label> compactZoneID(1024);
-    forAllConstIters(patchSize, iter)
+    if (Pstream::master())
     {
-        label sz = compactZoneID.size();
-        compactZoneID.insert(iter.key(), sz);
+        forAllConstIters(patchSize, iter)
+        {
+            compactZoneID.insert(iter.key(), compactZoneID.size());
+        }
     }
     Pstream::broadcast(compactZoneID);
 
@@ -431,7 +433,7 @@ void extractSurface
         label patchi = bMesh.findPatchID(iter.key());
         if (patchi != -1)
         {
-            patchToCompactZone[patchi] = iter();
+            patchToCompactZone[patchi] = iter.val();
         }
     }
 
@@ -663,7 +665,7 @@ void removeZeroSizedPatches(fvMesh& mesh)
             if
             (
                 isA<coupledPolyPatch>(pp)
-             || returnReduce(pp.size(), sumOp<label>())
+             || returnReduceOr(pp.size())
             )
             {
                 // Coupled (and unknown size) or uncoupled and used
@@ -1889,11 +1891,8 @@ int main(int argc, char *argv[])
         );
 
         // Use the maxLocalCells from the refinement parameters
-        bool preBalance = returnReduce
-        (
-            (mesh.nCells() >= refineParams.maxLocalCells()),
-            orOp<bool>()
-        );
+        const bool preBalance =
+            returnReduceOr(mesh.nCells() >= refineParams.maxLocalCells());
 
 
         if (!overwrite && !debugLevel)

applications/utilities/mesh/manipulation/checkMesh/checkTopology.C

@@ -187,7 +187,7 @@ void Foam::checkPatch
     {
         const labelList& mp = pp.meshPoints();
 
-        if (returnReduce(mp.size(), sumOp<label>()) > 0)
+        if (returnReduceOr(mp.size()))
         {
             boundBox bb(pp.points(), mp, true); // reduce
             Info<< ' ' << bb;
@@ -252,10 +252,10 @@ Foam::label Foam::checkTopology
             }
         }
         reduce(nEmpty, sumOp<label>());
-        label nTotCells = returnReduce(mesh.cells().size(), sumOp<label>());
+        const label nCells = returnReduce(mesh.cells().size(), sumOp<label>());
 
         // These are actually warnings, not errors.
-        if (nTotCells && (nEmpty % nTotCells))
+        if (nCells && (nEmpty % nCells))
         {
             Info<< " ***Total number of faces on empty patches"
                 << " is not divisible by the number of cells in the mesh."
@@ -335,7 +335,7 @@ Foam::label Foam::checkTopology
         {
             noFailedChecks++;
 
-            label nPoints = returnReduce(points.size(), sumOp<label>());
+            const label nPoints = returnReduce(points.size(), sumOp<label>());
 
             Info<< "  <<Writing " << nPoints
                 << " unused points to set " << points.name() << endl;
@@ -472,7 +472,7 @@ Foam::label Foam::checkTopology
             }
         }
 
-        label nOneCells = returnReduce(oneCells.size(), sumOp<label>());
+        const label nOneCells = returnReduce(oneCells.size(), sumOp<label>());
 
         if (nOneCells > 0)
         {
@@ -488,7 +488,7 @@ Foam::label Foam::checkTopology
             }
         }
 
-        label nTwoCells = returnReduce(twoCells.size(), sumOp<label>());
+        const label nTwoCells = returnReduce(twoCells.size(), sumOp<label>());
 
         if (nTwoCells > 0)
         {
@@ -588,11 +588,7 @@ Foam::label Foam::checkTopology
                     }
                 }
 
-                Pstream::listCombineAllGather
-                (
-                    regionDisconnected,
-                    andEqOp<bool>()
-                );
+                Pstream::listCombineReduce(regionDisconnected, andEqOp<bool>());
             }
 
 
@@ -639,7 +635,7 @@ Foam::label Foam::checkTopology
                 cellRegions[i].write();
             }
 
-            label nPoints = returnReduce(points.size(), sumOp<label>());
+            const label nPoints = returnReduce(points.size(), sumOp<label>());
             if (nPoints)
             {
                 Info<< "  <<Writing " << nPoints

applications/utilities/mesh/manipulation/createBaffles/createBaffles.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2021 OpenCFD Ltd.
+    Copyright (C) 2016-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -138,7 +138,7 @@ void filterPatches(fvMesh& mesh, const wordHashSet& addedPatchNames)
             if
             (
                 isA<coupledPolyPatch>(pp)
-             || returnReduce(pp.size(), sumOp<label>())
+             || returnReduceOr(pp.size())
              || addedPatchNames.found(pp.name())
             )
             {

applications/utilities/mesh/manipulation/createPatch/createPatch.C

@@ -686,7 +686,7 @@ void syncPoints
 
     //- Note: hasTransformation is only used for warning messages so
     //  reduction not strictly necessary.
-    //reduce(hasTransformation, orOp<bool>());
+    //Pstream::reduceOr(hasTransformation);
 
     // Synchronize multiple shared points.
     const globalMeshData& pd = mesh.globalData();
@@ -714,7 +714,7 @@ void syncPoints
         }
 
         // Combine - globally consistent
-        Pstream::listCombineAllGather(sharedPts, cop);
+        Pstream::listCombineReduce(sharedPts, cop);
 
         // Now we will all have the same information. Merge it back with
         // my local information.

applications/utilities/mesh/manipulation/orientFaceZone/orientFaceZone.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2016 OpenFOAM Foundation
+    Copyright (C) 2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -272,7 +273,7 @@ int main(int argc, char *argv[])
         }
 
 
-        if (returnReduce(changedEdges.size(), sumOp<label>()) == 0)
+        if (returnReduceAnd(changedEdges.empty()))
         {
             break;
         }

applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2021 OpenCFD Ltd.
+    Copyright (C) 2016-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -1103,9 +1103,7 @@ int main(int argc, char *argv[])
         // Update proc maps
         if (cellProcAddressing.headerOk())
         {
-            bool localOk = (cellProcAddressing.size() == mesh.nCells());
-
-            if (returnReduce(localOk, andOp<bool>()))
+            if (returnReduceAnd(cellProcAddressing.size() == mesh.nCells()))
             {
                 Info<< "Renumbering processor cell decomposition map "
                     << cellProcAddressing.name() << endl;
@@ -1129,9 +1127,7 @@ int main(int argc, char *argv[])
 
         if (faceProcAddressing.headerOk())
         {
-            bool localOk = (faceProcAddressing.size() == mesh.nFaces());
-
-            if (returnReduce(localOk, andOp<bool>()))
+            if (returnReduceAnd(faceProcAddressing.size() == mesh.nFaces()))
             {
                 Info<< "Renumbering processor face decomposition map "
                     << faceProcAddressing.name() << endl;
@@ -1171,9 +1167,7 @@ int main(int argc, char *argv[])
 
         if (pointProcAddressing.headerOk())
         {
-            bool localOk = (pointProcAddressing.size() == mesh.nPoints());
-
-            if (returnReduce(localOk, andOp<bool>()))
+            if (returnReduceAnd(pointProcAddressing.size() == mesh.nPoints()))
             {
                 Info<< "Renumbering processor point decomposition map "
                     << pointProcAddressing.name() << endl;
@@ -1197,12 +1191,13 @@ int main(int argc, char *argv[])
 
         if (boundaryProcAddressing.headerOk())
         {
-            bool localOk =
+            if
             (
-                boundaryProcAddressing.size()
-             == mesh.boundaryMesh().size()
-            );
-            if (returnReduce(localOk, andOp<bool>()))
+                returnReduceAnd
+                (
+                    boundaryProcAddressing.size() == mesh.boundaryMesh().size()
+                )
+            )
             {
                 // No renumbering needed
             }

applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C

@@ -847,7 +847,7 @@ void createAndWriteRegion
 
         if (!isA<processorPolyPatch>(pp))
         {
-            if (returnReduce(pp.size(), sumOp<label>()) > 0)
+            if (returnReduceOr(pp.size()))
             {
                 oldToNew[patchi] = newI;
                 if (!addedPatches.found(patchi))
@@ -1114,7 +1114,7 @@ label findCorrespondingRegion
         }
     }
 
-    Pstream::listCombineAllGather(cellsInZone, plusEqOp<label>());
+    Pstream::listCombineReduce(cellsInZone, plusEqOp<label>());
 
     // Pick region with largest overlap of zoneI
     label regionI = findMax(cellsInZone);

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());
-    Pstream::combineAllGather(classNames, uniqueEqOp<word>());
+    Pstream::combineReduce(classNames, uniqueEqOp<word>());
 
     // Check for correct type
     if (classNames[0] == T::typeName)
@@ -429,7 +429,7 @@ int main(int argc, char *argv[])
             )
         );
 
-        Pstream::combineAllGather(lagrangianDirs, uniqueEqOp<fileName>());
+        Pstream::combineReduce(lagrangianDirs, uniqueEqOp<fileName>());
 
         if (!lagrangianDirs.empty())
         {
@@ -466,7 +466,7 @@ int main(int argc, char *argv[])
                 )
             );
 
-            Pstream::combineAllGather(cloudDirs, uniqueEqOp<fileName>());
+            Pstream::combineReduce(cloudDirs, uniqueEqOp<fileName>());
 
             forAll(cloudDirs, i)
             {
@@ -492,7 +492,7 @@ int main(int argc, char *argv[])
                 );
 
                 // Combine with all other cloud objects
-                Pstream::combineAllGather(cloudFields, uniqueEqOp<word>());
+                Pstream::combineReduce(cloudFields, uniqueEqOp<word>());
 
                 for (const word& name : cloudFields)
                 {

applications/utilities/parallelProcessing/reconstructParMesh/reconstructParMesh.C

@@ -252,7 +252,7 @@ autoPtr<mapPolyMesh> mergeSharedPoints
     Info<< "mergeSharedPoints : detected " << pointToMaster.size()
         << " points that are to be merged." << endl;
 
-    if (returnReduce(pointToMaster.size(), sumOp<label>()) == 0)
+    if (returnReduceAnd(pointToMaster.empty()))
     {
         return nullptr;
     }

applications/utilities/parallelProcessing/redistributePar/parLagrangianDistributor.C

@@ -89,10 +89,11 @@ void Foam::parLagrangianDistributor::findClouds
     }
 
     // Synchronise cloud names
-    Pstream::combineGather(cloudNames, ListOps::uniqueEqOp<word>());
-    Pstream::broadcast(cloudNames);
+    Pstream::combineReduce(cloudNames, ListOps::uniqueEqOp<word>());
+    Foam::sort(cloudNames);  // Consistent order
 
-    objectNames.setSize(cloudNames.size());
+    objectNames.clear();
+    objectNames.resize(cloudNames.size());
 
     for (const fileName& localCloudName : localCloudDirs)
     {
@@ -124,11 +125,11 @@ void Foam::parLagrangianDistributor::findClouds
         }
     }
 
-    // Synchronise objectNames
-    forAll(objectNames, i)
+    // Synchronise objectNames (per cloud)
+    for (wordList& objNames : objectNames)
     {
-        Pstream::combineGather(objectNames[i], ListOps::uniqueEqOp<word>());
-        Pstream::broadcast(objectNames[i]);
+        Pstream::combineReduce(objNames, ListOps::uniqueEqOp<word>());
+        Foam::sort(objNames);  // Consistent order
     }
 }
 
@@ -291,7 +292,7 @@ Foam::parLagrangianDistributor::distributeLagrangianPositions
         nsTransPs[sendProcI] = subMap[sendProcI].size();
     }
     // Send sizes across. Note: blocks.
-    Pstream::combineAllGather(sizes, Pstream::listEq());
+    Pstream::combineReduce(sizes, Pstream::listEq());
 
     labelListList constructMap(Pstream::nProcs());
     label constructSize = 0;

applications/utilities/parallelProcessing/redistributePar/parLagrangianDistributorTemplates.C

@@ -51,14 +51,9 @@ Foam::wordList Foam::parLagrangianDistributor::filterObjects
       : objects.names<Container>(selectedFields)
     );
 
-    // Parallel synchronise
-    // - Combine names from all processors
-
-    Pstream::combineGather(fieldNames, ListOps::uniqueEqOp<word>());
-    Pstream::broadcast(fieldNames);
-
-    // Sort for consistent order on all processors
-    Foam::sort(fieldNames);
+    // Parallel synchronise - combine names from all processors
+    Pstream::combineReduce(fieldNames, ListOps::uniqueEqOp<word>());
+    Foam::sort(fieldNames);  // Consistent order
 
     return fieldNames;
 }
@@ -93,11 +88,8 @@ Foam::label Foam::parLagrangianDistributor::distributeFields
         if (!nFields)
         {
             // Performing an all-to-one (reconstruct)?
-            reconstruct = returnReduce
-            (
-                (!map.constructSize() || Pstream::master()),
-                andOp<bool>()
-            );
+            reconstruct =
+                returnReduceAnd(!map.constructSize() || Pstream::master());
         }
 
         if (verbose_)
@@ -200,11 +192,8 @@ Foam::label Foam::parLagrangianDistributor::distributeFieldFields
         if (!nFields)
         {
             // Performing an all-to-one (reconstruct)?
-            reconstruct = returnReduce
-            (
-                (!map.constructSize() || Pstream::master()),
-                andOp<bool>()
-            );
+            reconstruct =
+                returnReduceAnd(!map.constructSize() || Pstream::master());
         }
 
         if (verbose_)
@@ -349,11 +338,8 @@ Foam::label Foam::parLagrangianDistributor::distributeStoredFields
         if (!nFields)
         {
             // Performing an all-to-one (reconstruct)?
-            reconstruct = returnReduce
-            (
-                (!map.constructSize() || Pstream::master()),
-                andOp<bool>()
-            );
+            reconstruct =
+                returnReduceAnd(!map.constructSize() || Pstream::master());
         }
 
         if (verbose_)

applications/utilities/parallelProcessing/redistributePar/redistributePar.C

@@ -132,7 +132,7 @@ void createTimeDirs(const fileName& path)
         //Pstream::parRun(oldParRun);  // Restore parallel state
         masterTimeDirs = localTimeDirs;
     }
-    Pstream::scatter(masterTimeDirs);
+    Pstream::broadcast(masterTimeDirs);
     //DebugVar(masterTimeDirs);
     //DebugVar(localTimeDirs);
 
@@ -1168,7 +1168,7 @@ int main(int argc, char *argv[])
     bool nfs = true;
     {
         List<fileName> roots(1, args.rootPath());
-        Pstream::combineAllGather(roots, ListOps::uniqueEqOp<fileName>());
+        Pstream::combineReduce(roots, ListOps::uniqueEqOp<fileName>());
         nfs = (roots.size() == 1);
     }
 
@@ -1203,7 +1203,7 @@ int main(int argc, char *argv[])
         }
     }
     // If master changed to decompose mode make sure all nodes know about it
-    Pstream::scatter(decompose);
+    Pstream::broadcast(decompose);
 
 
     // If running distributed we have problem of new processors not finding
@@ -1288,7 +1288,7 @@ int main(int argc, char *argv[])
         // use the times list from the master processor
         // and select a subset based on the command-line options
         instantList timeDirs = timeSelector::select(runTime.times(), args);
-        Pstream::scatter(timeDirs);
+        Pstream::broadcast(timeDirs);
 
         if (timeDirs.empty())
         {
@@ -1547,7 +1547,7 @@ int main(int argc, char *argv[])
                 if
                 (
                     !volMeshHaveUndecomposed
-                 || !returnReduce(haveVolAddressing, andOp<bool>())
+                 || !returnReduceAnd(haveVolAddressing)
                 )
                 {
                     Info<< "No undecomposed mesh. Creating from: "
@@ -1615,7 +1615,7 @@ int main(int argc, char *argv[])
                  &&
                     (
                         !areaMeshHaveUndecomposed
-                     || !returnReduce(haveAreaAddressing, andOp<bool>())
+                     || !returnReduceAnd(haveAreaAddressing)
                     )
                 )
                 {
@@ -2067,7 +2067,7 @@ int main(int argc, char *argv[])
                 args
             )[0].value();
         }
-        Pstream::scatter(masterTime);
+        Pstream::broadcast(masterTime);
         Info<< "Setting time to that of master or undecomposed case : "
             << masterTime << endl;
         runTime.setTime(masterTime, 0);

applications/utilities/postProcessing/dataConversion/foamToEnsight/checkMeshMoving.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
+    Copyright (C) 2021-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@@ -78,6 +78,6 @@ if (timeDirs.size() > 1)
 }
 
 // Ensure consistency
-reduce(hasMovingMesh, orOp<bool>());
+Pstream::reduceOr(hasMovingMesh);
 
 // ************************************************************************* //

applications/utilities/postProcessing/dataConversion/foamToEnsight/convertLagrangian.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2021 OpenCFD Ltd.
+    Copyright (C) 2018-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@@ -38,7 +38,7 @@ if (doLagrangian)
         Info<< "Cloud " << cloudName << " (";
 
         const bool cloudExists =
-            returnReduce(currentCloudDirs.found(cloudName), orOp<bool>());
+            returnReduceOr(currentCloudDirs.found(cloudName));
 
         {
             autoPtr<ensightFile> os = ensCase.newCloud(cloudName);
@@ -82,7 +82,7 @@ if (doLagrangian)
                 const bool oldParRun = Pstream::parRun(false);
                 fieldExists = fieldObject.typeHeaderOk<IOField<scalar>>(false);
                 Pstream::parRun(oldParRun);  // Restore parallel state
-                reduce(fieldExists, orOp<bool>());
+                Pstream::reduceOr(fieldExists);
             }
 
             bool wrote = false;

applications/utilities/postProcessing/dataConversion/foamToEnsight/findCloudFields.H

@@ -101,7 +101,7 @@ if (timeDirs.size() && doLagrangian)
     {
         for (auto& cloudFields : regionCloudFields)
         {
-            Pstream::mapCombineAllGather
+            Pstream::mapCombineReduce
             (
                 cloudFields,
                 HashTableOps::plusEqOp<word>()

applications/utilities/postProcessing/dataConversion/foamToEnsight/readFields.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -64,9 +64,7 @@ Foam::label Foam::checkData
             }
         }
 
-        reduce(good, andOp<bool>());
-
-        if (good)
+        if (returnReduceAnd(good))
         {
             goodFields.insert(fieldName);
         }

applications/utilities/postProcessing/dataConversion/foamToVTK/convertLagrangian.H

@@ -46,11 +46,9 @@ if (doLagrangian)
     if (Pstream::parRun())
     {
         // Synchronise cloud names
-        Pstream::combineGather(cloudNames, ListOps::uniqueEqOp<word>());
-        Pstream::broadcast(cloudNames);
+        Pstream::combineReduce(cloudNames, ListOps::uniqueEqOp<word>());
     }
-    // Consistent order
-    Foam::sort(cloudNames);
+    Foam::sort(cloudNames);  // Consistent order
 
     for (const word& cloudName : cloudNames)
     {
@@ -66,7 +64,7 @@ if (doLagrangian)
             isCloud = true;
         }
 
-        if (!returnReduce(isCloud, orOp<bool>()))
+        if (!returnReduceOr(isCloud))
         {
             continue;
         }
@@ -78,7 +76,7 @@ if (doLagrangian)
         );
 
         // Are there cloud fields (globally)?
-        if (returnReduce(cloudObjs.empty(), andOp<bool>()))
+        if (returnReduceAnd(cloudObjs.empty()))
         {
             continue;
         }

applications/utilities/postProcessing/lagrangian/particleTracks/particleTracks.C

@@ -186,7 +186,7 @@ int main(int argc, char *argv[])
         const label maxNProcs = returnReduce(maxIds.size(), maxOp<label>());
         maxIds.resize(maxNProcs, -1);
 
-        Pstream::listCombineAllGather(maxIds, maxEqOp<label>());
+        Pstream::listCombineReduce(maxIds, maxEqOp<label>());
 
         // From ids to count
         const labelList numIds = maxIds + 1;

applications/utilities/postProcessing/lagrangian/particleTracks/particleTracksSamplerTemplates.C

@@ -77,8 +77,7 @@ Foam::label Foam::particleTracksSampler::setTrackFields
 
         if (Pstream::parRun())
         {
-            Pstream::combineGather(fieldNames, ListOps::uniqueEqOp<word>());
-            Pstream::broadcast(fieldNames);
+            Pstream::combineReduce(fieldNames, ListOps::uniqueEqOp<word>());
         }
 
         for (const word& fieldName : fieldNames)

applications/utilities/postProcessing/lagrangian/steadyParticleTracks/createControls.H

@@ -15,10 +15,5 @@ wordRes acceptFields(propsDict.get<wordRes>("fields"));
 wordRes excludeFields;
 propsDict.readIfPresent("exclude", excludeFields);
 
-const dictionary formatOptions
-(
-    propsDict.subOrEmptyDict("formatOptions", keyType::LITERAL)
-);
-
 
 // ************************************************************************* //

applications/utilities/postProcessing/miscellaneous/postChannel/channelIndex.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -68,7 +69,7 @@ void Foam::channelIndex::walkOppositeFaces
         blockedFace[facei] = true;
     }
 
-    while (returnReduce(frontFaces.size(), sumOp<label>()) > 0)
+    while (returnReduceOr(frontFaces.size()))
     {
         // Transfer across.
         boolList isFrontBndFace(nBnd, false);

applications/utilities/postProcessing/miscellaneous/postChannel/channelIndexTemplates.C

@@ -40,7 +40,7 @@ Foam::Field<T> Foam::channelIndex::regionSum(const Field<T>& cellField) const
     }
 
     // Global sum
-    Pstream::listCombineAllGather(regionField, plusEqOp<T>());
+    Pstream::listCombineReduce(regionField, plusEqOp<T>());
 
     return regionField;
 }

applications/utilities/preProcessing/foamUpgradeCyclics/foamUpgradeCyclics.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -238,7 +238,7 @@ void rewriteBoundary
 
     patches.reorder(oldToNew);
 
-    if (returnReduce(nOldCyclics, sumOp<label>()) > 0)
+    if (returnReduceOr(nOldCyclics))
     {
         if (dryrun)
         {
@@ -299,7 +299,7 @@ void rewriteField
 
     dictionary& boundaryField = fieldDict.subDict("boundaryField");
 
-    label nChanged = 0;
+    bool hasChange = false;
 
     forAllConstIters(thisNames, iter)
     {
@@ -337,13 +337,13 @@ void rewriteField
             );
             Info<< "    Adding entry " << nbrNames[patchName] << endl;
 
-            nChanged++;
+            hasChange = true;
         }
     }
 
     //Info<< "New boundaryField:" << boundaryField << endl;
 
-    if (returnReduce(nChanged, sumOp<label>()) > 0)
+    if (returnReduceOr(hasChange))
     {
         if (dryrun)
         {

applications/utilities/preProcessing/mapFieldsPar/mapLagrangian.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2018-2019 OpenCFD Ltd.
+    Copyright (C) 2018-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -113,10 +113,9 @@ void mapLagrangian(const meshToMesh& interp)
 
         if
         (
-            returnReduce
+            returnReduceOr
             (
-                (objects.found("coordinates") || objects.found("positions")),
-                orOp<bool>()
+                objects.found("coordinates") || objects.found("positions")
             )
         )
         {

applications/utilities/preProcessing/viewFactorsGen/shootRays.H

@@ -85,7 +85,7 @@ for (const int proci : Pstream::allProcs())
             }
         }
 
-    } while (returnReduce(i < myFc.size(), orOp<bool>()));
+    } while (returnReduceOr(i < myFc.size()));
 
     List<pointIndexHit> hitInfo(startIndex.size());
     surfacesMesh.findLine(start, end, hitInfo);
@@ -176,7 +176,7 @@ for (const int proci : Pstream::allProcs())
         firstLoop = false;
         iter ++;
 
-    } while (returnReduce(hitInfo.size(), orOp<bool>()) > 0 && iter < 10);
+    } while (returnReduceOr(hitInfo.size()) && iter < 10);
 
     startIndex.clear();
     end.clear();

applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C

@@ -531,15 +531,10 @@ int main(int argc, char *argv[])
         nFineFaces += patches[patchi].size();
     }
 
-    // total number of coarse faces
-    label totalNCoarseFaces = nCoarseFaces;
 
-    reduce(totalNCoarseFaces, sumOp<label>());
-
-    if (Pstream::master())
-    {
-        Info << "\nTotal number of coarse faces: "<< totalNCoarseFaces << endl;
-    }
+    Info<< "\nTotal number of coarse faces: "
+        << returnReduce(nCoarseFaces, sumOp<label>())
+        << endl;
 
     if (Pstream::master() && debug)
     {
@@ -829,15 +824,11 @@ int main(int argc, char *argv[])
         nCoarseFaces
     );
 
-    label totalPatches = coarsePatches.size();
-    reduce(totalPatches, maxOp<label>());
+    const label totalPatches =
+        returnReduce(coarsePatches.size(), maxOp<label>());
 
     // Matrix sum in j(Fij) for each i (if enclosure sum = 1)
-    scalarSquareMatrix sumViewFactorPatch
-    (
-        totalPatches,
-        0.0
-    );
+    scalarSquareMatrix sumViewFactorPatch(totalPatches, Zero);
 
     scalarList patchArea(totalPatches, Zero);
 

applications/utilities/surface/surfaceMeshExtract/surfaceMeshExtract.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2017-2021 OpenCFD Ltd.
+    Copyright (C) 2017-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -302,17 +302,19 @@ int main(int argc, char *argv[])
             Pstream::mapCombineGather(patchSize, plusEqOp<label>());
             Pstream::mapCombineGather(zoneSize, plusEqOp<label>());
 
-            // Allocate compact numbering for all patches/faceZones
-            forAllConstIters(patchSize, iter)
+            if (Pstream::master())
             {
-                compactZoneID.insert(iter.key(), compactZoneID.size());
-            }
+                // Allocate compact numbering for all patches/faceZones
+                forAllConstIters(patchSize, iter)
+                {
+                    compactZoneID.insert(iter.key(), compactZoneID.size());
+                }
 
-            forAllConstIters(zoneSize, iter)
-            {
-                compactZoneID.insert(iter.key(), compactZoneID.size());
+                forAllConstIters(zoneSize, iter)
+                {
+                    compactZoneID.insert(iter.key(), compactZoneID.size());
+                }
             }
-
             Pstream::broadcast(compactZoneID);
 
 
@@ -391,13 +393,9 @@ int main(int argc, char *argv[])
         // Gather all faces
         List<faceList> gatheredFaces(Pstream::nProcs());
         gatheredFaces[Pstream::myProcNo()] = allBoundary.localFaces();
-        forAll(gatheredFaces[Pstream::myProcNo()], i)
+        for (face& f : gatheredFaces[Pstream::myProcNo()])
         {
-            inplaceRenumber
-            (
-                pointToGlobal,
-                gatheredFaces[Pstream::myProcNo()][i]
-            );
+            inplaceRenumber(pointToGlobal, f);
         }
         Pstream::gatherList(gatheredFaces);
 

applications/utilities/surface/surfaceRedistributePar/surfaceRedistributePar.C

@@ -252,7 +252,7 @@ int main(int argc, char *argv[])
         dict.add("bounds", bbs);
 
         // Scatter patch information
-        Pstream::scatter(s.patches());
+        Pstream::broadcast(s.patches());
 
         // Construct distributedTrisurfaceMesh from components
         IOobject notReadIO(io);

src/OpenFOAM/containers/Lists/SortableList/ParSortableList.C

@@ -202,7 +202,7 @@ void Foam::ParSortableList<Type>::sort()
 
         getPivots(sortedPivots, pivots);
     }
-    Pstream::scatter(pivots);
+    Pstream::broadcast(pivots);
 
     if (debug)
     {

src/OpenFOAM/db/IOobjectList/IOobjectList.C

@@ -90,12 +90,10 @@ void Foam::IOobjectList::syncNames(wordList& objNames)
     if (Pstream::parRun())
     {
         // Synchronize names
-        Pstream::combineGather(objNames, ListOps::uniqueEqOp<word>());
-        Pstream::broadcast(objNames);
+        Pstream::combineReduce(objNames, ListOps::uniqueEqOp<word>());
     }
 
-    // Consistent order on all processors
-    Foam::sort(objNames);
+    Foam::sort(objNames);  // Consistent order
 }
 
 

src/OpenFOAM/db/IOstreams/Pstreams/PstreamBuffers.C

@@ -482,7 +482,7 @@ bool Foam::PstreamBuffers::finishedSends
         }
     }
 
-    reduce(changed, orOp<bool>());
+    UPstream::reduceOr(changed);
 
     if (changed)
     {

src/OpenFOAM/db/IOstreams/Pstreams/UOPstream.H

@@ -327,7 +327,7 @@ public:
         //  \return True on success
         static bool write
         (
-            const commsTypes commsType,
+            const UPstream::commsTypes commsType,
             const int toProcNo,
             const char* buf,
             const std::streamsize bufSize,

src/OpenFOAM/db/IOstreams/Pstreams/UPstream.C

@@ -434,7 +434,7 @@ Foam::UPstream::commsTypes Foam::UPstream::defaultCommsType
 
 namespace Foam
 {
-    // Register re-reader
+    //- Registered reader for UPstream::defaultCommsType
     class addcommsTypeToOpt
     :
         public ::Foam::simpleRegIOobject

src/OpenFOAM/db/IOstreams/Pstreams/UPstream.H

@@ -66,12 +66,12 @@ public:
     //- Types of communications
     enum class commsTypes : char
     {
-        blocking,       //!< "blocking"
-        scheduled,      //!< "scheduled"
-        nonBlocking     //!< "nonBlocking"
+        blocking,       //!< "blocking"    : (MPI_Bsend, MPI_Recv)
+        scheduled,      //!< "scheduled"   : (MPI_Send, MPI_Recv)
+        nonBlocking     //!< "nonBlocking" : (MPI_Isend, MPI_Irecv)
     };
 
-    //- Names of the communication types
+    //- Enumerated names for the communication types
     static const Enum<commsTypes> commsTypeNames;
 
 

src/OpenFOAM/db/dictionary/functionEntries/codeStream/codeStream.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2018-2021 OpenCFD Ltd.
+    Copyright (C) 2018-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -317,7 +317,7 @@ Foam::functionEntries::codeStream::getFunction
     bool haveLib = lib;
     if (!doingMasterOnlyReading(topDict))
     {
-        reduce(haveLib, andOp<bool>());
+        Pstream::reduceAnd(haveLib);
     }
 
     if (!haveLib)

src/OpenFOAM/db/dynamicLibrary/codedBase/codedBase.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2021 OpenCFD Ltd.
+    Copyright (C) 2016-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -332,7 +332,8 @@ void Foam::codedBase::createLibrary
             << " after waiting: have masterSize:" << masterSize
             << " and localSize:" << mySize << endl;
     }
-    reduce(create, orOp<bool>());
+
+    Pstream::reduceOr(create);  // MPI barrier
 }
 
 

src/OpenFOAM/db/functionObjects/functionObjectList/functionObjectList.C

@@ -1171,7 +1171,7 @@ bool Foam::functionObjectList::read()
                 FatalIOError.throwing(oldThrowingIOerr);
 
                 // Require valid functionObject on all processors
-                if (!returnReduce(bool(objPtr), andOp<bool>()))
+                if (!returnReduceAnd(bool(objPtr)))
                 {
                     objPtr.reset(nullptr);
                     ok = false;

src/OpenFOAM/expressions/exprDriver/exprDriverTemplates.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2010-2018 Bernhard Gschaider
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -309,8 +309,7 @@ bool Foam::expressions::exprDriver::isLocalVariable
         // Do size checking if requested
         if (good && expectedSize >= 0)
         {
-            good = (var.size() == expectedSize);
-            reduce(good, andOp<bool>());
+            good = returnReduceAnd(var.size() == expectedSize);
 
             if (debug && !good)
             {

src/OpenFOAM/expressions/fields/fieldExprDriverTemplates.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -52,7 +52,7 @@ Foam::expressions::fieldExpr::parseDriver::getField
 
         const label len = (hasPointData ? this->pointSize() : this->size());
 
-        if (returnReduce((vals.size() == len), andOp<bool>()))
+        if (returnReduceAnd(vals.size() == len))
         {
             // Return a copy of the field
             return tmp<Field<Type>>::New(vals);

src/OpenFOAM/fields/Fields/Field/FieldFunctions.C

@@ -579,7 +579,7 @@ template<class Type>                                                           \
 ReturnType gFunc(const UList<Type>& f, const label comm)                       \
 {                                                                              \
     ReturnType res = Func(f);                                                  \
-    reduce(res, rFunc##Op<ReturnType>(), Pstream::msgType(), comm);            \
+    reduce(res, rFunc##Op<ReturnType>(), UPstream::msgType(), comm);           \
     return res;                                                                \
 }                                                                              \
 TMP_UNARY_FUNCTION(ReturnType, gFunc)
@@ -611,7 +611,7 @@ typename scalarProduct<Type, Type>::type gSumProd
     typedef typename scalarProduct<Type, Type>::type prodType;
 
     prodType result = sumProd(f1, f2);
-    reduce(result, sumOp<prodType>(), Pstream::msgType(), comm);
+    reduce(result, sumOp<prodType>(), UPstream::msgType(), comm);
     return result;
 }
 
@@ -624,7 +624,7 @@ Type gSumCmptProd
 )
 {
     Type SumProd = sumCmptProd(f1, f2);
-    reduce(SumProd, sumOp<Type>(), Pstream::msgType(), comm);
+    reduce(SumProd, sumOp<Type>(), UPstream::msgType(), comm);
     return SumProd;
 }
 
@@ -637,7 +637,7 @@ Type gAverage
 {
     label n = f.size();
     Type s = sum(f);
-    sumReduce(s, n, Pstream::msgType(), comm);
+    sumReduce(s, n, UPstream::msgType(), comm);
 
     if (n > 0)
     {

src/OpenFOAM/fields/Fields/Field/FieldOps.C

@@ -172,7 +172,7 @@ Foam::Tuple2<T1,T2> Foam::FieldOps::findMinData
         result.second() = data[i];
     }
 
-    Pstream::combineAllGather(result, minFirstEqOp<T1>());
+    Pstream::combineReduce(result, minFirstEqOp<T1>());
     return result;
 }
 
@@ -193,7 +193,7 @@ Foam::Tuple2<T1,T2> Foam::FieldOps::findMaxData
         result.second() = data[i];
     }
 
-    Pstream::combineAllGather(result, maxFirstEqOp<T1>());
+    Pstream::combineReduce(result, maxFirstEqOp<T1>());
     return result;
 }
 

src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C

@@ -1077,18 +1077,16 @@ bool Foam::GeometricField<Type, PatchField, GeoMesh>::needReference() const
 
     bool needRef = true;
 
-    forAll(boundaryField_, patchi)
+    for (const auto& pf : boundaryField_)
     {
-        if (boundaryField_[patchi].fixesValue())
+        if (pf.fixesValue())
         {
             needRef = false;
             break;
         }
     }
 
-    reduce(needRef, andOp<bool>());
-
-    return needRef;
+    return returnReduceAnd(needRef);
 }
 
 

src/OpenFOAM/global/fileOperations/masterUncollatedFileOperation/masterUncollatedFileOperation.C

@@ -2032,7 +2032,7 @@ Foam::fileOperations::masterUncollatedFileOperation::readStream
             }
             // Reduce (not broadcast)
             // - if we have multiple master files (FUTURE)
-            reduce(bigSize, orOp<bool>());  //, UPstream::msgType(), comm_);
+            Pstream::reduceOr(bigSize);  //, comm_);
 
             const UPstream::commsTypes myCommsType
             (

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGProcAgglomerations/procFacesGAMGProcAgglomeration/procFacesGAMGProcAgglomeration.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -196,7 +196,7 @@ Foam::procFacesGAMGProcAgglomeration::processorAgglomeration
         fineToCoarse = labelUIndList(oldToNew, fineToCoarse)();
     }
 
-    Pstream::scatter(fineToCoarse, Pstream::msgType(), mesh.comm());
+    Pstream::broadcast(fineToCoarse, mesh.comm());
     UPstream::freeCommunicator(singleCellMeshComm);
 
     return tfineToCoarse;
@@ -208,9 +208,9 @@ bool Foam::procFacesGAMGProcAgglomeration::doProcessorAgglomeration
     const lduMesh& mesh
 ) const
 {
-    // Check the need for further agglomeration on all processors
+    // Check the need for further agglomeration on any processors
     bool doAgg = mesh.lduAddr().size() < nAgglomeratingCells_;
-    mesh.reduce(doAgg, orOp<bool>());
+    UPstream::reduceOr(doAgg, mesh.comm());
     return doAgg;
 }
 

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C

@@ -1934,7 +1934,7 @@ Foam::pointField Foam::globalMeshData::geometricSharedPoints() const
     pointField sharedPoints(mesh_.points(), sharedPointLabels());
 
     // Append from all processors, globally consistent
-    Pstream::combineAllGather(sharedPoints, ListOps::appendEqOp<point>());
+    Pstream::combineReduce(sharedPoints, ListOps::appendEqOp<point>());
 
     // Merge tolerance
     scalar tolDim = matchTol_ * mesh_.bounds().mag();
@@ -2730,7 +2730,7 @@ void Foam::globalMeshData::updateMesh()
     (
         mesh_.nFaces(),
         sumOp<label>(),
-        Pstream::msgType(),
+        UPstream::msgType(),
         comm
     );
 
@@ -2743,7 +2743,7 @@ void Foam::globalMeshData::updateMesh()
     (
         mesh_.nCells(),
         sumOp<label>(),
-        Pstream::msgType(),
+        UPstream::msgType(),
         comm
     );
 
@@ -2756,7 +2756,7 @@ void Foam::globalMeshData::updateMesh()
     (
         mesh_.nPoints(),
         sumOp<label>(),
-        Pstream::msgType(),
+        UPstream::msgType(),
         comm
     );
 

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalPoints.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -276,7 +276,7 @@ bool Foam::globalPoints::mergeInfo
     // Updates database of current information on meshpoints with nbrInfo.  Uses
     // mergeInfo above. Returns true if data kept for meshPointi changed.
 
-    label infoChanged = false;
+    bool infoChanged = false;
 
     // Get the index into the procPoints list.
     const auto iter = meshToProcPoint_.cfind(localPointi);
@@ -325,7 +325,7 @@ bool Foam::globalPoints::storeInitialInfo
     // Updates database of current information on meshpoints with nbrInfo.  Uses
     // mergeInfo above. Returns true if data kept for meshPointi changed.
 
-    label infoChanged = false;
+    bool infoChanged = false;
 
     // Get the index into the procPoints list.
     const auto iter = meshToProcPoint_.find(localPointi);
@@ -929,8 +929,7 @@ void Foam::globalPoints::calculateSharedPoints
         neighbourList = meshToProcPoint_;
     }
 
-    // Exchange until nothing changes on all processors.
-    bool changed = false;
+    // Exchange until nothing changes on any processors.
 
     do
     {
@@ -959,10 +958,7 @@ void Foam::globalPoints::calculateSharedPoints
             changedPoints
         );
 
-        changed = changedPoints.size() > 0;
-        reduce(changed, orOp<bool>());
-
-    } while (changed);
+    } while (returnReduceOr(changedPoints.size()));
 
 
     //Pout<< "**ALL** connected points:" << endl;

src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeBase.H

@@ -412,7 +412,7 @@ public:
             const globalIndex&,
             labelList& elements,
             List<Map<label>>& compactMap,
-            const int tag = Pstream::msgType(),
+            const int tag = UPstream::msgType(),
             const label comm = UPstream::worldComm
         );
 
@@ -426,7 +426,7 @@ public:
             const globalIndex&,
             labelListList& cellCells,
             List<Map<label>>& compactMap,
-            const int tag = Pstream::msgType(),
+            const int tag = UPstream::msgType(),
             const label comm = UPstream::worldComm
         );
 

src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C

@@ -1046,9 +1046,7 @@ bool Foam::polyBoundaryMesh::checkParallelSync(const bool report) const
     }
 
     // Reduce (not broadcast) to respect local out-of-order errors (first loop)
-    reduce(hasError, orOp<bool>());
-
-    return hasError;
+    return returnReduceOr(hasError);
 }
 
 
@@ -1091,7 +1089,7 @@ bool Foam::polyBoundaryMesh::checkDefinition(const bool report) const
         nextPatchStart += bm[patchi].size();
     }
 
-    reduce(hasError, orOp<bool>());
+    Pstream::reduceOr(hasError);
 
     if (debug || report)
     {

src/OpenFOAM/meshes/polyMesh/polyMesh.C

@@ -62,8 +62,8 @@ void Foam::polyMesh::calcDirections() const
     // Knock out empty and wedge directions. Note:they will be present on all
     // domains.
 
-    label nEmptyPatches = 0;
-    label nWedgePatches = 0;
+    bool hasEmptyPatches = false;
+    bool hasWedgePatches = false;
 
     vector emptyDirVec = Zero;
     vector wedgeDirVec = Zero;
@@ -79,7 +79,7 @@ void Foam::polyMesh::calcDirections() const
 
             if (pp.size())
             {
-                nEmptyPatches++;
+                hasEmptyPatches = true;
                 emptyDirVec += sum(cmptMag(fa));
             }
         }
@@ -93,16 +93,14 @@ void Foam::polyMesh::calcDirections() const
 
             if (pp.size())
             {
-                nWedgePatches++;
+                hasWedgePatches = true;
                 wedgeDirVec += cmptMag(wpp.centreNormal());
             }
         }
     }
 
-    reduce(nEmptyPatches, maxOp<label>());
-    reduce(nWedgePatches, maxOp<label>());
 
-    if (nEmptyPatches)
+    if (returnReduceOr(hasEmptyPatches))
     {
         reduce(emptyDirVec, sumOp<vector>());
 
@@ -126,7 +124,7 @@ void Foam::polyMesh::calcDirections() const
 
     geometricD_ = solutionD_;
 
-    if (nWedgePatches)
+    if (returnReduceOr(hasWedgePatches))
     {
         reduce(wedgeDirVec, sumOp<vector>());
 
@@ -339,21 +337,21 @@ Foam::polyMesh::polyMesh(const IOobject& io, const bool doInit)
         neighbour_.write();
     }
 
-    // Warn if global empty mesh
-    if (returnReduce(boundary_.empty(), orOp<bool>()))
+    if (returnReduceOr(boundary_.empty()))
     {
         WarningInFunction
-            << "mesh missing boundary on one or more domains" << endl;
+            << "Missing mesh boundary on one or more domains" << endl;
 
-        if (returnReduce(nPoints(), sumOp<label>()) == 0)
+        // Warn if global empty mesh
+        if (returnReduceAnd(!nPoints()))
         {
             WarningInFunction
-                << "no points in mesh" << endl;
+                << "No points in mesh" << endl;
         }
-        if (returnReduce(nCells(), sumOp<label>()) == 0)
+        if (returnReduceAnd(!nCells()))
         {
             WarningInFunction
-                << "no cells in mesh" << endl;
+                << "No cells in mesh" << endl;
         }
     }
 
@@ -802,14 +800,10 @@ void Foam::polyMesh::resetPrimitives
         boundary_.calcGeometry();
 
         // Warn if global empty mesh
-        if
-        (
-            (returnReduce(nPoints(), sumOp<label>()) == 0)
-         || (returnReduce(nCells(), sumOp<label>()) == 0)
-        )
+        if (returnReduceAnd(!nPoints()) || returnReduceAnd(!nCells()))
         {
             FatalErrorInFunction
-                << "no points or no cells in mesh" << endl;
+                << "No points or no cells in mesh" << endl;
         }
     }
 }

src/OpenFOAM/meshes/polyMesh/polyMeshFromShapeMesh.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016, 2020 OpenFOAM Foundation
-    Copyright (C) 2018-2021 OpenCFD Ltd.
+    Copyright (C) 2018-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -634,7 +634,6 @@ Foam::polyMesh::polyMesh
 
     label nAllPatches = boundaryFaces.size();
 
-
     label nDefaultFaces = nFaces - defaultPatchStart;
     if (syncPar)
     {

src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C

@@ -776,7 +776,7 @@ Foam::labelList Foam::polyMeshTetDecomposition::adjustTetBasePtIs
 
     syncTools::syncFaceList(mesh, tetBasePtIs, maxEqOp<label>());
 
-    if (report && returnReduce(nAdapted, sumOp<label>()))
+    if (report && returnReduceOr(nAdapted))
     {
         Pout<< "Adapted starting point of triangulation on "
             << nAdapted << " faces." << endl;

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/symmetryPlane/symmetryPlanePolyPatch.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2015 OpenFOAM Foundation
+    Copyright (C) 2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -46,9 +47,9 @@ void Foam::symmetryPlanePolyPatch::calcGeometry(PstreamBuffers&)
 {
     if (n_ == vector::rootMax)
     {
-        if (returnReduce(size(), sumOp<label>()))
+        if (returnReduceOr(size()))
         {
-            const vectorField& nf(faceNormals());
+            const vectorField& nf = faceNormals();
             n_ = gAverage(nf);
 
             if (debug)

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -50,9 +51,9 @@ void Foam::wedgePolyPatch::calcGeometry(PstreamBuffers&)
         return;
     }
 
-    if (returnReduce(size(), sumOp<label>()))
+    if (returnReduceOr(size()))
     {
-        const vectorField& nf(faceNormals());
+        const vectorField& nf = faceNormals();
         n_ = gAverage(nf);
 
         if (debug)

src/OpenFOAM/meshes/polyMesh/syncTools/syncToolsTemplates.C

@@ -322,7 +322,7 @@ void Foam::syncTools::syncPointMap
             }
 
             // Broadcast: send merged values to all
-            Pstream::scatter(sharedPointValues);
+            Pstream::broadcast(sharedPointValues);
         }
 
         // Merge sharedPointValues (keyed on sharedPointAddr) into
@@ -667,7 +667,7 @@ void Foam::syncTools::syncEdgeMap
         }
 
         // Broadcast: send merged values to all
-        Pstream::scatter(sharedEdgeValues);
+        Pstream::broadcast(sharedEdgeValues);
     }
 
 

src/OpenFOAM/meshes/polyMesh/zones/faceZone/faceZone.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2017-2021 OpenCFD Ltd.
+    Copyright (C) 2017-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -601,7 +601,7 @@ bool Foam::faceZone::checkParallelSync(const bool report) const
         }
     }
 
-    return returnReduce(hasError, orOp<bool>());
+    return returnReduceOr(hasError);
 }
 
 

src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchTools.H

@@ -221,9 +221,48 @@ public:
     //
     //  \param[in] mergeDist Geometric merge tolerance for Foam::mergePoints
     //  \param[in] pp The patch to merge
-    //  \param[out] mergedPoints
-    //  \param[out] mergedFaces
-    //  \param[out] pointMergeMap
+    //  \param[out] mergedPoints merged points (master only, empty elsewhere)
+    //  \param[out] mergedFaces merged faces (master only, empty elsewhere)
+    //  \param[out] pointAddr Points globalIndex gather addressing
+    //      (master only, empty elsewhere)
+    //  \param[out] faceAddr Faces globalIndex gather addressing
+    //      (master only, empty elsewhere)
+    //  \param[out] pointMergeMap An old-to-new mapping from original
+    //      point index to the index into merged points.
+    //  \param[in]  useLocal  gather/merge patch localFaces/localPoints
+    //      instead of faces/points
+    //
+    //  \note
+    //  - OpenFOAM-v2112 and earlier: geometric merge on all patch points.
+    //  - OpenFOAM-v2206 and later: geometric merge on patch boundary points.
+    template<class FaceList, class PointField>
+    static void gatherAndMerge
+    (
+        const scalar mergeDist,
+        const PrimitivePatch<FaceList, PointField>& pp,
+        Field
+        <
+            typename PrimitivePatch<FaceList, PointField>::point_type
+        >& mergedPoints,
+        List
+        <
+            typename PrimitivePatch<FaceList, PointField>::face_type
+        >& mergedFaces,
+        globalIndex& pointAddr,
+        globalIndex& faceAddr,
+        labelList& pointMergeMap = const_cast<labelList&>(labelList::null()),
+        const bool useLocal = false
+    );
+
+    //- Gather points and faces onto master and merge into single patch.
+    //  Note: Normally uses faces/points (not localFaces/localPoints)
+    //
+    //  \param[in] mergeDist Geometric merge tolerance for Foam::mergePoints
+    //  \param[in] pp The patch to merge
+    //  \param[out] mergedPoints merged points (master only, empty elsewhere)
+    //  \param[out] mergedFaces merged faces (master only, empty elsewhere)
+    //  \param[out] pointMergeMap An old-to-new mapping from original
+    //      point index to the index into merged points.
     //  \param[in] useLocal gather/merge patch localFaces/localPoints
     //      instead of faces/points
     //
@@ -243,19 +282,20 @@ public:
         <
             typename PrimitivePatch<FaceList, PointField>::face_type
         >& mergedFaces,
-        labelList& pointMergeMap,
+        labelList& pointMergeMap = const_cast<labelList&>(labelList::null()),
         const bool useLocal = false
     );
 
     //- Gather (mesh!) points and faces onto master and merge collocated
     //  points into a single patch. Uses coupled point mesh
     //  structure so does not need tolerances.
-    //  On master and slave returns:
+    //  On master and sub-ranks returns:
     //  - pointToGlobal : for every local point index the global point index
     //  - uniqueMeshPointLabels : my local mesh points
     //  - globalPoints : global numbering for the global points
     //  - globalFaces : global numbering for the faces
-    //  On master only:
+    //  .
+    //  On master only returns:
     //  - mergedFaces : the merged faces
     //  - mergedPoints : the merged points
     template<class FaceList>

src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsGatherAndMerge.C

@@ -46,17 +46,19 @@ void Foam::PatchTools::gatherAndMerge
     <
         typename PrimitivePatch<FaceList, PointField>::face_type
     >& mergedFaces,
+    globalIndex& pointAddr,
+    globalIndex& faceAddr,
     labelList& pointMergeMap,
     const bool useLocal
 )
 {
-    typedef typename PrimitivePatch<FaceList,PointField>::face_type FaceType;
+    typedef typename PrimitivePatch<FaceList, PointField>::face_type FaceType;
 
     // Faces from all ranks
-    const globalIndex faceAddr(pp.size(), globalIndex::gatherOnly{});
+    faceAddr = globalIndex(pp.size(), globalIndex::gatherOnly{});
 
     // Points from all ranks
-    const globalIndex pointAddr
+    pointAddr = globalIndex
     (
         (useLocal ? pp.localPoints().size() : pp.points().size()),
         globalIndex::gatherOnly{}
@@ -152,6 +154,40 @@ void Foam::PatchTools::gatherAndMerge
 }
 
 
+template<class FaceList, class PointField>
+void Foam::PatchTools::gatherAndMerge
+(
+    const scalar mergeDist,
+    const PrimitivePatch<FaceList, PointField>& pp,
+    Field
+    <
+        typename PrimitivePatch<FaceList, PointField>::point_type
+    >& mergedPoints,
+    List
+    <
+        typename PrimitivePatch<FaceList, PointField>::face_type
+    >& mergedFaces,
+    labelList& pointMergeMap,
+    const bool useLocal
+)
+{
+    globalIndex pointAddr;
+    globalIndex faceAddr;
+
+    PatchTools::gatherAndMerge<FaceList, PointField>
+    (
+        mergeDist,
+        pp,
+        mergedPoints,
+        mergedFaces,
+        pointAddr,
+        faceAddr,
+        pointMergeMap,
+        useLocal
+    );
+}
+
+
 template<class FaceList>
 void Foam::PatchTools::gatherAndMerge
 (

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -425,8 +425,7 @@ bool Foam::primitiveMesh::checkFaceOrthogonality
 
     if (debug || report)
     {
-        label neiSize = ortho.size();
-        reduce(neiSize, sumOp<label>());
+        const label neiSize = returnReduce(ortho.size(), sumOp<label>());
 
         if (neiSize > 0)
         {
@@ -929,9 +928,9 @@ bool Foam::primitiveMesh::checkUpperTriangular
     label internal = nInternalFaces();
 
     // Has error occurred?
-    bool error = false;
-    // Have multiple faces been detected?
-    label nMultipleCells = false;
+    bool hasError = false;
+    // Multiple faces detected?
+    label nMultipleCells = 0;
 
     // Loop through faceCells once more and make sure that for internal cell
     // the first label is smaller
@@ -939,7 +938,7 @@ bool Foam::primitiveMesh::checkUpperTriangular
     {
         if (own[facei] >= nei[facei])
         {
-            error  = true;
+            hasError = true;
 
             if (setPtr)
             {
@@ -1024,7 +1023,7 @@ bool Foam::primitiveMesh::checkUpperTriangular
             }
             else if (thisFace < prevFace)
             {
-                error = true;
+                hasError = true;
 
                 if (setPtr)
                 {
@@ -1042,7 +1041,7 @@ bool Foam::primitiveMesh::checkUpperTriangular
         }
     }
 
-    reduce(error, orOp<bool>());
+    Pstream::reduceOr(hasError);
     reduce(nMultipleCells, sumOp<label>());
 
     if ((debug || report) && nMultipleCells > 0)
@@ -1051,7 +1050,7 @@ bool Foam::primitiveMesh::checkUpperTriangular
             << " neighbouring cells with multiple inbetween faces." << endl;
     }
 
-    if (error)
+    if (hasError)
     {
         if (debug || report)
         {

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheckEdgeLength.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -69,8 +69,7 @@ bool Foam::primitiveMesh::checkEdgeLength
     reduce(minLenSqr, minOp<scalar>());
     reduce(maxLenSqr, maxOp<scalar>());
 
-    label nSmall = smallEdgeSet.size();
-    reduce(nSmall, sumOp<label>());
+    label nSmall = returnReduce(smallEdgeSet.size(), sumOp<label>());
 
     if (setPtr)
     {

src/OpenFOAM/parallel/globalIndex/globalIndex.H

@@ -108,8 +108,8 @@ public:
 
     // Constructors
 
-        //- Default construct
-        globalIndex() = default;
+        //- Default construct (empty)
+        globalIndex() noexcept = default;
 
         //- Copy construct from a list of offsets.
         //- No communication required
@@ -184,9 +184,15 @@ public:
         //- Global max of localSizes
         inline label maxSize() const;
 
+
+    // Access
+
         //- Const-access to the offsets
         inline const labelList& offsets() const noexcept;
 
+        //- Write-access to the offsets, for changing after construction
+        inline labelList& offsets() noexcept;
+
 
     // Dimensions
 
@@ -202,8 +208,8 @@ public:
 
     // Edit
 
-        //- Write-access to the offsets, for changing after construction
-        inline labelList& offsets() noexcept;
+        //- Reset to be empty (no offsets)
+        inline void clear();
 
         //- Reset from local size, using gather/broadcast
         //- with default/specified communicator if parallel.

src/OpenFOAM/parallel/globalIndex/globalIndexI.H

@@ -176,6 +176,12 @@ inline Foam::labelList& Foam::globalIndex::offsets() noexcept
 }
 
 
+inline void Foam::globalIndex::clear()
+{
+    offsets_.clear();
+}
+
+
 inline const Foam::labelUList Foam::globalIndex::localStarts() const
 {
     const label len = (offsets_.size() - 1);

src/OpenFOAM/primitives/Scalar/scalar/scalar.C

@@ -40,6 +40,24 @@ Foam::scalar Foam::readScalar(Istream& is)
 }
 
 
+Foam::scalar Foam::readScalarOrDefault(Istream& is, const scalar defaultValue)
+{
+    if (is.good())
+    {
+        token tok(is);
+
+        if (tok.isNumber())
+        {
+            return tok.scalarToken();
+        }
+
+        is.putBack(tok);
+    }
+
+    return defaultValue;
+}
+
+
 Foam::scalar Foam::readRawScalar(Istream& is)
 {
     scalar val(0);

src/OpenFOAM/primitives/Scalar/scalar/scalar.H

@@ -103,6 +103,9 @@ namespace Foam
     //- Read scalar from stream.
     scalar readScalar(Istream& is);
 
+    //- Read scalar from stream if present or return default value
+    scalar readScalarOrDefault(Istream& is, const scalar defaultValue);
+
     //- Read raw scalar from binary stream.
     //  \note No internal check for binary vs ascii,
     //        the caller knows what they are doing
@@ -172,6 +175,9 @@ namespace Foam
     //- Read scalar from stream.
     scalar readScalar(Istream& is);
 
+    //- Read scalar from stream if present or return default value
+    scalar readScalarOrDefault(Istream& is, const scalar defaultValue);
+
     //- Read raw scalar from binary stream.
     //  \note No internal check for binary vs ascii,
     //        the caller knows what they are doing

src/OpenFOAM/primitives/globalIndexAndTransform/globalIndexAndTransform.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -269,7 +269,7 @@ void Foam::globalIndexAndTransform::determineTransforms()
     }
 
     transforms_.transfer(localTransforms);
-    Pstream::scatter(transforms_);
+    Pstream::broadcast(transforms_);
 }
 
 

src/OpenFOAM/primitives/random/Random/Random.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.
@@ -143,14 +143,14 @@ Foam::label Foam::Random::position(const label& start, const label& end)
 template<>
 Foam::scalar Foam::Random::globalSample01()
 {
-    scalar value(-GREAT);
+    scalar value(0);
 
     if (Pstream::master())
     {
         value = scalar01();
     }
 
-    Pstream::scatter(value);
+    Pstream::broadcast(value);
 
     return value;
 }
@@ -159,14 +159,14 @@ Foam::scalar Foam::Random::globalSample01()
 template<>
 Foam::label Foam::Random::globalSample01()
 {
-    label value(labelMin);
+    label value(0);
 
     if (Pstream::master())
     {
         value = round(scalar01());
     }
 
-    Pstream::scatter(value);
+    Pstream::broadcast(value);
 
     return value;
 }
@@ -175,14 +175,14 @@ Foam::label Foam::Random::globalSample01()
 template<>
 Foam::scalar Foam::Random::globalGaussNormal()
 {
-    scalar value(-GREAT);
+    scalar value(0);
 
     if (Pstream::master())
     {
         value = GaussNormal<scalar>();
     }
 
-    Pstream::scatter(value);
+    Pstream::broadcast(value);
 
     return value;
 }
@@ -191,14 +191,14 @@ Foam::scalar Foam::Random::globalGaussNormal()
 template<>
 Foam::label Foam::Random::globalGaussNormal()
 {
-    label value(labelMin);
+    label value(0);
 
     if (Pstream::master())
     {
         value = GaussNormal<label>();
     }
 
-    Pstream::scatter(value);
+    Pstream::broadcast(value);
 
     return value;
 }
@@ -211,14 +211,14 @@ Foam::scalar Foam::Random::globalPosition
     const scalar& end
 )
 {
-    scalar value(-GREAT);
+    scalar value(0);
 
     if (Pstream::master())
     {
         value = position<scalar>(start, end);
     }
 
-    Pstream::scatter(value);
+    Pstream::broadcast(value);
 
     return value;
 }
@@ -231,14 +231,14 @@ Foam::label Foam::Random::globalPosition
     const label& end
 )
 {
-    label value(labelMin);
+    label value(0);
 
     if (Pstream::master())
     {
         value = position<label>(start, end);
     }
 
-    Pstream::scatter(value);
+    Pstream::broadcast(value);
 
     return value;
 }

src/OpenFOAM/primitives/random/Random/RandomTemplates.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -92,14 +92,14 @@ void Foam::Random::shuffle(UList<Type>& values)
 template<class Type>
 Type Foam::Random::globalSample01()
 {
-    Type value = -GREAT*pTraits<Type>::one;
+    Type value(Zero);
 
     if (Pstream::master())
     {
         value = sample01<Type>();
     }
 
-    Pstream::scatter(value);
+    Pstream::broadcast(value);
 
     return value;
 }
@@ -108,14 +108,14 @@ Type Foam::Random::globalSample01()
 template<class Type>
 Type Foam::Random::globalGaussNormal()
 {
-    Type value = -GREAT*pTraits<Type>::one;
+    Type value(Zero);
 
     if (Pstream::master())
     {
         value = GaussNormal<Type>();
     }
 
-    Pstream::scatter(value);
+    Pstream::broadcast(value);
 
     return value;
 }
@@ -124,14 +124,14 @@ Type Foam::Random::globalGaussNormal()
 template<class Type>
 Type Foam::Random::globalPosition(const Type& start, const Type& end)
 {
-    Type value = -GREAT*pTraits<Type>::one;
+    Type value(Zero);
 
     if (Pstream::master())
     {
         value = position<Type>(start, end);
     }
 
-    Pstream::scatter(value);
+    Pstream::broadcast(value);
 
     return value;
 }
@@ -140,14 +140,12 @@ Type Foam::Random::globalPosition(const Type& start, const Type& end)
 template<class Type>
 void Foam::Random::globalRandomise01(Type& value)
 {
-    value = -GREAT*pTraits<Type>::one;
-
     if (Pstream::master())
     {
         value = sample01<Type>();
     }
 
-    Pstream::scatter(value);
+    Pstream::broadcast(value);
 }
 
 

src/Pstream/dummy/UOPstreamWrite.C

@@ -41,7 +41,7 @@ bool Foam::UOPstream::bufferIPCsend()
 
 bool Foam::UOPstream::write
 (
-    const commsTypes commsType,
+    const UPstream::commsTypes commsType,
     const int toProcNo,
     const char* buf,
     const std::streamsize bufSize,

src/Pstream/mpi/UOPstreamWrite.C

@@ -52,7 +52,7 @@ bool Foam::UOPstream::bufferIPCsend()
 
 bool Foam::UOPstream::write
 (
-    const commsTypes commsType,
+    const UPstream::commsTypes commsType,
     const int toProcNo,
     const char* buf,
     const std::streamsize bufSize,
@@ -65,7 +65,7 @@ bool Foam::UOPstream::write
         Pout<< "UOPstream::write : starting write to:" << toProcNo
             << " tag:" << tag
             << " comm:" << communicator << " size:" << label(bufSize)
-            << " commsType:" << UPstream::commsTypeNames[commsType]
+            << " commType:" << UPstream::commsTypeNames[commsType]
             << Foam::endl;
     }
     if (UPstream::warnComm != -1 && communicator != UPstream::warnComm)
@@ -73,7 +73,7 @@ bool Foam::UOPstream::write
         Pout<< "UOPstream::write : starting write to:" << toProcNo
             << " tag:" << tag
             << " comm:" << communicator << " size:" << label(bufSize)
-            << " commsType:" << UPstream::commsTypeNames[commsType]
+            << " commType:" << UPstream::commsTypeNames[commsType]
             << " warnComm:" << UPstream::warnComm
             << Foam::endl;
         error::printStack(Pout);
@@ -154,7 +154,7 @@ bool Foam::UOPstream::write
         {
             Pout<< "UOPstream::write : started write to:" << toProcNo
                 << " tag:" << tag << " size:" << label(bufSize)
-                << " commsType:" << UPstream::commsTypeNames[commsType]
+                << " commType:" << UPstream::commsTypeNames[commsType]
                 << " request:" << PstreamGlobals::outstandingRequests_.size()
                 << Foam::endl;
         }

src/TurbulenceModels/incompressible/turbulentTransportModels/turbulentTransportModels.C

@@ -126,6 +126,9 @@ makeLESModel(dynamicKEqn);
 #include "dynamicLagrangian.H"
 makeLESModel(dynamicLagrangian);
 
+#include "sigma.H"
+makeLESModel(sigma);
+
 #include "SpalartAllmarasDES.H"
 makeLESModel(SpalartAllmarasDES);
 

src/TurbulenceModels/schemes/DEShybrid/DEShybrid.H

@@ -7,6 +7,7 @@
 -------------------------------------------------------------------------------
     Copyright (C) 2015 OpenFOAM Foundation
     Copyright (C) 2015-2022 OpenCFD Ltd.
+    Copyright (C) 2022 Upstream CFD GmbH
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,7 +29,8 @@ Class
     Foam::DEShybrid
 
 Description
-    Hybrid convection scheme of Travin et al. for hybrid RAS/LES calculations.
+    Improved hybrid convection scheme of Travin et al. for hybrid RAS/LES
+    calculations with enhanced Grey Area Mitigation (GAM) behaviour.
 
     The scheme provides a blend between two convection schemes, based on local
     properties including the wall distance, velocity gradient and eddy
@@ -53,25 +55,26 @@ Description
 
     First published in:
     \verbatim
-        A. Travin, M. Shur, M. Strelets, P. Spalart (2000).
-        Physical and numerical upgrades in the detached-eddy simulation of
-        complex turbulent flows.
-        In Proceedings of the 412th Euromech Colloquium on LES and Complex
-        Transition and Turbulent Flows, Munich, Germany
+        Travin, A., Shur, M., Strelets, M., & Spalart, P. R. (2000).
+        Physical and numerical upgrades in the detached-eddy
+        simulation of complex turbulent flows.
+        In LES of Complex Transitional and Turbulent Flows.
+        Proceedings of the Euromech Colloquium 412. Munich, Germany
     \endverbatim
 
-    Original publication contained a typo for C_H3 constant. Corrected version
-    with minor changes for 2 lower limiters published in:
+    Original publication contained a typo for \c C_H3 constant.
+    Corrected version with minor changes for 2 lower limiters published in:
     \verbatim
-        P. Spalart, M. Shur, M. Strelets, A. Travin (2012).
-        Sensitivity of Landing-Gear Noise Predictions by Large-Eddy
-        Simulation to Numerics and Resolution.
-        AIAA Paper 2012-1174, 50th AIAA Aerospace Sciences Meeting,
-        Nashville / TN, Jan. 2012
+        Spalart, P., Shur, M., Strelets, M., & Travin, A. (2012).
+        Sensitivity of landing-gear noise predictions by large-eddy
+        simulation to numerics and resolution.
+        In 50th AIAA Aerospace Sciences Meeting Including the
+        New Horizons Forum and Aerospace Exposition. Nashville, US.
+        DOI:10.2514/6.2012-1174
     \endverbatim
 
-    Example of the DEShybrid scheme specification using linear within the LES
-    region and linearUpwind within the RAS region:
+    Example of the \c DEShybrid scheme specification using \c linear
+    within the LES region and \c linearUpwind within the RAS region:
     \verbatim
     divSchemes
     {
@@ -80,13 +83,14 @@ Description
         div(phi,U)      Gauss DEShybrid
             linear                    // scheme 1
             linearUpwind grad(U)      // scheme 2
-            hmax                      // LES delta name, e.g. 'delta', 'hmax'
-            0.65                      // DES coefficient, typically = 0.65
+            delta                     // LES delta name, e.g. 'delta', 'hmax'
+            0.65                      // CDES coefficient
             30                        // Reference velocity scale
             2                         // Reference length scale
             0                         // Minimum sigma limit (0-1)
             1                         // Maximum sigma limit (0-1)
-            1.0e-03;                  // Limiter of B function, typically 1e-03
+            1.0e-03                   // Limiter of B function, typically 1e-03
+            1.0;                      // nut limiter (if > 1, GAM extension is active)
         .
         .
     }
@@ -97,12 +101,12 @@ Notes
     be used in the detached/vortex shedding regions.
   - Scheme 2 should be an upwind/deferred correction/TVD scheme which will
     be used in the free-stream/Euler/boundary layer regions.
-  - the scheme is compiled into a separate library, and not available to
+  - The scheme is compiled into a separate library, and not available to
     solvers by default.  In order to use the scheme, add the library as a
     run-time loaded library in the \$FOAM\_CASE/system/controlDict
     dictionary, e.g.:
     \verbatim
-        libs ("libturbulenceModelSchemes.so");
+        libs (turbulenceModelSchemes);
     \endverbatim
 
 SourceFiles
@@ -110,15 +114,14 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef DEShybrid_H
-#define DEShybrid_H
+#ifndef Foam_DEShybrid_H
+#define Foam_DEShybrid_H
 
 #include "surfaceInterpolationScheme.H"
 #include "surfaceInterpolate.H"
 #include "fvcGrad.H"
 #include "blendedSchemeBase.H"
-#include "turbulentTransportModel.H"
-#include "turbulentFluidThermoModel.H"
+#include "turbulenceModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -135,6 +138,9 @@ class DEShybrid
     public surfaceInterpolationScheme<Type>,
     public blendedSchemeBase<Type>
 {
+    typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
+    typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
+
     // Private Data
 
         //- Scheme 1
@@ -146,7 +152,7 @@ class DEShybrid
         //- Name of the LES delta field
         word deltaName_;
 
-        //- DES Coefficient
+        //- DES coefficient
         scalar CDES_;
 
         //- Reference velocity scale [m/s]
@@ -164,10 +170,14 @@ class DEShybrid
         //- Limiter of B function
         scalar OmegaLim_;
 
+        //- Limiter for modified GAM behaviour
+        scalar nutLim_;
+
         //- Scheme constants
         scalar CH1_;
         scalar CH2_;
         scalar CH3_;
+        scalar Cs_;
 
         //- No copy construct
         DEShybrid(const DEShybrid&) = delete;
@@ -178,51 +188,108 @@ class DEShybrid
 
     // Private Member Functions
 
+        //- Check the scheme coefficients
+        void checkValues()
+        {
+            if (U0_.value() <= 0)
+            {
+                FatalErrorInFunction
+                    << "U0 coefficient must be > 0. "
+                    << "Current value: " << U0_ << exit(FatalError);
+            }
+            if (L0_.value() <= 0)
+            {
+                FatalErrorInFunction
+                    << "L0 coefficient must be > 0. "
+                    << "Current value: " << L0_ << exit(FatalError);
+            }
+            if (sigmaMin_ < 0)
+            {
+                FatalErrorInFunction
+                    << "sigmaMin coefficient must be >= 0. "
+                    << "Current value: " << sigmaMin_ << exit(FatalError);
+            }
+            if (sigmaMax_ < 0)
+            {
+                FatalErrorInFunction
+                    << "sigmaMax coefficient must be >= 0. "
+                    << "Current value: " << sigmaMax_ << exit(FatalError);
+            }
+            if (sigmaMin_ > 1)
+            {
+                FatalErrorInFunction
+                    << "sigmaMin coefficient must be <= 1. "
+                    << "Current value: " << sigmaMin_ << exit(FatalError);
+            }
+            if (sigmaMax_ > 1)
+            {
+                FatalErrorInFunction
+                    << "sigmaMax coefficient must be <= 1. "
+                    << "Current value: " << sigmaMax_ << exit(FatalError);
+            }
+
+            if (debug)
+            {
+                Info<< type() << "coefficients:" << nl
+                    << "    delta : " << deltaName_ << nl
+                    << "    CDES : " << CDES_ << nl
+                    << "    U0 : " << U0_.value() << nl
+                    << "    L0 : " << L0_.value() << nl
+                    << "    sigmaMin : " << sigmaMin_ << nl
+                    << "    sigmaMax : " << sigmaMax_ << nl
+                    << "    OmegaLim : " << OmegaLim_ << nl
+                    << "    nutLim : " << nutLim_ << nl
+                    << "    CH1 : " << CH1_ << nl
+                    << "    CH2 : " << CH2_ << nl
+                    << "    CH3 : " << CH3_ << nl
+                    << "    Cs : " << Cs_ << nl
+                    << endl;
+            }
+        }
+
+
         //- Calculate the blending factor
         tmp<surfaceScalarField> calcBlendingFactor
         (
-            const GeometricField<Type, fvPatchField, volMesh>& vf,
-            const volScalarField& nuEff,
+            const VolFieldType& vf,
+            const volScalarField& nut,
+            const volScalarField& nu,
             const volVectorField& U,
             const volScalarField& delta
         ) const
         {
-            tmp<volTensorField> gradU(fvc::grad(U));
-            const volScalarField S(sqrt(2.0)*mag(symm(gradU())));
-            const volScalarField Omega(sqrt(2.0)*mag(skew(gradU())));
+            tmp<volTensorField> tgradU = fvc::grad(U);
+            const volTensorField& gradU = tgradU.cref();
+            const volScalarField S(sqrt(2.0)*mag(symm(gradU)));
+            const volScalarField Omega(sqrt(2.0)*mag(skew(tgradU)));
             const dimensionedScalar tau0_ = L0_/U0_;
 
-            const volScalarField B
-            (
+            tmp<volScalarField> tB =
                 CH3_*Omega*max(S, Omega)
-               /max(0.5*(sqr(S) + sqr(Omega)), sqr(OmegaLim_/tau0_))
-            );
+               /max(0.5*(sqr(S) + sqr(Omega)), sqr(OmegaLim_/tau0_));
 
-            const volScalarField K
-            (
-                max(Foam::sqrt(0.5*(sqr(S) + sqr(Omega))), 0.1/tau0_)
-            );
+            tmp<volScalarField> tg = tanh(pow4(tB));
 
-            const volScalarField lTurb
-            (
+            tmp<volScalarField> tK =
+                max(Foam::sqrt(0.5*(sqr(S) + sqr(Omega))), 0.1/tau0_);
+
+            tmp<volScalarField> tlTurb =
                 Foam::sqrt
                 (
                     max
                     (
-                        nuEff/(pow(0.09, 1.5)*K),
-                        dimensionedScalar("l0", sqr(dimLength), 0)
+                        (max(nut, min(sqr(Cs_*delta)*S, nutLim_*nut)) + nu)
+                       /(pow(0.09, 1.5)*tK),
+                        dimensionedScalar(sqr(dimLength), Zero)
                     )
-                )
-            );
-
-            const volScalarField g(tanh(pow4(B)));
+                );
 
             const volScalarField A
             (
                 CH2_*max
                 (
                     scalar(0),
-                    CDES_*delta/max(lTurb*g, SMALL*L0_) - 0.5
+                    CDES_*delta/max(tlTurb*tg, SMALL*L0_) - 0.5
                 )
             );
 
@@ -297,14 +364,8 @@ public:
         DEShybrid(const fvMesh& mesh, Istream& is)
         :
             surfaceInterpolationScheme<Type>(mesh),
-            tScheme1_
-            (
-                surfaceInterpolationScheme<Type>::New(mesh, is)
-            ),
-            tScheme2_
-            (
-                surfaceInterpolationScheme<Type>::New(mesh, is)
-            ),
+            tScheme1_(surfaceInterpolationScheme<Type>::New(mesh, is)),
+            tScheme2_(surfaceInterpolationScheme<Type>::New(mesh, is)),
             deltaName_(is),
             CDES_(readScalar(is)),
             U0_("U0", dimLength/dimTime, readScalar(is)),
@@ -312,46 +373,13 @@ public:
             sigmaMin_(readScalar(is)),
             sigmaMax_(readScalar(is)),
             OmegaLim_(readScalar(is)),
+            nutLim_(readScalarOrDefault(is, scalar(1))),
             CH1_(3.0),
             CH2_(1.0),
-            CH3_(2.0)
+            CH3_(2.0),
+            Cs_(0.18)
         {
-            if (U0_.value() <= 0)
-            {
-                FatalErrorInFunction
-                    << "U0 coefficient must be > 0. "
-                    << "Current value: " << U0_ << exit(FatalError);
-            }
-            if (L0_.value() <= 0)
-            {
-                FatalErrorInFunction
-                    << "L0 coefficient must be > 0. "
-                    << "Current value: " << L0_ << exit(FatalError);
-            }
-            if (sigmaMin_ < 0)
-            {
-                FatalErrorInFunction
-                    << "sigmaMin coefficient must be >= 0. "
-                    << "Current value: " << sigmaMin_ << exit(FatalError);
-            }
-            if (sigmaMax_ < 0)
-            {
-                FatalErrorInFunction
-                    << "sigmaMax coefficient must be >= 0. "
-                    << "Current value: " << sigmaMax_ << exit(FatalError);
-            }
-            if (sigmaMin_ > 1)
-            {
-                FatalErrorInFunction
-                    << "sigmaMin coefficient must be <= 1. "
-                    << "Current value: " << sigmaMin_ << exit(FatalError);
-            }
-            if (sigmaMax_ > 1)
-            {
-                FatalErrorInFunction
-                    << "sigmaMax coefficient must be <= 1. "
-                    << "Current value: " << sigmaMax_ << exit(FatalError);
-            }
+            checkValues();
         }
 
         //- Construct from mesh, faceFlux and Istream
@@ -378,46 +406,13 @@ public:
             sigmaMin_(readScalar(is)),
             sigmaMax_(readScalar(is)),
             OmegaLim_(readScalar(is)),
+            nutLim_(readScalarOrDefault(is, scalar(1))),
             CH1_(3.0),
             CH2_(1.0),
-            CH3_(2.0)
+            CH3_(2.0),
+            Cs_(0.18)
         {
-            if (U0_.value() <= 0)
-            {
-                FatalErrorInFunction
-                    << "U0 coefficient must be > 0. "
-                    << "Current value: " << U0_ << exit(FatalError);
-            }
-            if (L0_.value() <= 0)
-            {
-                FatalErrorInFunction
-                    << "L0 coefficient must be > 0. "
-                    << "Current value: " << U0_ << exit(FatalError);
-            }
-            if (sigmaMin_ < 0)
-            {
-                FatalErrorInFunction
-                    << "sigmaMin coefficient must be >= 0. "
-                    << "Current value: " << sigmaMin_ << exit(FatalError);
-            }
-            if (sigmaMax_ < 0)
-            {
-                FatalErrorInFunction
-                    << "sigmaMax coefficient must be >= 0. "
-                    << "Current value: " << sigmaMax_ << exit(FatalError);
-            }
-            if (sigmaMin_ > 1)
-            {
-                FatalErrorInFunction
-                    << "sigmaMin coefficient must be <= 1. "
-                    << "Current value: " << sigmaMin_ << exit(FatalError);
-            }
-            if (sigmaMax_ > 1)
-            {
-                FatalErrorInFunction
-                    << "sigmaMax coefficient must be <= 1. "
-                    << "Current value: " << sigmaMax_ << exit(FatalError);
-            }
+            checkValues();
         }
 
 
@@ -431,33 +426,24 @@ public:
         {
             const fvMesh& mesh = this->mesh();
 
-            typedef compressible::turbulenceModel cmpModel;
-            typedef incompressible::turbulenceModel icoModel;
+            // Retrieve LES delta from the mesh database
+            const auto& delta =
+                mesh.lookupObject<const volScalarField>(deltaName_);
 
-            // Lookup the LES delta from the mesh database
-            const volScalarField& delta = this->mesh().template
-                 lookupObject<const volScalarField>(deltaName_);
+            // Retrieve turbulence model from the mesh database
+            const auto* modelPtr =
+                mesh.cfindObject<turbulenceModel>
+                (
+                    turbulenceModel::propertiesName
+                );
 
-            // Could avoid the compressible/incompressible case by looking
-            // up all fields from the database - but retrieving from model
-            // ensures consistent fields are being employed e.g. for multiphase
-            // where group name is used
-
-            if (mesh.foundObject<icoModel>(icoModel::propertiesName))
+            if (modelPtr)
             {
-                const icoModel& model =
-                    mesh.lookupObject<icoModel>(icoModel::propertiesName);
-
-                return calcBlendingFactor(vf, model.nuEff(), model.U(), delta);
-            }
-            else if (mesh.foundObject<cmpModel>(cmpModel::propertiesName))
-            {
-                const cmpModel& model =
-                    mesh.lookupObject<cmpModel>(cmpModel::propertiesName);
+                const auto& model = *modelPtr;
 
                 return calcBlendingFactor
                 (
-                    vf, model.muEff()/model.rho(), model.U(), delta
+                    vf, model.nut(), model.nu(), model.U(), delta
                 );
             }
 
@@ -471,12 +457,9 @@ public:
 
 
         //- Return the interpolation weighting factors
-        tmp<surfaceScalarField> weights
-        (
-            const GeometricField<Type, fvPatchField, volMesh>& vf
-        ) const
+        tmp<surfaceScalarField> weights(const VolFieldType& vf) const
         {
-            surfaceScalarField bf(blendingFactor(vf));
+            const surfaceScalarField bf(blendingFactor(vf));
 
             return
                 (scalar(1) - bf)*tScheme1_().weights(vf)
@@ -485,14 +468,10 @@ public:
 
 
         //- Return the face-interpolate of the given cell field
-        //  with explicit correction
-        tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
-        interpolate
-        (
-            const GeometricField<Type, fvPatchField, volMesh>& vf
-        ) const
+        //- with explicit correction
+        tmp<SurfaceFieldType> interpolate(const VolFieldType& vf) const
         {
-            surfaceScalarField bf(blendingFactor(vf));
+            const surfaceScalarField bf(blendingFactor(vf));
 
             return
                 (scalar(1) - bf)*tScheme1_().interpolate(vf)
@@ -508,14 +487,10 @@ public:
 
 
         //- Return the explicit correction to the face-interpolate
-        //  for the given field
-        virtual tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
-        correction
-        (
-            const GeometricField<Type, fvPatchField, volMesh>& vf
-        ) const
+        //- for the given field
+        virtual tmp<SurfaceFieldType> correction(const VolFieldType& vf) const
         {
-            surfaceScalarField bf(blendingFactor(vf));
+            const surfaceScalarField bf(blendingFactor(vf));
 
             if (tScheme1_().corrected())
             {

src/TurbulenceModels/schemes/Make/options

@@ -1,19 +1,9 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/transportModels \
-    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude
 
 LIB_LIBS = \
     -lfiniteVolume \
     -lmeshTools \
-    -lcompressibleTransportModels \
-    -lturbulenceModels \
-    -lincompressibleTurbulenceModels \
-    -lincompressibleTransportModels \
-    -lcompressibleTurbulenceModels \
-    -lfluidThermophysicalModels
+    -lturbulenceModels

src/TurbulenceModels/turbulenceModels/Base/SpalartAllmaras/SpalartAllmarasBase.C

@@ -0,0 +1,499 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2016-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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "SpalartAllmarasBase.H"
+#include "wallDist.H"
+#include "bound.H"
+#include "fvOptions.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField> SpalartAllmarasBase<BasicEddyViscosityModel>::chi() const
+{
+    return nuTilda_/this->nu();
+}
+
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField> SpalartAllmarasBase<BasicEddyViscosityModel>::fv1
+(
+    const volScalarField& chi
+) const
+{
+    const volScalarField chi3("chi3", pow3(chi));
+    return chi3/(chi3 + pow3(Cv1_));
+}
+
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField> SpalartAllmarasBase<BasicEddyViscosityModel>::fv2
+(
+    const volScalarField& chi,
+    const volScalarField& fv1
+) const
+{
+    return scalar(1) - chi/(scalar(1) + chi*fv1);
+}
+
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField> SpalartAllmarasBase<BasicEddyViscosityModel>::ft2
+(
+    const volScalarField& chi
+) const
+{
+    if (ft2_)
+    {
+        return Ct3_*exp(-Ct4_*sqr(chi));
+    }
+
+    return tmp<volScalarField>::New
+    (
+        IOobject
+        (
+            "ft2",
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        this->mesh_,
+        dimensionedScalar(dimless, Zero)
+    );
+}
+
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField> SpalartAllmarasBase<BasicEddyViscosityModel>::Omega
+(
+    const volTensorField& gradU
+) const
+{
+    return sqrt(2.0)*mag(skew(gradU));
+}
+
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField> SpalartAllmarasBase<BasicEddyViscosityModel>::r
+(
+    const volScalarField& nur,
+    const volScalarField& Stilda,
+    const volScalarField& dTilda
+) const
+{
+    const dimensionedScalar eps(Stilda.dimensions(), SMALL);
+
+    tmp<volScalarField> tr =
+        min(nur/(max(Stilda, eps)*sqr(kappa_*dTilda)), scalar(10));
+
+    tr.ref().boundaryFieldRef() == 0;
+
+    return tr;
+}
+
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField::Internal> SpalartAllmarasBase<BasicEddyViscosityModel>::fw
+(
+    const volScalarField& Stilda,
+    const volScalarField& dTilda
+) const
+{
+    const volScalarField::Internal r(this->r(nuTilda_, Stilda, dTilda)()());
+    const volScalarField::Internal g(r + Cw2_*(pow6(r) - r));
+
+    return g*pow((1 + pow6(Cw3_))/(pow6(g) + pow6(Cw3_)), 1.0/6.0);
+}
+
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField> SpalartAllmarasBase<BasicEddyViscosityModel>::Stilda
+(
+    const volScalarField& chi,
+    const volScalarField& fv1,
+    const volTensorField& gradU,
+    const volScalarField& dTilda
+) const
+{
+    const volScalarField Omega(this->Omega(gradU));
+
+    return
+        max
+        (
+            Omega + fv2(chi, fv1)*nuTilda_/sqr(kappa_*dTilda),
+            Cs_*Omega
+        );
+}
+
+
+template<class BasicEddyViscosityModel>
+void SpalartAllmarasBase<BasicEddyViscosityModel>::correctNut
+(
+    const volScalarField& fv1
+)
+{
+    this->nut_ = nuTilda_*fv1;
+    this->nut_.correctBoundaryConditions();
+    fv::options::New(this->mesh_).correct(this->nut_);
+}
+
+
+template<class BasicEddyViscosityModel>
+void SpalartAllmarasBase<BasicEddyViscosityModel>::correctNut()
+{
+    correctNut(fv1(this->chi()));
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class BasicEddyViscosityModel>
+SpalartAllmarasBase<BasicEddyViscosityModel>::SpalartAllmarasBase
+(
+    const word& type,
+    const alphaField& alpha,
+    const rhoField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& alphaRhoPhi,
+    const surfaceScalarField& phi,
+    const transportModel& transport,
+    const word& propertiesName
+)
+:
+    BasicEddyViscosityModel
+    (
+        type,
+        alpha,
+        rho,
+        U,
+        alphaRhoPhi,
+        phi,
+        transport,
+        propertiesName
+    ),
+
+    sigmaNut_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "sigmaNut",
+            this->coeffDict_,
+            0.66666
+        )
+    ),
+    kappa_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "kappa",
+            this->coeffDict_,
+            0.41
+        )
+    ),
+    Cb1_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Cb1",
+            this->coeffDict_,
+            0.1355
+        )
+    ),
+    Cb2_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Cb2",
+            this->coeffDict_,
+            0.622
+        )
+    ),
+    Cw1_(Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_),
+    Cw2_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Cw2",
+            this->coeffDict_,
+            0.3
+        )
+    ),
+    Cw3_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Cw3",
+            this->coeffDict_,
+            2.0
+        )
+    ),
+    Cv1_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Cv1",
+            this->coeffDict_,
+            7.1
+        )
+    ),
+    Cs_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Cs",
+            this->coeffDict_,
+            0.3
+        )
+    ),
+    ck_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "ck",
+            this->coeffDict_,
+            0.07
+        )
+    ),
+    ft2_
+    (
+        Switch::getOrAddToDict
+        (
+            "ft2",
+            this->coeffDict_,
+            false
+        )
+    ),
+    Ct3_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Ct3",
+            this->coeffDict_,
+            1.2
+        )
+    ),
+    Ct4_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Ct4",
+            this->coeffDict_,
+            0.5
+        )
+    ),
+
+    nuTilda_
+    (
+        IOobject
+        (
+            "nuTilda",
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh_
+    ),
+
+    y_(wallDist::New(this->mesh_).y())
+{
+    if (ft2_)
+    {
+        Info<< "ft2 term: active" << nl;
+    }
+    else
+    {
+        Info<< "ft2 term: inactive" << nl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class BasicEddyViscosityModel>
+bool SpalartAllmarasBase<BasicEddyViscosityModel>::read()
+{
+    if (BasicEddyViscosityModel::read())
+    {
+        sigmaNut_.readIfPresent(this->coeffDict());
+        kappa_.readIfPresent(this->coeffDict());
+
+        Cb1_.readIfPresent(this->coeffDict());
+        Cb2_.readIfPresent(this->coeffDict());
+        Cw1_ = Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_;
+        Cw2_.readIfPresent(this->coeffDict());
+        Cw3_.readIfPresent(this->coeffDict());
+        Cv1_.readIfPresent(this->coeffDict());
+        Cs_.readIfPresent(this->coeffDict());
+        ck_.readIfPresent(this->coeffDict());
+
+        ft2_.readIfPresent("ft2", this->coeffDict());
+        Ct3_.readIfPresent(this->coeffDict());
+        Ct4_.readIfPresent(this->coeffDict());
+
+        if (ft2_)
+        {
+            Info<< "    ft2 term: active" << nl;
+        }
+        else
+        {
+            Info<< "    ft2 term: inactive" << nl;
+        }
+
+        return true;
+    }
+
+    return false;
+}
+
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField>
+SpalartAllmarasBase<BasicEddyViscosityModel>::DnuTildaEff() const
+{
+    return tmp<volScalarField>::New
+    (
+        IOobject::groupName("DnuTildaEff", this->alphaRhoPhi_.group()),
+        (nuTilda_ + this->nu())/sigmaNut_
+    );
+}
+
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField> SpalartAllmarasBase<BasicEddyViscosityModel>::k() const
+{
+    // (B:Eq. 4.50)
+    const scalar Cmu = 0.09;
+    const auto fv1 = this->fv1(chi());
+
+    return tmp<volScalarField>::New
+    (
+        IOobject::groupName("k", this->alphaRhoPhi_.group()),
+        cbrt(fv1)*nuTilda_*::sqrt(scalar(2)/Cmu)*mag(symm(fvc::grad(this->U_)))
+    );
+}
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField>
+SpalartAllmarasBase<BasicEddyViscosityModel>::epsilon() const
+{
+    // (B:Eq. 4.50)
+    const scalar Cmu = 0.09;
+    const auto fv1 = this->fv1(chi());
+    const dimensionedScalar nutSMALL(nuTilda_.dimensions(), SMALL);
+
+    return tmp<volScalarField>::New
+    (
+        IOobject::groupName("epsilon", this->alphaRhoPhi_.group()),
+        sqrt(fv1)*sqr(::sqrt(Cmu)*this->k())/(nuTilda_ + this->nut_ + nutSMALL)
+    );
+}
+
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField> SpalartAllmarasBase<BasicEddyViscosityModel>::omega() const
+{
+    // (P:p. 384)
+    const scalar betaStar = 0.09;
+    const dimensionedScalar k0(sqr(dimLength/dimTime), SMALL);
+
+    return tmp<volScalarField>::New
+    (
+        IOobject::groupName("omega", this->alphaRhoPhi_.group()),
+        this->epsilon()/(betaStar*(this->k() + k0))
+    );
+}
+
+
+template<class BasicEddyViscosityModel>
+void SpalartAllmarasBase<BasicEddyViscosityModel>::correct()
+{
+    if (!this->turbulence_)
+    {
+        return;
+    }
+
+    {
+        // Local references
+        const alphaField& alpha = this->alpha_;
+        const rhoField& rho = this->rho_;
+        const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
+        const volVectorField& U = this->U_;
+        fv::options& fvOptions(fv::options::New(this->mesh_));
+
+        BasicEddyViscosityModel::correct();
+
+        const volScalarField chi(this->chi());
+        const volScalarField fv1(this->fv1(chi));
+        const volScalarField ft2(this->ft2(chi));
+
+        tmp<volTensorField> tgradU = fvc::grad(U);
+        volScalarField dTilda(this->dTilda(chi, fv1, tgradU()));
+        volScalarField Stilda(this->Stilda(chi, fv1, tgradU(), dTilda));
+        tgradU.clear();
+
+        tmp<fvScalarMatrix> nuTildaEqn
+        (
+            fvm::ddt(alpha, rho, nuTilda_)
+          + fvm::div(alphaRhoPhi, nuTilda_)
+          - fvm::laplacian(alpha*rho*DnuTildaEff(), nuTilda_)
+          - Cb2_/sigmaNut_*alpha()*rho()*magSqr(fvc::grad(nuTilda_)()())
+         ==
+            Cb1_*alpha()*rho()*Stilda()*nuTilda_()*(scalar(1) - ft2())
+          - fvm::Sp
+            (
+                (Cw1_*fw(Stilda, dTilda) - Cb1_/sqr(kappa_)*ft2())
+               *alpha()*rho()*nuTilda_()/sqr(dTilda()),
+                nuTilda_
+            )
+          + fvOptions(alpha, rho, nuTilda_)
+        );
+
+        nuTildaEqn.ref().relax();
+        fvOptions.constrain(nuTildaEqn.ref());
+        solve(nuTildaEqn);
+        fvOptions.correct(nuTilda_);
+        bound(nuTilda_, dimensionedScalar(nuTilda_.dimensions(), Zero));
+        nuTilda_.correctBoundaryConditions();
+    }
+
+    correctNut();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/Base/SpalartAllmaras/SpalartAllmarasBase.H

@@ -0,0 +1,246 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2015-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/>.
+
+Class
+    Foam::LESModels::SpalartAllmarasBase
+
+Group
+    grpDESTurbulence
+
+Description
+    Base class to handle various characteristics for \c SpalartAllmaras based
+    LES/DES turbulence models for incompressible and compressible flows.
+
+    References:
+    \verbatim
+        Standard model:
+            Spalart, P.R., & Allmaras, S.R. (1994).
+            A one-equation turbulence model for aerodynamic flows.
+            La Recherche Aerospatiale, 1, 5-21.
+
+        Standard model:
+            Spalart, P. R., Jou, W. H., Strelets, M., & Allmaras, S. R. (1997).
+            Comments on the feasibility of LES for wings, and on a hybrid
+            RANS/LES approach.
+            Advances in DNS/LES, 1, 4-8.
+
+        Estimation expression for k and epsilon (tag:B), Eq. 4.50:
+            Bourgoin, A. (2019).
+            Bathymetry induced turbulence modelling the
+            Alderney Race site: regional approach with TELEMAC-LES.
+            Normandie Université.
+
+        Estimation expressions for omega (tag:P):
+            Pope, S. B. (2000).
+            Turbulent flows.
+            Cambridge, UK: Cambridge Univ. Press
+            DOI:10.1017/CBO9780511840531
+    \endverbatim
+
+SourceFiles
+    SpalartAllmarasBase.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_SpalartAllmarasBase_H
+#define Foam_SpalartAllmarasBase_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class SpalartAllmarasBase Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class BasicEddyViscosityModel>
+class SpalartAllmarasBase
+:
+    public BasicEddyViscosityModel
+{
+    // Private Member Functions
+
+        //- No copy construct
+        SpalartAllmarasBase(const SpalartAllmarasBase&) = delete;
+
+        //- No copy assignment
+        void operator=(const SpalartAllmarasBase&) = delete;
+
+
+protected:
+
+    // Protected Data
+
+        // Model constants
+
+            dimensionedScalar sigmaNut_;
+            dimensionedScalar kappa_;
+
+            dimensionedScalar Cb1_;
+            dimensionedScalar Cb2_;
+            dimensionedScalar Cw1_;
+            dimensionedScalar Cw2_;
+            dimensionedScalar Cw3_;
+            dimensionedScalar Cv1_;
+            dimensionedScalar Cs_;
+            dimensionedScalar ck_;
+
+            Switch ft2_;
+            dimensionedScalar Ct3_;
+            dimensionedScalar Ct4_;
+
+
+        // Fields
+
+            //- Modified kinematic viscosity [m^2/s]
+            volScalarField nuTilda_;
+
+            //- Wall distance
+            //  Note: different to wall distance in parent RASModel
+            //  which is for near-wall cells only
+            const volScalarField& y_;
+
+
+    // Protected Member Functions
+
+        tmp<volScalarField> chi() const;
+
+        tmp<volScalarField> fv1(const volScalarField& chi) const;
+
+        tmp<volScalarField> fv2
+        (
+            const volScalarField& chi,
+            const volScalarField& fv1
+        ) const;
+
+        tmp<volScalarField> ft2(const volScalarField& chi) const;
+
+        tmp<volScalarField> Omega(const volTensorField& gradU) const;
+
+        tmp<volScalarField> r
+        (
+            const volScalarField& nur,
+            const volScalarField& Stilda,
+            const volScalarField& dTilda
+        ) const;
+
+        tmp<volScalarField::Internal> fw
+        (
+            const volScalarField& Stilda,
+            const volScalarField& dTilda
+        ) const;
+
+        virtual tmp<volScalarField> Stilda
+        (
+            const volScalarField& chi,
+            const volScalarField& fv1,
+            const volTensorField& gradU,
+            const volScalarField& dTilda
+        ) const;
+
+        //- Length scale
+        virtual tmp<volScalarField> dTilda
+        (
+            const volScalarField& chi,
+            const volScalarField& fv1,
+            const volTensorField& gradU
+        ) const = 0;
+
+        void correctNut(const volScalarField& fv1);
+        virtual void correctNut();
+
+
+public:
+
+    typedef typename BasicEddyViscosityModel::alphaField alphaField;
+    typedef typename BasicEddyViscosityModel::rhoField rhoField;
+    typedef typename BasicEddyViscosityModel::transportModel transportModel;
+
+
+    // Constructors
+
+        //- Construct from components
+        SpalartAllmarasBase
+        (
+            const word& type,
+            const alphaField& alpha,
+            const rhoField& rho,
+            const volVectorField& U,
+            const surfaceScalarField& alphaRhoPhi,
+            const surfaceScalarField& phi,
+            const transportModel& transport,
+            const word& propertiesName = turbulenceModel::propertiesName
+        );
+
+
+    //- Destructor
+    virtual ~SpalartAllmarasBase() = default;
+
+
+    // Member Functions
+
+        //- Re-read model coefficients if they have changed
+        virtual bool read();
+
+        //- Return the effective diffusivity for nuTilda
+        tmp<volScalarField> DnuTildaEff() const;
+
+        //- Return the (estimated) turbulent kinetic energy
+        virtual tmp<volScalarField> k() const;
+
+        //- Return the (estimated) turbulent kinetic energy dissipation rate
+        virtual tmp<volScalarField> epsilon() const;
+
+        //- Return the (estimated) specific dissipation rate
+        virtual tmp<volScalarField> omega() const;
+
+        //- Return the modified kinematic viscosity
+        tmp<volScalarField> nuTilda() const
+        {
+            return nuTilda_;
+        }
+
+        //- Correct nuTilda and related properties
+        virtual void correct();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "SpalartAllmarasBase.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/Base/kOmegaSST/kOmegaSSTBase.C

@@ -6,7 +6,8 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2022 Upstream CFD GmbH
+    Copyright (C) 2016-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -47,7 +48,7 @@ tmp<volScalarField> kOmegaSSTBase<BasicEddyViscosityModel>::F1
     tmp<volScalarField> CDkOmegaPlus = max
     (
         CDkOmega,
-        dimensionedScalar("1.0e-10", dimless/sqr(dimTime), 1.0e-10)
+        dimensionedScalar(dimless/sqr(dimTime), 1.0e-10)
     );
 
     tmp<volScalarField> arg1 = min
@@ -132,6 +133,17 @@ void kOmegaSSTBase<BasicEddyViscosityModel>::correctNut()
 }
 
 
+template<class BasicEddyViscosityModel>
+Foam::tmp<Foam::volScalarField> kOmegaSSTBase<BasicEddyViscosityModel>::S2
+(
+    const volScalarField& F1,
+    const volTensorField& gradU
+) const
+{
+    return 2*magSqr(symm(gradU));
+}
+
+
 template<class BasicEddyViscosityModel>
 tmp<volScalarField::Internal> kOmegaSSTBase<BasicEddyViscosityModel>::Pk
 (
@@ -143,17 +155,32 @@ tmp<volScalarField::Internal> kOmegaSSTBase<BasicEddyViscosityModel>::Pk
 
 
 template<class BasicEddyViscosityModel>
-tmp<volScalarField::Internal>
-kOmegaSSTBase<BasicEddyViscosityModel>::epsilonByk
+tmp<volScalarField::Internal> kOmegaSSTBase<BasicEddyViscosityModel>::epsilonByk
 (
-    const volScalarField& F1,
-    const volTensorField& gradU
+    const volScalarField& /* F1 not used */,
+    const volTensorField& /* gradU not used */
 ) const
 {
     return betaStar_*omega_();
 }
 
 
+template<class BasicEddyViscosityModel>
+tmp<volScalarField::Internal> kOmegaSSTBase<BasicEddyViscosityModel>::GbyNu0
+(
+    const volTensorField& gradU,
+    const volScalarField& F1,
+    const volScalarField& S2
+) const
+{
+    return tmp<volScalarField::Internal>::New
+    (
+        IOobject::scopedName(this->type(), "GbyNu"),
+        gradU() && dev(twoSymm(gradU()))
+    );
+}
+
+
 template<class BasicEddyViscosityModel>
 tmp<volScalarField::Internal> kOmegaSSTBase<BasicEddyViscosityModel>::GbyNu
 (
@@ -165,8 +192,7 @@ tmp<volScalarField::Internal> kOmegaSSTBase<BasicEddyViscosityModel>::GbyNu
     return min
     (
         GbyNu0,
-        (c1_/a1_)*betaStar_*omega_()
-       *max(a1_*omega_(), b1_*F2*sqrt(S2))
+        (c1_/a1_)*betaStar_*omega_()*max(a1_*omega_(), b1_*F2*sqrt(S2))
     );
 }
 
@@ -174,13 +200,10 @@ tmp<volScalarField::Internal> kOmegaSSTBase<BasicEddyViscosityModel>::GbyNu
 template<class BasicEddyViscosityModel>
 tmp<fvScalarMatrix> kOmegaSSTBase<BasicEddyViscosityModel>::kSource() const
 {
-    return tmp<fvScalarMatrix>
+    return tmp<fvScalarMatrix>::New
     (
-        new fvScalarMatrix
-        (
-            k_,
-            dimVolume*this->rho_.dimensions()*k_.dimensions()/dimTime
-        )
+        k_,
+        dimVolume*this->rho_.dimensions()*k_.dimensions()/dimTime
     );
 }
 
@@ -188,13 +211,10 @@ tmp<fvScalarMatrix> kOmegaSSTBase<BasicEddyViscosityModel>::kSource() const
 template<class BasicEddyViscosityModel>
 tmp<fvScalarMatrix> kOmegaSSTBase<BasicEddyViscosityModel>::omegaSource() const
 {
-    return tmp<fvScalarMatrix>
+    return tmp<fvScalarMatrix>::New
     (
-        new fvScalarMatrix
-        (
-            omega_,
-            dimVolume*this->rho_.dimensions()*omega_.dimensions()/dimTime
-        )
+        omega_,
+        dimVolume*this->rho_.dimensions()*omega_.dimensions()/dimTime
     );
 }
 
@@ -207,13 +227,10 @@ tmp<fvScalarMatrix> kOmegaSSTBase<BasicEddyViscosityModel>::Qsas
     const volScalarField::Internal& beta
 ) const
 {
-    return tmp<fvScalarMatrix>
+    return tmp<fvScalarMatrix>::New
     (
-        new fvScalarMatrix
-        (
-            omega_,
-            dimVolume*this->rho_.dimensions()*omega_.dimensions()/dimTime
-        )
+        omega_,
+        dimVolume*this->rho_.dimensions()*omega_.dimensions()/dimTime
     );
 }
 
@@ -389,6 +406,7 @@ kOmegaSSTBase<BasicEddyViscosityModel>::kOmegaSSTBase
         ),
         this->mesh_
     ),
+
     decayControl_
     (
         Switch::getOrAddToDict
@@ -498,31 +516,30 @@ void kOmegaSSTBase<BasicEddyViscosityModel>::correct()
 
     BasicEddyViscosityModel::correct();
 
-    volScalarField::Internal divU(fvc::div(fvc::absolute(this->phi(), U)));
+    const volScalarField::Internal divU
+    (
+        fvc::div(fvc::absolute(this->phi(), U))
+    );
+
+    const volScalarField CDkOmega
+    (
+        (2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_
+    );
+
+    const volScalarField F1(this->F1(CDkOmega));
+    const volScalarField F23(this->F23());
 
     tmp<volTensorField> tgradU = fvc::grad(U);
-    volScalarField S2(2*magSqr(symm(tgradU())));
-    volScalarField::Internal GbyNu0
-    (
-        this->type() + ":GbyNu",
-        (tgradU() && dev(twoSymm(tgradU())))
-    );
+    const volScalarField S2(this->S2(F1, tgradU()));
+    volScalarField::Internal GbyNu0(this->GbyNu0(tgradU(), F1, S2));
     volScalarField::Internal G(this->GName(), nut*GbyNu0);
 
     // Update omega and G at the wall
     omega_.boundaryFieldRef().updateCoeffs();
 
-    volScalarField CDkOmega
-    (
-        (2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_
-    );
-
-    volScalarField F1(this->F1(CDkOmega));
-    volScalarField F23(this->F23());
-
     {
-        volScalarField::Internal gamma(this->gamma(F1));
-        volScalarField::Internal beta(this->beta(F1));
+        const volScalarField::Internal gamma(this->gamma(F1));
+        const volScalarField::Internal beta(this->beta(F1));
 
         GbyNu0 = GbyNu(GbyNu0, F23(), S2());
 
@@ -555,28 +572,30 @@ void kOmegaSSTBase<BasicEddyViscosityModel>::correct()
         bound(omega_, this->omegaMin_);
     }
 
-    // Turbulent kinetic energy equation
-    tmp<fvScalarMatrix> kEqn
-    (
-        fvm::ddt(alpha, rho, k_)
-      + fvm::div(alphaRhoPhi, k_)
-      - fvm::laplacian(alpha*rho*DkEff(F1), k_)
-     ==
-        alpha()*rho()*Pk(G)
-      - fvm::SuSp((2.0/3.0)*alpha()*rho()*divU, k_)
-      - fvm::Sp(alpha()*rho()*epsilonByk(F1, tgradU()), k_)
-      + alpha()*rho()*betaStar_*omegaInf_*kInf_
-      + kSource()
-      + fvOptions(alpha, rho, k_)
-    );
+    {
+        // Turbulent kinetic energy equation
+        tmp<fvScalarMatrix> kEqn
+        (
+            fvm::ddt(alpha, rho, k_)
+          + fvm::div(alphaRhoPhi, k_)
+          - fvm::laplacian(alpha*rho*DkEff(F1), k_)
+         ==
+            alpha()*rho()*Pk(G)
+          - fvm::SuSp((2.0/3.0)*alpha()*rho()*divU, k_)
+          - fvm::Sp(alpha()*rho()*epsilonByk(F1, tgradU()), k_)
+          + alpha()*rho()*betaStar_*omegaInf_*kInf_
+          + kSource()
+          + fvOptions(alpha, rho, k_)
+        );
 
-    tgradU.clear();
+        tgradU.clear();
 
-    kEqn.ref().relax();
-    fvOptions.constrain(kEqn.ref());
-    solve(kEqn);
-    fvOptions.correct(k_);
-    bound(k_, this->kMin_);
+        kEqn.ref().relax();
+        fvOptions.constrain(kEqn.ref());
+        solve(kEqn);
+        fvOptions.correct(k_);
+        bound(k_, this->kMin_);
+    }
 
     correctNut(S2);
 }

src/TurbulenceModels/turbulenceModels/Base/kOmegaSST/kOmegaSSTBase.H

@@ -6,7 +6,8 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2017-2021 OpenCFD Ltd.
+    Copyright (C) 2022 Upstream CFD GmbH
+    Copyright (C) 2017-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -141,7 +142,7 @@ class kOmegaSSTBase
 
 protected:
 
-    // Protected data
+    // Protected Data
 
         // Model coefficients
 
@@ -174,7 +175,10 @@ protected:
             //  which is for near-wall cells only
             const volScalarField& y_;
 
+            //- Turbulent kinetic energy field [m^2/s^2]
             volScalarField k_;
+
+            //- Specific dissipation rate field [1/s]
             volScalarField omega_;
 
 
@@ -188,6 +192,7 @@ protected:
 
     // Protected Member Functions
 
+        //- Set decay control with kInf and omegaInf
         void setDecayControl(const dictionary& dict);
 
         virtual tmp<volScalarField> F1(const volScalarField& CDkOmega) const;
@@ -195,6 +200,7 @@ protected:
         virtual tmp<volScalarField> F3() const;
         virtual tmp<volScalarField> F23() const;
 
+        //- Return the blended field
         tmp<volScalarField> blend
         (
             const volScalarField& F1,
@@ -205,6 +211,7 @@ protected:
             return F1*(psi1 - psi2) + psi2;
         }
 
+        //- Return the internal blended field
         tmp<volScalarField::Internal> blend
         (
             const volScalarField::Internal& F1,
@@ -253,6 +260,13 @@ protected:
 
         virtual void correctNut();
 
+        //- Return square of strain rate
+        virtual tmp<volScalarField> S2
+        (
+            const volScalarField& F1,
+            const volTensorField& gradU
+        ) const;
+
         //- Return k production rate
         virtual tmp<volScalarField::Internal> Pk
         (
@@ -266,6 +280,14 @@ protected:
             const volTensorField& gradU
         ) const;
 
+        //- Return (G/nu)_0
+        virtual tmp<volScalarField::Internal> GbyNu0
+        (
+            const volTensorField& gradU,
+            const volScalarField& F1,
+            const volScalarField& S2
+        ) const;
+
         //- Return G/nu
         virtual tmp<volScalarField::Internal> GbyNu
         (
@@ -321,22 +343,20 @@ public:
         //- Return the effective diffusivity for k
         tmp<volScalarField> DkEff(const volScalarField& F1) const
         {
-            return tmp<volScalarField>
+            return tmp<volScalarField>::New
             (
-                new volScalarField("DkEff", alphaK(F1)*this->nut_ + this->nu())
+                "DkEff",
+                alphaK(F1)*this->nut_ + this->nu()
             );
         }
 
         //- Return the effective diffusivity for omega
         tmp<volScalarField> DomegaEff(const volScalarField& F1) const
         {
-            return tmp<volScalarField>
+            return tmp<volScalarField>::New
             (
-                new volScalarField
-                (
-                    "DomegaEff",
-                    alphaOmega(F1)*this->nut_ + this->nu()
-                )
+                "DomegaEff",
+                alphaOmega(F1)*this->nut_ + this->nu()
             );
         }