BUG: replced with AMDClang general includes

- broadcast           : (replaces scatter)
  - combineReduce       == combineGather + broadcast
  - listCombineReduce   == listCombineGather + broadcast
  - mapCombineReduce    == mapCombineGather + broadcast
  - allGatherList       == gatherList + scatterList

  Before settling on a more consistent naming convention,
  some intermediate namings were used in OpenFOAM-v2206:

    - combineReduce       (2206: combineAllGather)
    - listCombineReduce   (2206: listCombineAllGather)
    - mapCombineReduce    (2206: mapCombineAllGather)

applications/solvers/combustion/PDRFoam/setDeltaT.H

@@ -36,11 +36,11 @@ Description
 if (adjustTimeStep)
 {
     scalar maxDeltaTFact = maxCo/(CoNum + StCoNum + SMALL);
-    scalar deltaTFact = min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
+    scalar deltaTFact = Foam::min(Foam::min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
 
     runTime.setDeltaT
     (
-        min
+        Foam::min
         (
             deltaTFact*runTime.deltaTValue(),
             maxDeltaT

applications/solvers/combustion/chemFoam/setDeltaT.H

@@ -1,5 +1,5 @@
 if (adjustTimeStep)
 {
-    runTime.setDeltaT(min(dtChem, maxDeltaT));
+    runTime.setDeltaT(Foam::min(dtChem, maxDeltaT));
     Info<< "deltaT = " <<  runTime.deltaT().value() << endl;
 }

applications/solvers/combustion/fireFoam/setMultiRegionDeltaT.H

@@ -54,9 +54,9 @@ if (adjustTimeStep)
 
     runTime.setDeltaT
     (
-        min
+        Foam::min
         (
-            dt0*min(min(TFactorFluid, min(TFactorFilm, TFactorSolid)), 1.2),
+            dt0*Foam::min(Foam::min(TFactorFluid, Foam::min(TFactorFilm, TFactorSolid)), 1.2),
             maxDeltaT
         )
     );

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/compressibleMultiRegionCourantNo.H

@@ -18,13 +18,13 @@
             const surfaceScalarField& phi2 =
                 phaseSystemFluid[regioni].phase2().phiRef();
 
-            sumPhi = max
+            sumPhi = Foam::max
             (
                 sumPhi,
                 fvc::surfaceSum(mag(phi1))().primitiveField()
             );
 
-            sumPhi = max
+            sumPhi = Foam::max
             (
                 sumPhi,
                 fvc::surfaceSum(mag(phi2))().primitiveField()
@@ -43,7 +43,7 @@
               / fluidRegions[regioni].V().field()
             )*runTime.deltaTValue(),
 
-            CoNum = max(UrCoNum, CoNum);
+            CoNum = Foam::max(UrCoNum, CoNum);
         }
     }
 

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H

@@ -2,7 +2,7 @@
 
     forAll(fluidRegions, regioni)
     {
-        CoNum = max
+        CoNum = Foam::max
         (
             compressibleCourantNo
             (

applications/solvers/heatTransfer/chtMultiRegionFoam/include/setInitialMultiRegionDeltaT.H

@@ -47,10 +47,10 @@ if (adjustTimeStep)
 
         runTime.setDeltaT
         (
-            min
+            Foam::min
             (
-                min(maxCo/CoNum, maxDi/DiNum)*runTime.deltaT().value(),
-                min(runTime.deltaTValue(), maxDeltaT)
+                Foam::min(maxCo/CoNum, maxDi/DiNum)*runTime.deltaT().value(),
+                Foam::min(runTime.deltaTValue(), maxDeltaT)
             )
         );
         Info<< "deltaT = " <<  runTime.deltaT().value() << endl;

applications/solvers/heatTransfer/chtMultiRegionFoam/include/setMultiRegionDeltaT.H

@@ -49,17 +49,17 @@ if (adjustTimeStep)
     scalar maxDeltaTSolid = maxDi/(DiNum + SMALL);
 
     scalar deltaTFluid =
-        min
+        Foam::min
         (
-            min(maxDeltaTFluid, 1.0 + 0.1*maxDeltaTFluid),
+            Foam::min(maxDeltaTFluid, 1.0 + 0.1*maxDeltaTFluid),
             1.2
         );
 
     runTime.setDeltaT
     (
-        min
+        Foam::min
         (
-            min(deltaTFluid, maxDeltaTSolid)*runTime.deltaT().value(),
+            Foam::min(deltaTFluid, maxDeltaTSolid)*runTime.deltaT().value(),
             maxDeltaT
         )
     );

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffusionNo.H

@@ -22,7 +22,7 @@ forAll(solidRegions, i)
     tmp<volScalarField> trho = thermo.rho();
     const volScalarField& rho = trho();
 
-    DiNum = max
+    DiNum = Foam::max
     (
         solidRegionDiffNo
         (

applications/solvers/heatTransfer/solidFoam/solidDiffusionNo.H

@@ -17,7 +17,7 @@ scalar DiNum = -GREAT;
     tmp<volScalarField> trho = thermo.rho();
     const volScalarField& rho = trho();
 
-    DiNum = max
+    DiNum = Foam::max
     (
         solidRegionDiffNo
         (

applications/solvers/lagrangian/reactingParcelFoam/setMultiRegionDeltaT.H

@@ -36,13 +36,13 @@ Description
 if (adjustTimeStep)
 {
     const scalar maxDeltaTFact =
-        min(maxCo/(CoNum + SMALL), maxCo/(surfaceFilm.CourantNumber() + SMALL));
+        Foam::min(maxCo/(CoNum + SMALL), maxCo/(surfaceFilm.CourantNumber() + SMALL));
     const scalar deltaTFact =
-        min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
+        Foam::min(Foam::min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
 
     runTime.setDeltaT
     (
-        min
+        Foam::min
         (
             deltaTFact*runTime.deltaTValue(),
             maxDeltaT

applications/solvers/multiphase/MPPICInterFoam/alphaEqn.H

@@ -214,7 +214,7 @@
                 phiCN,
                 alphaPhi10,
                 Sp,
-                (Su + divU*min(alpha1(), scalar(1)))(),
+                (Su + divU*Foam::min(alpha1(), scalar(1)))(),
                 oneField(),
                 zeroField()
             );

applications/solvers/multiphase/VoF/alphaEqn.H

@@ -214,7 +214,7 @@
                 phiCN,
                 alphaPhi10,
                 Sp,
-                (Su + divU*min(alpha1(), scalar(1)))(),
+                (Su + divU*Foam::min(alpha1(), scalar(1)))(),
                 oneField(),
                 zeroField()
             );

applications/solvers/multiphase/VoF/setDeltaT.H

@@ -36,13 +36,13 @@ Description
 if (adjustTimeStep)
 {
     scalar maxDeltaTFact =
-        min(maxCo/(CoNum + SMALL), maxAlphaCo/(alphaCoNum + SMALL));
+        Foam::min(maxCo/(CoNum + SMALL), maxAlphaCo/(alphaCoNum + SMALL));
 
-    scalar deltaTFact = min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
+    scalar deltaTFact = Foam::min(Foam::min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
 
     runTime.setDeltaT
     (
-        min
+        Foam::min
         (
             deltaTFact*runTime.deltaTValue(),
             maxDeltaT

applications/solvers/multiphase/cavitatingFoam/setDeltaT.H

@@ -36,13 +36,13 @@ Description
 if (adjustTimeStep)
 {
     scalar maxDeltaTFact =
-        min(maxCo/(CoNum + SMALL), maxAcousticCo/(acousticCoNum + SMALL));
+        Foam::min(maxCo/(CoNum + SMALL), maxAcousticCo/(acousticCoNum + SMALL));
 
-    scalar deltaTFact = min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
+    scalar deltaTFact = Foam::min(Foam::min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
 
     runTime.setDeltaT
     (
-        min
+        Foam::min
         (
             deltaTFact*runTime.deltaTValue(),
             maxDeltaT

applications/solvers/multiphase/cavitatingFoam/setInitialDeltaT.H

@@ -37,11 +37,11 @@ if (adjustTimeStep)
     if (CoNum > SMALL)
     {
         scalar maxDeltaTFact =
-            min(maxCo/(CoNum + SMALL), maxAcousticCo/(acousticCoNum + SMALL));
+            Foam::min(maxCo/(CoNum + SMALL), maxAcousticCo/(acousticCoNum + SMALL));
 
         runTime.setDeltaT
         (
-            min
+            Foam::min
             (
                 maxDeltaTFact*runTime.deltaTValue(),
                 maxDeltaT

applications/solvers/multiphase/compressibleInterFoam/alphaSuSp.H

@@ -26,12 +26,12 @@ forAll(dgdt, celli)
 {
     if (dgdt[celli] > 0.0)
     {
-        Sp[celli] -= dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
-        Su[celli] += dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
+        Sp[celli] -= dgdt[celli]/Foam::max(1.0 - alpha1[celli], 1e-4);
+        Su[celli] += dgdt[celli]/Foam::max(1.0 - alpha1[celli], 1e-4);
     }
     else if (dgdt[celli] < 0.0)
     {
-        Sp[celli] += dgdt[celli]/max(alpha1[celli], 1e-4);
+        Sp[celli] += dgdt[celli]/Foam::max(alpha1[celli], 1e-4);
     }
 }
 

applications/solvers/multiphase/compressibleInterFoam/compressibleInterIsoFoam/alphaEqn.H

@@ -1,6 +1,6 @@
 // Update alpha1
 #include "alphaSuSp.H"
-advector.advect(Sp,(Su + divU*min(alpha1(), scalar(1)))());
+advector.advect(Sp,(Su + divU*Foam::min(alpha1(), scalar(1)))());
 
 // Update rhoPhi
 rhoPhi = advector.getRhoPhi(rho1, rho2);
@@ -10,6 +10,6 @@ alpha2 = 1.0 - alpha1;
 
 Info<< "Phase-1 volume fraction = "
     << alpha1.weightedAverage(mesh.Vsc()).value()
-    << "  Min(" << alpha1.name() << ") = " << min(alpha1).value()
-    << "  Max(" << alpha1.name() << ") - 1 = " << max(alpha1).value() - 1
+    << "  Min(" << alpha1.name() << ") = " << Foam::min(alpha1).value()
+    << "  Max(" << alpha1.name() << ") - 1 = " << Foam::max(alpha1).value() - 1
     << endl;

applications/solvers/multiphase/compressibleInterFoam/compressibleInterIsoFoam/alphaSuSp.H

@@ -26,12 +26,12 @@ forAll(dgdt, celli)
 {
     if (dgdt[celli] > 0.0)
     {
-        Sp[celli] -= dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
-        Su[celli] += dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
+        Sp[celli] -= dgdt[celli]/Foam::max(1.0 - alpha1[celli], 1e-4);
+        Su[celli] += dgdt[celli]/Foam::max(1.0 - alpha1[celli], 1e-4);
     }
     else if (dgdt[celli] < 0.0)
     {
-        Sp[celli] += dgdt[celli]/max(alpha1[celli], 1e-4);
+        Sp[celli] += dgdt[celli]/Foam::max(alpha1[celli], 1e-4);
     }
 }
 

applications/solvers/multiphase/compressibleInterFoam/overCompressibleInterDyMFoam/alphaSuSp.H

@@ -26,12 +26,12 @@ forAll(dgdt, celli)
 {
     if (dgdt[celli] > 0.0)
     {
-        Sp[celli] -= dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
-        Su[celli] += dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
+        Sp[celli] -= dgdt[celli]/Foam::max(1.0 - alpha1[celli], 1e-4);
+        Su[celli] += dgdt[celli]/Foam::max(1.0 - alpha1[celli], 1e-4);
     }
     else if (dgdt[celli] < 0.0)
     {
-        Sp[celli] += dgdt[celli]/max(alpha1[celli], 1e-4);
+        Sp[celli] += dgdt[celli]/Foam::max(alpha1[celli], 1e-4);
     }
 }
 

applications/solvers/multiphase/icoReactingMultiphaseInterFoam/setDeltaT.H

@@ -36,13 +36,13 @@ Description
 if (adjustTimeStep)
 {
     scalar maxDeltaTFact =
-        min
+        Foam::min
         (
             maxCo/(CoNum + SMALL),
-            min
+            Foam::min
             (
                 maxAlphaCo/(alphaCoNum + SMALL),
-                min
+                Foam::min
                 (
                     maxAlphaDdt/(ddtAlphaNum + SMALL),
                     maxDi/(DiNum + SMALL)
@@ -50,11 +50,11 @@ if (adjustTimeStep)
             )
         );
 
-    scalar deltaTFact = min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
+    scalar deltaTFact = Foam::min(Foam::min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
 
     runTime.setDeltaT
     (
-        min
+        Foam::min
         (
             deltaTFact*runTime.deltaTValue(),
             maxDeltaT

applications/solvers/multiphase/reactingTwoPhaseEulerFoam/CourantNos.H

@@ -8,5 +8,5 @@
 
     Info<< "Max Ur Courant Number = " << UrCoNum << endl;
 
-    CoNum = max(CoNum, UrCoNum);
+    CoNum = Foam::max(CoNum, UrCoNum);
 }

applications/solvers/multiphase/twoPhaseEulerFoam/CourantNos.H

@@ -8,5 +8,5 @@
 
     Info<< "Max Ur Courant Number = " << UrCoNum << endl;
 
-    CoNum = max(CoNum, UrCoNum);
+    CoNum = Foam::max(CoNum, UrCoNum);
 }

applications/utilities/mesh/advanced/snappyRefineMesh/snappyRefineMesh.C

@@ -142,26 +142,26 @@ scalar getEdgeStats(const primitiveMesh& mesh, const direction excludeCmpt)
 
         if (mag(eVec & x) > 1-edgeTol)
         {
-            minX = min(minX, eMag);
-            maxX = max(maxX, eMag);
+            minX = Foam::min(minX, eMag);
+            maxX = Foam::max(maxX, eMag);
             nX++;
         }
         else if (mag(eVec & y) > 1-edgeTol)
         {
-            minY = min(minY, eMag);
-            maxY = max(maxY, eMag);
+            minY = Foam::min(minY, eMag);
+            maxY = Foam::max(maxY, eMag);
             nY++;
         }
         else if (mag(eVec & z) > 1-edgeTol)
         {
-            minZ = min(minZ, eMag);
-            maxZ = max(maxZ, eMag);
+            minZ = Foam::min(minZ, eMag);
+            maxZ = Foam::max(maxZ, eMag);
             nZ++;
         }
         else
         {
-            minOther = min(minOther, eMag);
-            maxOther = max(maxOther, eMag);
+            minOther = Foam::min(minOther, eMag);
+            maxOther = Foam::max(maxOther, eMag);
         }
     }
 
@@ -179,19 +179,19 @@ scalar getEdgeStats(const primitiveMesh& mesh, const direction excludeCmpt)
 
     if (excludeCmpt == 0)
     {
-        return min(minY, min(minZ, minOther));
+        return Foam::min(minY, Foam::min(minZ, minOther));
     }
     else if (excludeCmpt == 1)
     {
-        return min(minX, min(minZ, minOther));
+        return Foam::min(minX, Foam::min(minZ, minOther));
     }
     else if (excludeCmpt == 2)
     {
-        return min(minX, min(minY, minOther));
+        return Foam::min(minX, Foam::min(minY, minOther));
     }
     else
     {
-        return min(minX, min(minY, min(minZ, minOther)));
+        return Foam::min(minX, Foam::min(minY, Foam::min(minZ, minOther)));
     }
 }
 
@@ -771,7 +771,7 @@ int main(int argc, char *argv[])
     {
         Info<< "Read existing refinement level from file "
             << refLevel.objectPath() << nl
-            << "   min level : " << min(refLevel) << nl
+            << "   min level : " << Foam::min(refLevel) << nl
             << "   max level : " << maxLevel << nl
             << endl;
     }

applications/utilities/mesh/conversion/cfx4ToFoam/cfx4ToFoam.C

@@ -271,7 +271,7 @@ int main(int argc, char *argv[])
                     if (blockPFacePointi != blockPFacePointi2)
                     {
                         sqrMergeTol =
-                            min
+                            Foam::min
                             (
                                 sqrMergeTol,
                                 magSqr
@@ -338,16 +338,16 @@ int main(int argc, char *argv[])
                         blockNFacePoints[blockNFacePointi]
                       + blockOffsets[blockNlabel];
 
-                    label minPN = min(PpointLabel, NpointLabel);
+                    label minPN = Foam::min(PpointLabel, NpointLabel);
 
                     if (pointMergeList[PpointLabel] != -1)
                     {
-                        minPN = min(minPN, pointMergeList[PpointLabel]);
+                        minPN = Foam::min(minPN, pointMergeList[PpointLabel]);
                     }
 
                     if (pointMergeList[NpointLabel] != -1)
                     {
-                        minPN = min(minPN, pointMergeList[NpointLabel]);
+                        minPN = Foam::min(minPN, pointMergeList[NpointLabel]);
                     }
 
                     pointMergeList[PpointLabel]
@@ -411,7 +411,7 @@ int main(int argc, char *argv[])
 
                         pointMergeList[PpointLabel]
                       = pointMergeList[NpointLabel]
-                      = min
+                      = Foam::min
                         (
                             pointMergeList[PpointLabel],
                             pointMergeList[NpointLabel]
@@ -757,7 +757,7 @@ int main(int argc, char *argv[])
     );
 
     // Set the precision of the points data to 10
-    IOstream::defaultPrecision(max(10u, IOstream::defaultPrecision()));
+    IOstream::defaultPrecision(Foam::max(10u, IOstream::defaultPrecision()));
 
     Info<< "Writing polyMesh" << endl;
     pShapeMesh.removeFiles();

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

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

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

@@ -1,11 +0,0 @@
-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

@@ -1,194 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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

@@ -1,118 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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/conversion/gambitToFoam/gambitToFoam.L

@@ -350,7 +350,7 @@ mtype                 {space}"MTYPE:"{space}
         // Find out how many labels are expected. If less or equal to
         // seven, read them all and finish with it. If there is more,
         // set read of the next line
-        label labelsToRead = min(8, nVertices);
+        label labelsToRead = Foam::min(8, nVertices);
         label labelI = 0;
         for (; labelI < labelsToRead; labelI++)
         {
@@ -387,7 +387,7 @@ mtype                 {space}"MTYPE:"{space}
 
         labelList& curLabels = cellLabels[curNumberOfCells];
 
-        label labelsToRead = min
+        label labelsToRead = Foam::min
         (
             (nCellContinuationLines + 1)*7,
             curLabels.size()
@@ -869,7 +869,7 @@ int main(int argc, char *argv[])
     );
 
     // Set the precision of the points data to 10
-    IOstream::defaultPrecision(max(10u, IOstream::defaultPrecision()));
+    IOstream::defaultPrecision(Foam::max(10u, IOstream::defaultPrecision()));
 
     Info<< "Writing polyMesh" << endl;
     pShapeMesh.removeFiles();

applications/utilities/mesh/conversion/ideasUnvToFoam/ideasUnvToFoam.C

@@ -269,7 +269,7 @@ void readCells
     label maxUnvPoint = 0;
     forAll(unvPointID, pointi)
     {
-        maxUnvPoint = max(maxUnvPoint, unvPointID[pointi]);
+        maxUnvPoint = Foam::max(maxUnvPoint, unvPointID[pointi]);
     }
     labelList unvToFoam(invert(maxUnvPoint+1, unvPointID));
 
@@ -784,7 +784,7 @@ int main(int argc, char *argv[])
     label maxUnvPoint = 0;
     forAll(unvPointID, pointi)
     {
-        maxUnvPoint = max(maxUnvPoint, unvPointID[pointi]);
+        maxUnvPoint = Foam::max(maxUnvPoint, unvPointID[pointi]);
     }
     labelList unvToFoam(invert(maxUnvPoint+1, unvPointID));
 
@@ -1281,7 +1281,7 @@ int main(int argc, char *argv[])
     }
 
     // Set the precision of the points data to 10
-    IOstream::defaultPrecision(max(10u, IOstream::defaultPrecision()));
+    IOstream::defaultPrecision(Foam::max(10u, IOstream::defaultPrecision()));
 
     mesh.write();
 

applications/utilities/mesh/conversion/kivaToFoam/checkPatch.H

@@ -8,10 +8,10 @@
         {
             if (kivaVersion == kiva3v)
             {
-                regionIndex = max
+                regionIndex = Foam::max
                 (
-                    max(idface[quadFace[0]], idface[quadFace[1]]),
-                    max(idface[quadFace[2]], idface[quadFace[3]])
+                    Foam::max(idface[quadFace[0]], idface[quadFace[1]]),
+                    Foam::max(idface[quadFace[2]], idface[quadFace[3]])
                 );
 
                 if (regionIndex > 0)

applications/utilities/mesh/conversion/kivaToFoam/readKivaGrid.H

@@ -148,7 +148,7 @@ for (label i=0; i<nPoints; i++)
                     end = pointMap[end];
                 }
 
-                label minLabel = min(start, end);
+                label minLabel = Foam::min(start, end);
 
                 pointMap[start] = pointMap[end] = minLabel;
             }
@@ -331,7 +331,7 @@ if
     {
         forAll(pf, pfi)
         {
-            minz = min(minz, points[pf[pfi]].z());
+            minz = Foam::min(minz, points[pf[pfi]].z());
         }
     }
 
@@ -344,7 +344,7 @@ if
         scalar minfz = GREAT;
         forAll(pf, pfi)
         {
-            minfz = min(minfz, points[pf[pfi]].z());
+            minfz = Foam::min(minfz, points[pf[pfi]].z());
         }
 
         if (minfz > minz)
@@ -371,7 +371,7 @@ if
         scalar minfz = GREAT;
         forAll(pf, pfi)
         {
-            minfz = min(minfz, points[pf[pfi]].z());
+            minfz = Foam::min(minfz, points[pf[pfi]].z());
         }
 
         if (minfz < zHeadMin)
@@ -570,7 +570,7 @@ polyMesh pShapeMesh
 );
 
 // Set the precision of the points data to 10
-IOstream::defaultPrecision(max(10u, IOstream::defaultPrecision()));
+IOstream::defaultPrecision(Foam::max(10u, IOstream::defaultPrecision()));
 
 Info << "Writing polyMesh" << endl;
 pShapeMesh.removeFiles();

applications/utilities/mesh/conversion/netgenNeutralToFoam/netgenNeutralToFoam.C

@@ -188,7 +188,7 @@ int main(int argc, char *argv[])
         }
 
 
-        maxPatch = max(maxPatch, patchi);
+        maxPatch = Foam::max(maxPatch, patchi);
 
         triFace tri(readLabel(str)-1, readLabel(str)-1, readLabel(str)-1);
 
@@ -319,7 +319,7 @@ int main(int argc, char *argv[])
     );
 
     // Set the precision of the points data to 10
-    IOstream::defaultPrecision(max(10u, IOstream::defaultPrecision()));
+    IOstream::defaultPrecision(Foam::max(10u, IOstream::defaultPrecision()));
 
     Info<< "Writing mesh ..." << endl;
     mesh.removeFiles();

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

@@ -556,8 +556,8 @@ void calcEdgeMinMaxZone
             forAll(eFaces, i)
             {
                 label zoneI = mappedZoneID[eFaces[i]];
-                minZoneID[edgeI] = min(minZoneID[edgeI], zoneI);
-                maxZoneID[edgeI] = max(maxZoneID[edgeI], zoneI);
+                minZoneID[edgeI] = Foam::min(minZoneID[edgeI], zoneI);
+                maxZoneID[edgeI] = Foam::max(maxZoneID[edgeI], zoneI);
             }
         }
     }
@@ -813,8 +813,8 @@ void addCoupledPatches
             forAll(eFaces, i)
             {
                 label proci = procID[eFaces[i]];
-                minProcID[edgeI] = min(minProcID[edgeI], proci);
-                maxProcID[edgeI] = max(maxProcID[edgeI], proci);
+                minProcID[edgeI] = Foam::min(minProcID[edgeI], proci);
+                maxProcID[edgeI] = Foam::max(maxProcID[edgeI], proci);
             }
         }
     }
@@ -1291,7 +1291,7 @@ void extrudeGeometricProperties
             label celli = regionMesh.faceOwner()[facei];
             if (regionMesh.isInternalFace(facei))
             {
-                celli = max(celli, regionMesh.faceNeighbour()[facei]);
+                celli = Foam::max(celli, regionMesh.faceNeighbour()[facei]);
             }
 
             // Calculate layer from cell numbering (see createShellMesh)
@@ -2192,8 +2192,8 @@ int main(int argc, char *argv[])
 
             if (zone0 != zone1) // || (cos(angle) > blabla))
             {
-                label minZone = min(zone0,zone1);
-                label maxZone = max(zone0,zone1);
+                label minZone = Foam::min(zone0,zone1);
+                label maxZone = Foam::max(zone0,zone1);
                 label index = minZone*zoneNames.size()+maxZone;
 
                 ePatches.setSize(eFaces.size());

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

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

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

@@ -27,7 +27,7 @@ void maxFaceToCell
         const cell& cFaces = cells[cellI];
         forAll(cFaces, i)
         {
-            cellFld[cellI] = max(cellFld[cellI], faceData[cFaces[i]]);
+            cellFld[cellI] = Foam::max(cellFld[cellI], faceData[cFaces[i]]);
         }
     }
 
@@ -57,7 +57,7 @@ void minFaceToCell
         const cell& cFaces = cells[cellI];
         forAll(cFaces, i)
         {
-            cellFld[cellI] = min(cellFld[cellI], faceData[cFaces[i]]);
+            cellFld[cellI] = Foam::min(cellFld[cellI], faceData[cFaces[i]]);
         }
     }
 
@@ -88,8 +88,8 @@ void minFaceToCell
     // Internal faces
     forAll(own, facei)
     {
-        cellFld[own[facei]] = min(cellFld[own[facei]], faceData[facei]);
-        cellFld[nei[facei]] = min(cellFld[nei[facei]], faceData[facei]);
+        cellFld[own[facei]] = Foam::min(cellFld[own[facei]], faceData[facei]);
+        cellFld[nei[facei]] = Foam::min(cellFld[nei[facei]], faceData[facei]);
     }
 
     // Patch faces
@@ -100,7 +100,7 @@ void minFaceToCell
 
         forAll(fc, i)
         {
-            cellFld[fc[i]] = min(cellFld[fc[i]], fvp[i]);
+            cellFld[fc[i]] = Foam::min(cellFld[fc[i]], fvp[i]);
         }
     }
 
@@ -180,7 +180,7 @@ void Foam::writeFields
         (
             radToDeg
             (
-                Foam::acos(min(scalar(1), max(scalar(-1), faceOrthogonality)))
+                Foam::acos(Foam::min(scalar(1), Foam::max(scalar(-1), faceOrthogonality)))
             )
         );
 
@@ -540,7 +540,7 @@ void Foam::writeFields
                         ownCc,
                         fc
                     ).quality();
-                    ownVol = min(ownVol, tetQual);
+                    ownVol = Foam::min(ownVol, tetQual);
                 }
             }
             if (mesh.isInternalFace(facei))
@@ -556,7 +556,7 @@ void Foam::writeFields
                         fc,
                         neiCc
                     ).quality();
-                    neiVol = min(neiVol, tetQual);
+                    neiVol = Foam::min(neiVol, tetQual);
                 }
             }
         }
@@ -608,8 +608,8 @@ void Foam::writeFields
         // Internal faces
         forAll(own, facei)
         {
-            cellFld[own[facei]] = min(cellFld[own[facei]], ownPyrVol[facei]);
-            cellFld[nei[facei]] = min(cellFld[nei[facei]], neiPyrVol[facei]);
+            cellFld[own[facei]] = Foam::min(cellFld[own[facei]], ownPyrVol[facei]);
+            cellFld[nei[facei]] = Foam::min(cellFld[nei[facei]], neiPyrVol[facei]);
         }
 
         // Patch faces
@@ -620,7 +620,7 @@ void Foam::writeFields
             forAll(fc, i)
             {
                 const label meshFacei = fvp.patch().start();
-                cellFld[fc[i]] = min(cellFld[fc[i]], ownPyrVol[meshFacei]);
+                cellFld[fc[i]] = Foam::min(cellFld[fc[i]], ownPyrVol[meshFacei]);
             }
         }
 
@@ -631,7 +631,7 @@ void Foam::writeFields
         if (writeFaceFields)
         {
             scalarField minFacePyrVol(neiPyrVol);
-            minFacePyrVol = min
+            minFacePyrVol = Foam::min
             (
                 minFacePyrVol,
                 SubField<scalar>(ownPyrVol, mesh.nInternalFaces())

applications/utilities/mesh/manipulation/refineMesh/refineMesh.C

@@ -100,26 +100,26 @@ void printEdgeStats(const polyMesh& mesh)
 
             if (mag(eVec & x) > 1-edgeTol)
             {
-                minX = min(minX, eMag);
-                maxX = max(maxX, eMag);
+                minX = Foam::min(minX, eMag);
+                maxX = Foam::max(maxX, eMag);
                 nX++;
             }
             else if (mag(eVec & y) > 1-edgeTol)
             {
-                minY = min(minY, eMag);
-                maxY = max(maxY, eMag);
+                minY = Foam::min(minY, eMag);
+                maxY = Foam::max(maxY, eMag);
                 nY++;
             }
             else if (mag(eVec & z) > 1-edgeTol)
             {
-                minZ = min(minZ, eMag);
-                maxZ = max(maxZ, eMag);
+                minZ = Foam::min(minZ, eMag);
+                maxZ = Foam::max(maxZ, eMag);
                 nZ++;
             }
             else
             {
-                minOther = min(minOther, eMag);
-                maxOther = max(maxOther, eMag);
+                minOther = Foam::min(minOther, eMag);
+                maxOther = Foam::max(maxOther, eMag);
             }
         }
     }

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

@@ -145,10 +145,10 @@ void getBand
 
         // Note: mag not necessary for correct (upper-triangular) ordering.
         label diff = nei-own;
-        cellBandwidth[nei] = max(cellBandwidth[nei], diff);
+        cellBandwidth[nei] = Foam::max(cellBandwidth[nei], diff);
     }
 
-    bandwidth = max(cellBandwidth);
+    bandwidth = Foam::max(cellBandwidth);
 
     // Do not use field algebra because of conversion label to scalar
     profile = 0.0;
@@ -340,7 +340,7 @@ labelList getRegionFaceOrder
     }
 
     // Do region interfaces
-    label nRegions = max(cellToRegion)+1;
+    label nRegions = Foam::max(cellToRegion)+1;
     {
         // Sort in increasing region
         SortableList<label> sortKey(mesh.nFaces(), labelMax);
@@ -353,8 +353,8 @@ labelList getRegionFaceOrder
             if (ownRegion != neiRegion)
             {
                 sortKey[facei] =
-                    min(ownRegion, neiRegion)*nRegions
-                   +max(ownRegion, neiRegion);
+                    Foam::min(ownRegion, neiRegion)*nRegions
+                   +Foam::max(ownRegion, neiRegion);
             }
         }
         sortKey.sort();
@@ -566,7 +566,7 @@ labelList regionRenumber
 
     labelList cellOrder(cellToRegion.size());
 
-    label nRegions = max(cellToRegion)+1;
+    label nRegions = Foam::max(cellToRegion)+1;
 
     labelListList regionToCells(invertOneToMany(nRegions, cellToRegion));
 

applications/utilities/mesh/manipulation/setSet/setSet.C

@@ -320,10 +320,10 @@ bool doCommand
     const globalMeshData& parData = mesh.globalData();
 
     label typSize =
-        max
+        Foam::max
         (
             parData.nTotalCells(),
-            max
+            Foam::max
             (
                 parData.nTotalFaces(),
                 parData.nTotalPoints()
@@ -375,7 +375,7 @@ bool doCommand
 
                 topoSet& currentSet = currentSetPtr();
                 // Presize it according to current mesh data.
-                currentSet.resize(max(currentSet.size(), typSize));
+                currentSet.resize(Foam::max(currentSet.size(), typSize));
             }
         }
 

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

@@ -286,8 +286,8 @@ void addToInterface
 {
     edge interface
     (
-        min(ownRegion, neiRegion),
-        max(ownRegion, neiRegion)
+        Foam::min(ownRegion, neiRegion),
+        Foam::max(ownRegion, neiRegion)
     );
 
     auto iter = regionsToSize.find(interface);
@@ -544,8 +544,8 @@ void getInterfaceSizes
 
             edge interface
             (
-                min(ownRegion, neiRegion),
-                max(ownRegion, neiRegion)
+                Foam::min(ownRegion, neiRegion),
+                Foam::max(ownRegion, neiRegion)
             );
 
             faceToInterface[facei] = regionsToInterface[interface][zoneID];
@@ -567,8 +567,8 @@ void getInterfaceSizes
 
             edge interface
             (
-                min(ownRegion, neiRegion),
-                max(ownRegion, neiRegion)
+                Foam::min(ownRegion, neiRegion),
+                Foam::max(ownRegion, neiRegion)
             );
 
             faceToInterface[facei] = regionsToInterface[interface][zoneID];

applications/utilities/mesh/manipulation/subsetMesh/subsetMesh.C

@@ -345,7 +345,7 @@ void subsetTopoSets
 
         Info<< "Subsetting " << set.type() << " " << set.name() << endl;
 
-        labelHashSet subset(2*min(set.size(), map.size()));
+        labelHashSet subset(2*Foam::min(set.size(), map.size()));
 
         forAll(map, i)
         {

applications/utilities/parallelProcessing/reconstructParMesh/reconstructParMesh.C

@@ -977,7 +977,7 @@ int main(int argc, char *argv[])
                         for
                         (
                             label addedI=next;
-                            addedI<min(proci+step, nProcs);
+                            addedI<Foam::min(proci+step, nProcs);
                             addedI++
                         )
                         {

applications/utilities/parallelProcessing/redistributePar/redistributePar.C

@@ -334,11 +334,11 @@ void printMeshData(const polyMesh& mesh)
                 << nBndFaces-nProcFaces << endl;
         }
 
-        maxProcCells = max(maxProcCells, nLocalCells);
+        maxProcCells = Foam::max(maxProcCells, nLocalCells);
         totProcFaces += nProcFaces;
         totProcPatches += nei.size();
-        maxProcFaces = max(maxProcFaces, nProcFaces);
-        maxProcPatches = max(maxProcPatches, nei.size());
+        maxProcFaces = Foam::max(maxProcFaces, nProcFaces);
+        maxProcPatches = Foam::max(maxProcPatches, nei.size());
     }
 
     // Summary stats

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

@@ -131,7 +131,7 @@ int main(int argc, char *argv[])
     args.readIfPresent("format", setFormat);
 
     args.readIfPresent("stride", sampleFrequency);
-    sampleFrequency = max(1, sampleFrequency);  // sanity
+    sampleFrequency = Foam::max(1, sampleFrequency);  // sanity
 
     // Setup the writer
     auto writerPtr =
@@ -179,7 +179,7 @@ int main(int argc, char *argv[])
                     maxIds.resize(origProc+1, -1);
                 }
 
-                maxIds[origProc] = max(maxIds[origProc], origId);
+                maxIds[origProc] = Foam::max(maxIds[origProc], origId);
             }
         }
 

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

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

applications/utilities/preProcessing/PDR/pdrFields/PDRarraysCalc.C

@@ -1056,7 +1056,7 @@ void calc_drag_etc
     const scalar expon =
     (
         br > 0.0
-      ? min(max((surr_br / br - 0.25) * 4.0 / 3.0, scalar(0)), scalar(1))
+      ? Foam::min(Foam::max((surr_br / br - 0.25) * 4.0 / 3.0, scalar(0)), scalar(1))
       : 0.0
     );
 
@@ -1114,16 +1114,16 @@ void Foam::PDRarrays::blockageSummary() const
                 totVolBlock += v_block(ijk) * pdrBlock.V(ijk);
                 totArea += surf(ijk);
 
-                totCount += max(0, obs_count(ijk));
+                totCount += Foam::max(0, obs_count(ijk));
 
-                totDrag.x() += max(0, drag_s(ijk).xx());
-                totDrag.y() += max(0, drag_s(ijk).yy());
-                totDrag.z() += max(0, drag_s(ijk).zz());
+                totDrag.x() += Foam::max(0, drag_s(ijk).xx());
+                totDrag.y() += Foam::max(0, drag_s(ijk).yy());
+                totDrag.z() += Foam::max(0, drag_s(ijk).zz());
 
                 for (direction cmpt=0; cmpt < vector::nComponents; ++cmpt)
                 {
-                    totBlock[cmpt] += max(0, area_block_s(ijk)[cmpt]);
-                    totBlock[cmpt] += max(0, area_block_r(ijk)[cmpt]);
+                    totBlock[cmpt] += Foam::max(0, area_block_s(ijk)[cmpt]);
+                    totBlock[cmpt] += Foam::max(0, area_block_r(ijk)[cmpt]);
                 }
             }
         }

applications/utilities/preProcessing/PDR/pdrFields/PDRutilsOverlap.C

@@ -302,7 +302,7 @@ void Foam::PDRutils::circle_overlap
                     scalar da = ac - 0.5 * (a1 + a2);
                     scalar db = bc - 0.5 * (b1 + b2);
                     scalar dc = std::hypot(da, db);
-                    scalar rat1 = min(max((dc / sqrt(area) - 0.3) * 1.4, 0), 1);
+                    scalar rat1 = Foam::min(Foam::max((dc / sqrt(area) - 0.3) * 1.4, 0), 1);
                     scalar drg0 = c_drag(ia,ib).xx();
                     scalar drg1 = c_drag(ia,ib).yy();
                     scalar drg = std::hypot(drg0, drg1);
@@ -449,8 +449,8 @@ scalar block_overlap
             {
                 PDRobstacle over;
 
-                over.pt = max(blk1.pt, blk2.pt);
-                over.span = min(max1, max2) - over.pt;
+                over.pt = Foam::max(blk1.pt, blk2.pt);
+                over.span = Foam::min(max1, max2) - over.pt;
 
                 assert(cmptProduct(over.span) > 0.0);
 
@@ -603,11 +603,11 @@ scalar block_cylinder_overlap
 
                     over.x() = a_centre - 0.5 * a_lblk;
                     over.y() = b_centre - 0.5 * b_lblk;
-                    over.z() = max(blk1.z(), cyl2.z());
+                    over.z() = Foam::max(blk1.z(), cyl2.z());
 
                     over.span.x() = a_lblk;
                     over.span.y() = b_lblk;
-                    over.span.z() = min(max1.z(), cyl2.z() + cyl2.len()) - over.z();
+                    over.span.z() = Foam::min(max1.z(), cyl2.z() + cyl2.len()) - over.z();
                     assert(over.x() > -200.0);
                     assert(over.x() < 2000.0);
                 }
@@ -668,11 +668,11 @@ scalar block_cylinder_overlap
 
                     over.z() = a_centre - a_lblk * 0.5;
                     over.x() = b_centre - b_lblk * 0.5;
-                    over.y() = max(blk1.y(), cyl2.y());
+                    over.y() = Foam::max(blk1.y(), cyl2.y());
 
                     over.span.z() = a_lblk;
                     over.span.x() = b_lblk;
-                    over.span.y() = min(max1.y(), cyl2.y() + cyl2.len()) - over.y();
+                    over.span.y() = Foam::min(max1.y(), cyl2.y() + cyl2.len()) - over.y();
                 }
                 break;
 
@@ -734,11 +734,11 @@ scalar block_cylinder_overlap
 
                     over.y() = a_centre - a_lblk * 0.5;
                     over.z() = b_centre - b_lblk * 0.5;
-                    over.x() = max(blk1.x(), cyl2.x());
+                    over.x() = Foam::max(blk1.x(), cyl2.x());
 
                     over.span.y() = a_lblk;
                     over.span.z() = b_lblk;
-                    over.span.x() = min(max1.x(), cyl2.x() + cyl2.len()) - over.x();
+                    over.span.x() = Foam::min(max1.x(), cyl2.x() + cyl2.len()) - over.x();
                 }
                 break;
             }

applications/utilities/preProcessing/engineSwirl/engineSwirl.C

@@ -71,9 +71,9 @@ int main(int argc, char *argv[])
         {
             scalar b = j1(swirlProfile*r/cylinderRadius).value();
             scalar vEff = omega*b;
-            r = max(r, SMALL);
+            r = Foam::max(r, SMALL);
             U[celli] = ((vEff/r)*(c & yT))*xT + (-(vEff/r)*(c & xT))*yT;
-            Umax = max(Umax, mag(U[celli]));
+            Umax = Foam::max(Umax, mag(U[celli]));
         }
     }
 

applications/utilities/surface/surfaceCheck/surfaceCheck.C

@@ -889,7 +889,7 @@ int main(int argc, char *argv[])
             writeParts
             (
                 surf,
-                min(outputThreshold, numZones),
+                Foam::min(outputThreshold, numZones),
                 faceZone,
                 surfFilePath,
                 surfFileStem
@@ -955,7 +955,7 @@ int main(int argc, char *argv[])
             writeParts
             (
                 surf,
-                min(outputThreshold, numNormalZones),
+                Foam::min(outputThreshold, numNormalZones),
                 normalZone,
                 surfFilePath,
                 surfFileStem + "_normal"

applications/utilities/surface/surfaceHookUp/surfaceHookUp.C

@@ -290,7 +290,7 @@ int main(int argc, char *argv[])
     const IOdictionary dict(dictIO);
 
     const scalar dist(args.get<scalar>(1));
-    const scalar matchTolerance(max(1e-6*dist, SMALL));
+    const scalar matchTolerance(Foam::max(1e-6*dist, SMALL));
     const label maxIters = 100;
 
     Info<< "Hooking distance = " << dist << endl;

applications/utilities/surface/surfaceInflate/surfaceInflate.C

@@ -275,7 +275,7 @@ label detectIntersectionPoints
 
     // 1. Extrusion offset vectors intersecting new surface location
     {
-        scalar tol = max(tolerance, 10*s.tolerance());
+        scalar tol = Foam::max(tolerance, 10*s.tolerance());
 
         // Collect all the edge vectors. Slightly shorten the edges to prevent
         // finding lots of intersections. The fast triangle intersection routine
@@ -295,7 +295,7 @@ label detectIntersectionPoints
             &&  !localFaces[hits[pointI].index()].found(pointI)
             )
             {
-                scale[pointI] = max(0.0, scale[pointI]-0.2);
+                scale[pointI] = Foam::max(0.0, scale[pointI]-0.2);
 
                 isPointOnHitEdge.set(pointI);
                 nHits++;
@@ -328,7 +328,7 @@ label detectIntersectionPoints
                         << pt
                         << endl;
 
-                    scale[e[0]] = max(0.0, scale[e[0]]-0.2);
+                    scale[e[0]] = Foam::max(0.0, scale[e[0]]-0.2);
                     nHits++;
                 }
                 if (isPointOnHitEdge.set(e[1]))
@@ -340,7 +340,7 @@ label detectIntersectionPoints
                         << pt
                         << endl;
 
-                    scale[e[1]] = max(0.0, scale[e[1]]-0.2);
+                    scale[e[1]] = Foam::max(0.0, scale[e[1]]-0.2);
                     nHits++;
                 }
             }
@@ -416,7 +416,7 @@ void minSmooth
         const edge& e = edges[edgeI];
         scalar w = mag(points[mp[e[0]]]-points[mp[e[1]]]);
 
-        edgeWeights[edgeI] = 1.0/(max(w, SMALL));
+        edgeWeights[edgeI] = 1.0/(Foam::max(w, SMALL));
     }
 
     tmp<scalarField> tavgFld = avg(s, fld, edgeWeights);
@@ -427,7 +427,7 @@ void minSmooth
     {
         if (isAffectedPoint.test(pointI))
         {
-            newFld[pointI] = min
+            newFld[pointI] = Foam::min
             (
                 fld[pointI],
                 0.5*fld[pointI] + 0.5*avgFld[pointI]

applications/utilities/thermophysical/adiabaticFlameT/adiabaticFlameT.C

@@ -180,8 +180,8 @@ int main(int argc, char *argv[])
 
         scalar o2 = (1.0/equiv)*stoicO2;
         scalar n2 = (0.79/0.21)*o2;
-        scalar fres = max(1.0 - 1.0/equiv, 0.0);
-        scalar ores = max(1.0/equiv - 1.0, 0.0);
+        scalar fres = Foam::max(1.0 - 1.0/equiv, 0.0);
+        scalar ores = Foam::max(1.0/equiv - 1.0, 0.0);
         scalar fburnt = 1.0 - fres;
 
         thermo reactants

applications/utilities/thermophysical/equilibriumFlameT/equilibriumFlameT.C

@@ -196,12 +196,12 @@ int main(int argc, char *argv[])
         // Number of moles of species for one mole of fuel
         scalar o2 = (1.0/equiv)*stoicO2;
         scalar n2 = (0.79/0.21)*o2;
-        scalar fres = max(1.0 - 1.0/equiv, 0.0);
+        scalar fres = Foam::max(1.0 - 1.0/equiv, 0.0);
         scalar fburnt = 1.0 - fres;
 
         // Initial guess for number of moles of product species
         // ignoring product dissociation
-        scalar oresInit = max(1.0/equiv - 1.0, 0.0)*stoicO2;
+        scalar oresInit = Foam::max(1.0/equiv - 1.0, 0.0)*stoicO2;
         scalar co2Init = fburnt*stoicCO2;
         scalar h2oInit = fburnt*stoicH2O;
 
@@ -231,18 +231,18 @@ int main(int argc, char *argv[])
             if (j > 0)
             {
                 co = co2*
-                    min
+                    Foam::min
                     (
                         CO2BreakUp.Kn(P, equilibriumFlameTemperature, N)
-                       /::sqrt(max(ores, 0.001)),
+                       /::sqrt(Foam::max(ores, 0.001)),
                         1.0
                     );
 
                 h2 = h2o*
-                    min
+                    Foam::min
                     (
                         H2OBreakUp.Kn(P, equilibriumFlameTemperature, N)
-                       /::sqrt(max(ores, 0.001)),
+                       /::sqrt(Foam::max(ores, 0.001)),
                         1.0
                     );
 

src/OpenFOAM/containers/Lists/List/List.C

@@ -49,7 +49,7 @@ void Foam::List<T>::doResize(const label len)
         // With sign-check to avoid spurious -Walloc-size-larger-than
         T* nv = new T[len];
 
-        const label overlap = min(this->size_, len);
+        const label overlap = Foam::min(this->size_, len);
 
         if (overlap)
         {

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

@@ -327,7 +327,7 @@ public:
         //  \return True on success
         static bool write
         (
-            const UPstream::commsTypes commsType,
+            const 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
 {
-    //- Registered reader for UPstream::defaultCommsType
+    // Register re-reader
     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"    : (MPI_Bsend, MPI_Recv)
-        scheduled,      //!< "scheduled"   : (MPI_Send, MPI_Recv)
-        nonBlocking     //!< "nonBlocking" : (MPI_Isend, MPI_Irecv)
+        blocking,       //!< "blocking"
+        scheduled,      //!< "scheduled"
+        nonBlocking     //!< "nonBlocking"
     };
 
-    //- Enumerated names for the communication types
+    //- Names of the communication types
     static const Enum<commsTypes> commsTypeNames;
 
 

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

@@ -221,48 +221,9 @@ public:
     //
     //  \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] 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[out] mergedPoints
+    //  \param[out] mergedFaces
+    //  \param[out] pointMergeMap
     //  \param[in] useLocal gather/merge patch localFaces/localPoints
     //      instead of faces/points
     //
@@ -282,20 +243,19 @@ public:
         <
             typename PrimitivePatch<FaceList, PointField>::face_type
         >& mergedFaces,
-        labelList& pointMergeMap = const_cast<labelList&>(labelList::null()),
+        labelList& pointMergeMap,
         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 sub-ranks returns:
+    //  On master and slave 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 returns:
+    //  On master only:
     //  - mergedFaces : the merged faces
     //  - mergedPoints : the merged points
     template<class FaceList>

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

@@ -46,19 +46,17 @@ 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
-    faceAddr = globalIndex(pp.size(), globalIndex::gatherOnly{});
+    const globalIndex faceAddr(pp.size(), globalIndex::gatherOnly{});
 
     // Points from all ranks
-    pointAddr = globalIndex
+    const globalIndex pointAddr
     (
         (useLocal ? pp.localPoints().size() : pp.points().size()),
         globalIndex::gatherOnly{}
@@ -154,40 +152,6 @@ 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/parallel/globalIndex/globalIndex.H

@@ -108,8 +108,8 @@ public:
 
     // Constructors
 
-        //- Default construct (empty)
-        globalIndex() noexcept = default;
+        //- Default construct
+        globalIndex() = default;
 
         //- Copy construct from a list of offsets.
         //- No communication required
@@ -184,15 +184,9 @@ 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
 
@@ -208,8 +202,8 @@ public:
 
     // Edit
 
-        //- Reset to be empty (no offsets)
-        inline void clear();
+        //- Write-access to the offsets, for changing after construction
+        inline labelList& offsets() noexcept;
 
         //- Reset from local size, using gather/broadcast
         //- with default/specified communicator if parallel.

src/OpenFOAM/parallel/globalIndex/globalIndexI.H

@@ -176,12 +176,6 @@ 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,24 +40,6 @@ 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,9 +103,6 @@ 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
@@ -175,9 +172,6 @@ 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/Pstream/dummy/UOPstreamWrite.C

@@ -41,7 +41,7 @@ bool Foam::UOPstream::bufferIPCsend()
 
 bool Foam::UOPstream::write
 (
-    const UPstream::commsTypes commsType,
+    const 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 UPstream::commsTypes commsType,
+    const 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)
-            << " commType:" << UPstream::commsTypeNames[commsType]
+            << " commsType:" << 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)
-            << " commType:" << UPstream::commsTypeNames[commsType]
+            << " commsType:" << 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)
-                << " commType:" << UPstream::commsTypeNames[commsType]
+                << " commsType:" << UPstream::commsTypeNames[commsType]
                 << " request:" << PstreamGlobals::outstandingRequests_.size()
                 << Foam::endl;
         }

src/TurbulenceModels/incompressible/turbulentTransportModels/turbulentTransportModels.C

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

src/TurbulenceModels/schemes/DEShybrid/DEShybrid.H

@@ -7,7 +7,6 @@
 -------------------------------------------------------------------------------
     Copyright (C) 2015 OpenFOAM Foundation
     Copyright (C) 2015-2022 OpenCFD Ltd.
-    Copyright (C) 2022 Upstream CFD GmbH
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -29,8 +28,7 @@ Class
     Foam::DEShybrid
 
 Description
-    Improved hybrid convection scheme of Travin et al. for hybrid RAS/LES
-    calculations with enhanced Grey Area Mitigation (GAM) behaviour.
+    Hybrid convection scheme of Travin et al. for hybrid RAS/LES calculations.
 
     The scheme provides a blend between two convection schemes, based on local
     properties including the wall distance, velocity gradient and eddy
@@ -55,26 +53,25 @@ Description
 
     First published in:
     \verbatim
-        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
+        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
     \endverbatim
 
-    Original publication contained a typo for \c C_H3 constant.
-    Corrected version with minor changes for 2 lower limiters published in:
+    Original publication contained a typo for C_H3 constant. Corrected version
+    with minor changes for 2 lower limiters published in:
     \verbatim
-        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
+        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
     \endverbatim
 
-    Example of the \c DEShybrid scheme specification using \c linear
-    within the LES region and \c linearUpwind within the RAS region:
+    Example of the DEShybrid scheme specification using linear within the LES
+    region and linearUpwind within the RAS region:
     \verbatim
     divSchemes
     {
@@ -83,14 +80,13 @@ Description
         div(phi,U)      Gauss DEShybrid
             linear                    // scheme 1
             linearUpwind grad(U)      // scheme 2
-            delta                     // LES delta name, e.g. 'delta', 'hmax'
-            0.65                      // CDES coefficient
+            hmax                      // LES delta name, e.g. 'delta', 'hmax'
+            0.65                      // DES coefficient, typically = 0.65
             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.0;                      // nut limiter (if > 1, GAM extension is active)
+            1.0e-03;                  // Limiter of B function, typically 1e-03
         .
         .
     }
@@ -101,12 +97,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 (turbulenceModelSchemes);
+        libs ("libturbulenceModelSchemes.so");
     \endverbatim
 
 SourceFiles
@@ -114,14 +110,15 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_DEShybrid_H
-#define Foam_DEShybrid_H
+#ifndef DEShybrid_H
+#define DEShybrid_H
 
 #include "surfaceInterpolationScheme.H"
 #include "surfaceInterpolate.H"
 #include "fvcGrad.H"
 #include "blendedSchemeBase.H"
-#include "turbulenceModel.H"
+#include "turbulentTransportModel.H"
+#include "turbulentFluidThermoModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -138,9 +135,6 @@ class DEShybrid
     public surfaceInterpolationScheme<Type>,
     public blendedSchemeBase<Type>
 {
-    typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
-    typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
-
     // Private Data
 
         //- Scheme 1
@@ -152,7 +146,7 @@ class DEShybrid
         //- Name of the LES delta field
         word deltaName_;
 
-        //- DES coefficient
+        //- DES Coefficient
         scalar CDES_;
 
         //- Reference velocity scale [m/s]
@@ -170,14 +164,10 @@ 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;
@@ -188,108 +178,51 @@ 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 VolFieldType& vf,
-            const volScalarField& nut,
-            const volScalarField& nu,
+            const GeometricField<Type, fvPatchField, volMesh>& vf,
+            const volScalarField& nuEff,
             const volVectorField& U,
             const volScalarField& delta
         ) const
         {
-            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)));
+            tmp<volTensorField> gradU(fvc::grad(U));
+            const volScalarField S(sqrt(2.0)*mag(symm(gradU())));
+            const volScalarField Omega(sqrt(2.0)*mag(skew(gradU())));
             const dimensionedScalar tau0_ = L0_/U0_;
 
-            tmp<volScalarField> tB =
+            const volScalarField B
+            (
                 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_))
+            );
 
-            tmp<volScalarField> tg = tanh(pow4(tB));
+            const volScalarField K
+            (
+                max(Foam::sqrt(0.5*(sqr(S) + sqr(Omega))), 0.1/tau0_)
+            );
 
-            tmp<volScalarField> tK =
-                max(Foam::sqrt(0.5*(sqr(S) + sqr(Omega))), 0.1/tau0_);
-
-            tmp<volScalarField> tlTurb =
+            const volScalarField lTurb
+            (
                 Foam::sqrt
                 (
                     max
                     (
-                        (max(nut, min(sqr(Cs_*delta)*S, nutLim_*nut)) + nu)
-                       /(pow(0.09, 1.5)*tK),
-                        dimensionedScalar(sqr(dimLength), Zero)
+                        nuEff/(pow(0.09, 1.5)*K),
+                        dimensionedScalar("l0", sqr(dimLength), 0)
                     )
-                );
+                )
+            );
+
+            const volScalarField g(tanh(pow4(B)));
 
             const volScalarField A
             (
                 CH2_*max
                 (
                     scalar(0),
-                    CDES_*delta/max(tlTurb*tg, SMALL*L0_) - 0.5
+                    CDES_*delta/max(lTurb*g, SMALL*L0_) - 0.5
                 )
             );
 
@@ -364,8 +297,14 @@ 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)),
@@ -373,13 +312,46 @@ public:
             sigmaMin_(readScalar(is)),
             sigmaMax_(readScalar(is)),
             OmegaLim_(readScalar(is)),
-            nutLim_(readScalarOrDefault(is, scalar(1))),
             CH1_(3.0),
             CH2_(1.0),
-            CH3_(2.0),
-            Cs_(0.18)
+            CH3_(2.0)
         {
-            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);
+            }
         }
 
         //- Construct from mesh, faceFlux and Istream
@@ -406,13 +378,46 @@ public:
             sigmaMin_(readScalar(is)),
             sigmaMax_(readScalar(is)),
             OmegaLim_(readScalar(is)),
-            nutLim_(readScalarOrDefault(is, scalar(1))),
             CH1_(3.0),
             CH2_(1.0),
-            CH3_(2.0),
-            Cs_(0.18)
+            CH3_(2.0)
         {
-            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: " << 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);
+            }
         }
 
 
@@ -426,24 +431,33 @@ public:
         {
             const fvMesh& mesh = this->mesh();
 
-            // Retrieve LES delta from the mesh database
-            const auto& delta =
-                mesh.lookupObject<const volScalarField>(deltaName_);
+            typedef compressible::turbulenceModel cmpModel;
+            typedef incompressible::turbulenceModel icoModel;
 
-            // Retrieve turbulence model from the mesh database
-            const auto* modelPtr =
-                mesh.cfindObject<turbulenceModel>
-                (
-                    turbulenceModel::propertiesName
-                );
+            // Lookup the LES delta from the mesh database
+            const volScalarField& delta = this->mesh().template
+                 lookupObject<const volScalarField>(deltaName_);
 
-            if (modelPtr)
+            // 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))
             {
-                const auto& model = *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);
 
                 return calcBlendingFactor
                 (
-                    vf, model.nut(), model.nu(), model.U(), delta
+                    vf, model.muEff()/model.rho(), model.U(), delta
                 );
             }
 
@@ -457,9 +471,12 @@ public:
 
 
         //- Return the interpolation weighting factors
-        tmp<surfaceScalarField> weights(const VolFieldType& vf) const
+        tmp<surfaceScalarField> weights
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& vf
+        ) const
         {
-            const surfaceScalarField bf(blendingFactor(vf));
+            surfaceScalarField bf(blendingFactor(vf));
 
             return
                 (scalar(1) - bf)*tScheme1_().weights(vf)
@@ -468,10 +485,14 @@ public:
 
 
         //- Return the face-interpolate of the given cell field
-        //- with explicit correction
-        tmp<SurfaceFieldType> interpolate(const VolFieldType& vf) const
+        //  with explicit correction
+        tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
+        interpolate
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& vf
+        ) const
         {
-            const surfaceScalarField bf(blendingFactor(vf));
+            surfaceScalarField bf(blendingFactor(vf));
 
             return
                 (scalar(1) - bf)*tScheme1_().interpolate(vf)
@@ -487,10 +508,14 @@ public:
 
 
         //- Return the explicit correction to the face-interpolate
-        //- for the given field
-        virtual tmp<SurfaceFieldType> correction(const VolFieldType& vf) const
+        //  for the given field
+        virtual tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
+        correction
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& vf
+        ) const
         {
-            const surfaceScalarField bf(blendingFactor(vf));
+            surfaceScalarField bf(blendingFactor(vf));
 
             if (tScheme1_().corrected())
             {

src/TurbulenceModels/schemes/Make/options

@@ -1,9 +1,19 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/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
 
 LIB_LIBS = \
     -lfiniteVolume \
     -lmeshTools \
-    -lturbulenceModels
+    -lcompressibleTransportModels \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
+    -lincompressibleTransportModels \
+    -lcompressibleTurbulenceModels \
+    -lfluidThermophysicalModels

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

@@ -1,499 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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

@@ -1,246 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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,8 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2022 Upstream CFD GmbH
-    Copyright (C) 2016-2022 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -48,7 +47,7 @@ tmp<volScalarField> kOmegaSSTBase<BasicEddyViscosityModel>::F1
     tmp<volScalarField> CDkOmegaPlus = max
     (
         CDkOmega,
-        dimensionedScalar(dimless/sqr(dimTime), 1.0e-10)
+        dimensionedScalar("1.0e-10", dimless/sqr(dimTime), 1.0e-10)
     );
 
     tmp<volScalarField> arg1 = min
@@ -133,17 +132,6 @@ 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
 (
@@ -155,32 +143,17 @@ tmp<volScalarField::Internal> kOmegaSSTBase<BasicEddyViscosityModel>::Pk
 
 
 template<class BasicEddyViscosityModel>
-tmp<volScalarField::Internal> kOmegaSSTBase<BasicEddyViscosityModel>::epsilonByk
+tmp<volScalarField::Internal>
+kOmegaSSTBase<BasicEddyViscosityModel>::epsilonByk
 (
-    const volScalarField& /* F1 not used */,
-    const volTensorField& /* gradU not used */
+    const volScalarField& F1,
+    const volTensorField& gradU
 ) 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
 (
@@ -192,7 +165,8 @@ 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))
     );
 }
 
@@ -200,10 +174,13 @@ tmp<volScalarField::Internal> kOmegaSSTBase<BasicEddyViscosityModel>::GbyNu
 template<class BasicEddyViscosityModel>
 tmp<fvScalarMatrix> kOmegaSSTBase<BasicEddyViscosityModel>::kSource() const
 {
-    return tmp<fvScalarMatrix>::New
+    return tmp<fvScalarMatrix>
     (
-        k_,
-        dimVolume*this->rho_.dimensions()*k_.dimensions()/dimTime
+        new fvScalarMatrix
+        (
+            k_,
+            dimVolume*this->rho_.dimensions()*k_.dimensions()/dimTime
+        )
     );
 }
 
@@ -211,10 +188,13 @@ tmp<fvScalarMatrix> kOmegaSSTBase<BasicEddyViscosityModel>::kSource() const
 template<class BasicEddyViscosityModel>
 tmp<fvScalarMatrix> kOmegaSSTBase<BasicEddyViscosityModel>::omegaSource() const
 {
-    return tmp<fvScalarMatrix>::New
+    return tmp<fvScalarMatrix>
     (
-        omega_,
-        dimVolume*this->rho_.dimensions()*omega_.dimensions()/dimTime
+        new fvScalarMatrix
+        (
+            omega_,
+            dimVolume*this->rho_.dimensions()*omega_.dimensions()/dimTime
+        )
     );
 }
 
@@ -227,10 +207,13 @@ tmp<fvScalarMatrix> kOmegaSSTBase<BasicEddyViscosityModel>::Qsas
     const volScalarField::Internal& beta
 ) const
 {
-    return tmp<fvScalarMatrix>::New
+    return tmp<fvScalarMatrix>
     (
-        omega_,
-        dimVolume*this->rho_.dimensions()*omega_.dimensions()/dimTime
+        new fvScalarMatrix
+        (
+            omega_,
+            dimVolume*this->rho_.dimensions()*omega_.dimensions()/dimTime
+        )
     );
 }
 
@@ -406,7 +389,6 @@ kOmegaSSTBase<BasicEddyViscosityModel>::kOmegaSSTBase
         ),
         this->mesh_
     ),
-
     decayControl_
     (
         Switch::getOrAddToDict
@@ -516,30 +498,31 @@ void kOmegaSSTBase<BasicEddyViscosityModel>::correct()
 
     BasicEddyViscosityModel::correct();
 
-    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());
+    volScalarField::Internal divU(fvc::div(fvc::absolute(this->phi(), U)));
 
     tmp<volTensorField> tgradU = fvc::grad(U);
-    const volScalarField S2(this->S2(F1, tgradU()));
-    volScalarField::Internal GbyNu0(this->GbyNu0(tgradU(), F1, S2));
+    volScalarField S2(2*magSqr(symm(tgradU())));
+    volScalarField::Internal GbyNu0
+    (
+        this->type() + ":GbyNu",
+        (tgradU() && dev(twoSymm(tgradU())))
+    );
     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());
+
     {
-        const volScalarField::Internal gamma(this->gamma(F1));
-        const volScalarField::Internal beta(this->beta(F1));
+        volScalarField::Internal gamma(this->gamma(F1));
+        volScalarField::Internal beta(this->beta(F1));
 
         GbyNu0 = GbyNu(GbyNu0, F23(), S2());
 
@@ -572,30 +555,28 @@ 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,8 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2022 Upstream CFD GmbH
-    Copyright (C) 2017-2022 OpenCFD Ltd.
+    Copyright (C) 2017-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -142,7 +141,7 @@ class kOmegaSSTBase
 
 protected:
 
-    // Protected Data
+    // Protected data
 
         // Model coefficients
 
@@ -175,10 +174,7 @@ 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_;
 
 
@@ -192,7 +188,6 @@ protected:
 
     // Protected Member Functions
 
-        //- Set decay control with kInf and omegaInf
         void setDecayControl(const dictionary& dict);
 
         virtual tmp<volScalarField> F1(const volScalarField& CDkOmega) const;
@@ -200,7 +195,6 @@ protected:
         virtual tmp<volScalarField> F3() const;
         virtual tmp<volScalarField> F23() const;
 
-        //- Return the blended field
         tmp<volScalarField> blend
         (
             const volScalarField& F1,
@@ -211,7 +205,6 @@ protected:
             return F1*(psi1 - psi2) + psi2;
         }
 
-        //- Return the internal blended field
         tmp<volScalarField::Internal> blend
         (
             const volScalarField::Internal& F1,
@@ -260,13 +253,6 @@ 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
         (
@@ -280,14 +266,6 @@ 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
         (
@@ -343,20 +321,22 @@ public:
         //- Return the effective diffusivity for k
         tmp<volScalarField> DkEff(const volScalarField& F1) const
         {
-            return tmp<volScalarField>::New
+            return tmp<volScalarField>
             (
-                "DkEff",
-                alphaK(F1)*this->nut_ + this->nu()
+                new volScalarField("DkEff", alphaK(F1)*this->nut_ + this->nu())
             );
         }
 
         //- Return the effective diffusivity for omega
         tmp<volScalarField> DomegaEff(const volScalarField& F1) const
         {
-            return tmp<volScalarField>::New
+            return tmp<volScalarField>
             (
-                "DomegaEff",
-                alphaOmega(F1)*this->nut_ + this->nu()
+                new volScalarField
+                (
+                    "DomegaEff",
+                    alphaOmega(F1)*this->nut_ + this->nu()
+                )
             );
         }
 

src/TurbulenceModels/turbulenceModels/DES/DESModel/DESModel.C

@@ -6,8 +6,6 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2015 OpenFOAM Foundation
-    Copyright (C) 2022 Upstream CFD GmbH
-    Copyright (C) 2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -61,127 +59,15 @@ DESModel<BasicTurbulenceModel>::DESModel
         phi,
         transport,
         propertiesName
-    ),
-    Ctrans_(dimless, Zero)
-{
-    // Note: Ctrans is model-specific and initialised in derived classes
-}
+    )
+{}
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 template<class BasicTurbulenceModel>
-bool DESModel<BasicTurbulenceModel>::read()
-{
-    if (LESeddyViscosity<BasicTurbulenceModel>::read())
-    {
-        Ctrans_.readIfPresent(this->coeffDict());
-
-        return true;
-    }
-
-    return false;
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> DESModel<BasicTurbulenceModel>::Ssigma
-(
-    const volTensorField& gradU,
-    const dimensionedScalar& coeff
-)
-{
-    // Limiter
-    const dimensionedScalar eps0(dimless, SMALL);
-    const dimensionedScalar eps2(dimless/sqr(dimTime), SMALL);
-    const dimensionedScalar eps4(dimless/pow4(dimTime), SMALL);
-    const dimensionedScalar max2(dimless/sqr(dimTime), GREAT);
-    const dimensionedTensor maxTen2
-    (
-        dimless/sqr(dimTime),
-        tensor::max
-    );
-    const dimensionedTensor minTen2
-    (
-        dimless/sqr(dimTime),
-        tensor::min
-    );
-
-    const volTensorField G(max(min(gradU.T() & gradU, maxTen2), minTen2));
-
-    // Tensor invariants
-    const volScalarField I1(tr(G));
-    const volScalarField I2(0.5*(sqr(I1) - tr(G & G)));
-    tmp<volScalarField> tI3 = det(G);
-
-    const volScalarField alpha1(max(sqr(I1)/9.0 - I2/3.0, eps4));
-
-    tmp<volScalarField> talpha2 =
-        pow3(min(I1, max2))/27.0 - I1*I2/6.0 + 0.5*tI3;
-
-    const volScalarField alpha3
-    (
-        1.0/3.0
-       *acos
-        (
-            max
-            (
-                scalar(-1) + eps0,
-                min(scalar(1) - eps0, talpha2/pow(alpha1, 3.0/2.0))
-            )
-        )
-    );
-
-    const scalar piBy3 = constant::mathematical::pi/3.0;
-    const volScalarField sigma1
-    (
-        sqrt(max(I1/3.0 + 2.0*sqrt(alpha1)*cos(alpha3), eps2))
-    );
-    const volScalarField sigma2
-    (
-        sqrt(max(I1/3.0 - 2.0*sqrt(alpha1)*cos(piBy3 + alpha3), eps2))
-    );
-    const volScalarField sigma3
-    (
-        sqrt(max(I1/3.0 - 2.0*sqrt(alpha1)*cos(piBy3 - alpha3), eps2))
-    );
-
-    return
-        coeff
-       *sigma3
-       *(sigma1 - sigma2)
-       *(sigma2 - sigma3)
-       /max(sqr(sigma1), eps2);
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> DESModel<BasicTurbulenceModel>::Ssigma
-(
-    const volTensorField& gradU
-) const
-{
-    return Ssigma(gradU, Ctrans_);
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> DESModel<BasicTurbulenceModel>::fd() const
-{
-    return tmp<volScalarField>::New
-    (
-        IOobject
-        (
-            "fd",
-            this->mesh_.time().timeName(),
-            this->mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        this->mesh_,
-        dimensionedScalar(dimless, Zero)
-    );
-}
+DESModel<BasicTurbulenceModel>::~DESModel()
+{}
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/TurbulenceModels/turbulenceModels/DES/DESModel/DESModel.H

@@ -6,8 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2015 OpenFOAM Foundation
-    Copyright (C) 2022 Upstream CFD GmbH
-    Copyright (C) 2016, 2022 OpenCFD Ltd.
+    Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -58,6 +57,9 @@ class DESModel
     public DESModelBase,
     public LESeddyViscosity<BasicTurbulenceModel>
 {
+
+private:
+
     // Private Member Functions
 
         //- No copy construct
@@ -67,14 +69,6 @@ class DESModel
         void operator=(const DESModel&) = delete;
 
 
-protected:
-
-    // Protected Data
-
-        //- Model-specific transition constant
-        dimensionedScalar Ctrans_;
-
-
 public:
 
     typedef typename BasicTurbulenceModel::alphaField alphaField;
@@ -98,30 +92,13 @@ public:
 
 
     //- Destructor
-    virtual ~DESModel() = default;
+    virtual ~DESModel();
 
 
     // Public Member Functions
 
-        //- Re-read model coefficients if they have changed
-        virtual bool read();
-
         //- Return the LES field indicator
         virtual tmp<volScalarField> LESRegion() const = 0;
-
-        //- Return modified strain rate
-        static tmp<volScalarField> Ssigma
-        (
-            const volTensorField& gradU,
-            const dimensionedScalar& coeff
-        );
-
-        //- Return modified strain rate
-        //  Note: uses Ctrans_ coefficient
-        virtual tmp<volScalarField> Ssigma(const volTensorField& gradU) const;
-
-        //- Return the shielding function
-        virtual tmp<volScalarField> fd() const;
 };
 
 

src/TurbulenceModels/turbulenceModels/DES/DESModel/DESModelBase.C

@@ -23,6 +23,16 @@ License
     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::DESModelBase
+
+Description
+    Base class for DES models providing an interfaces to the LESRegion
+    function.
+
+SourceFiles
+    DESModelBase.C
+
 \*---------------------------------------------------------------------------*/
 
 #include "DESModelBase.H"

src/TurbulenceModels/turbulenceModels/DES/DESModel/DESModelBase.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015 OpenFOAM Foundation
-    Copyright (C) 2019-2022 OpenCFD Ltd.
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,7 +28,8 @@ Class
     Foam::DESModelBase
 
 Description
-    Base class for DES models providing an interfaces to DES fields.
+    Base class for DES models providing an interfaces to the LESRegion
+    function.
 
 SourceFiles
     DESModelBase.C
@@ -69,9 +70,6 @@ public:
 
         //- Return the LES field indicator
         virtual tmp<volScalarField> LESRegion() const = 0;
-
-        //- Return the shielding function
-        virtual tmp<volScalarField> fd() const = 0;
 };
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/TurbulenceModels/turbulenceModels/DES/SpalartAllmarasDDES/SpalartAllmarasDDES.C

@@ -6,8 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2022 Upstream CFD GmbH
-    Copyright (C) 2015-2022 OpenCFD Ltd.
+    Copyright (C) 2015-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -38,60 +37,46 @@ namespace LESModels
 
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
+template<class BasicTurbulenceModel>
+tmp<volScalarField> SpalartAllmarasDDES<BasicTurbulenceModel>::rd
+(
+    const volScalarField& magGradU
+) const
+{
+    tmp<volScalarField> tr
+    (
+        min
+        (
+            this->nuEff()
+           /(
+                max
+                (
+                    magGradU,
+                    dimensionedScalar("SMALL", magGradU.dimensions(), SMALL)
+                )
+               *sqr(this->kappa_*this->y_)
+            ),
+            scalar(10)
+        )
+    );
+    tr.ref().boundaryFieldRef() == 0.0;
+
+    return tr;
+}
+
+
 template<class BasicTurbulenceModel>
 tmp<volScalarField> SpalartAllmarasDDES<BasicTurbulenceModel>::fd
 (
     const volScalarField& magGradU
 ) const
 {
-    return
-        1
-      - tanh(pow(this->Cd1_*this->r(this->nuEff(), magGradU, this->y_), Cd2_));
+    return 1 - tanh(pow(Cd1_*rd(magGradU), Cd2_));
 }
 
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmarasDDES<BasicTurbulenceModel>::Stilda
-(
-    const volScalarField& chi,
-    const volScalarField& fv1,
-    const volTensorField& gradU,
-    const volScalarField& dTilda
-) const
-{
-    if (this->useSigma_)
-    {
-        const volScalarField& lRAS(this->y_);
-        const volScalarField fv2(this->fv2(chi, fv1));
-        const volScalarField lLES(this->lengthScaleLES(chi, fv1));
-        const volScalarField Omega(this->Omega(gradU));
-        const volScalarField Ssigma(this->Ssigma(gradU));
-        const volScalarField SsigmaDES
-        (
-            Omega - fd(mag(gradU))*pos(lRAS - lLES)*(Omega - Ssigma)
-        );
-
-        return
-            max
-            (
-                SsigmaDES + fv2*this->nuTilda_/sqr(this->kappa_*dTilda),
-                this->Cs_*SsigmaDES
-            );
-    }
-
-    return
-        SpalartAllmarasBase<DESModel<BasicTurbulenceModel>>::Stilda
-        (
-            chi,
-            fv1,
-            gradU,
-            dTilda
-        );
-}
-
-
 template<class BasicTurbulenceModel>
 tmp<volScalarField> SpalartAllmarasDDES<BasicTurbulenceModel>::dTilda
 (
@@ -101,12 +86,17 @@ tmp<volScalarField> SpalartAllmarasDDES<BasicTurbulenceModel>::dTilda
 ) const
 {
     const volScalarField& lRAS(this->y_);
-    const volScalarField lLES(this->lengthScaleLES(chi, fv1));
-    const dimensionedScalar l0(dimLength, Zero);
+    const volScalarField lLES(this->psi(chi, fv1)*this->CDES_*this->delta());
 
     return max
     (
-        lRAS - fd(mag(gradU))*max(lRAS - lLES, l0),
+        lRAS
+      - fd(mag(gradU))
+       *max
+        (
+            lRAS - lLES,
+            dimensionedScalar(dimLength, Zero)
+        ),
         dimensionedScalar("small", dimLength, SMALL)
     );
 }
@@ -141,19 +131,12 @@ SpalartAllmarasDDES<BasicTurbulenceModel>::SpalartAllmarasDDES
 
     Cd1_
     (
-        this->useSigma_
-          ? dimensioned<scalar>::getOrAddToDict
-            (
-                "Cd1Sigma",
-                this->coeffDict_,
-                10
-            )
-          : dimensioned<scalar>::getOrAddToDict
-            (
-                "Cd1",
-                this->coeffDict_,
-                8
-            )
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Cd1",
+            this->coeffDict_,
+            8
+        )
     ),
     Cd2_
     (
@@ -189,13 +172,6 @@ bool SpalartAllmarasDDES<BasicTurbulenceModel>::read()
 }
 
 
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmarasDDES<BasicTurbulenceModel>::fd() const
-{
-    return fd(mag(fvc::grad(this->U_)));
-}
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace LESModels

src/TurbulenceModels/turbulenceModels/DES/SpalartAllmarasDDES/SpalartAllmarasDDES.H

@@ -6,8 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2022 Upstream CFD GmbH
-    Copyright (C) 2019-2022 OpenCFD Ltd.
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -33,16 +32,15 @@ Group
 
 Description
     SpalartAllmaras DDES turbulence model for incompressible and compressible
-    flows.
+    flows
 
     Reference:
     \verbatim
-        Spalart, P. R., Deck, S., Shur, M. L., Squires,
-        K. D., Strelets, M. K., & Travin, A. (2006).
-        A new version of detached-eddy simulation,
-        resistant to ambiguous grid densities.
+        Spalart, P. R., Deck, S., Shur, M. L., Squires, K. D., Strelets, M. K.,
+        & Travin, A. (2006).
+        A new version of detached-eddy simulation, resistant to ambiguous grid
+        densities.
         Theoretical and computational fluid dynamics, 20(3), 181-195.
-        DOI:10.1007/s00162-006-0015-0
     \endverbatim
 
 SourceFiles
@@ -50,10 +48,10 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_SpalartAllmarasDDES_H
-#define Foam_SpalartAllmarasDDES_H
+#ifndef SpalartAllmarasDDES_H
+#define SpalartAllmarasDDES_H
 
-#include "SpalartAllmarasBase.H"
+#include "SpalartAllmarasDES.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -73,9 +71,10 @@ class SpalartAllmarasDDES
 {
     // Private Member Functions
 
-        //- Return the shielding function
         tmp<volScalarField> fd(const volScalarField& magGradU) const;
 
+        tmp<volScalarField> rd(const volScalarField& magGradU) const;
+
         //- No copy construct
         SpalartAllmarasDDES(const SpalartAllmarasDDES&) = delete;
 
@@ -85,7 +84,7 @@ class SpalartAllmarasDDES
 
 protected:
 
-    // Protected Data
+    // Protected data
 
         // Model coefficients
 
@@ -95,16 +94,7 @@ protected:
 
     // Protected Member Functions
 
-        //- Return the production term
-        virtual tmp<volScalarField> Stilda
-        (
-            const volScalarField& chi,
-            const volScalarField& fv1,
-            const volTensorField& gradU,
-            const volScalarField& dTilda
-        ) const;
-
-        //- Return the length scale
+        //- Length scale
         virtual tmp<volScalarField> dTilda
         (
             const volScalarField& chi,
@@ -148,9 +138,6 @@ public:
 
         //- Read from dictionary
         virtual bool read();
-
-        //- Return the shielding function
-        virtual tmp<volScalarField> fd() const;
 };
 
 

src/TurbulenceModels/turbulenceModels/DES/SpalartAllmarasDES/SpalartAllmarasDES.C

@@ -6,8 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2022 Upstream CFD GmbH
-    Copyright (C) 2016-2022 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -39,6 +38,120 @@ namespace LESModels
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
+template<class BasicTurbulenceModel>
+tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::chi() const
+{
+    return nuTilda_/this->nu();
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::fv1
+(
+    const volScalarField& chi
+) const
+{
+    const volScalarField chi3("chi3", pow3(chi));
+    return chi3/(chi3 + pow3(Cv1_));
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::fv2
+(
+    const volScalarField& chi,
+    const volScalarField& fv1
+) const
+{
+    return 1.0 - chi/(1.0 + chi*fv1);
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::ft2
+(
+    const volScalarField& chi
+) const
+{
+    return Ct3_*exp(-Ct4_*sqr(chi));
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::Omega
+(
+    const volTensorField& gradU
+) const
+{
+    return sqrt(2.0)*mag(skew(gradU));
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::Stilda
+(
+    const volScalarField& chi,
+    const volScalarField& fv1,
+    const volScalarField& Omega,
+    const volScalarField& dTilda
+) const
+{
+    return
+    (
+        max
+        (
+            Omega
+          + fv2(chi, fv1)*nuTilda_/sqr(kappa_*dTilda),
+            Cs_*Omega
+        )
+    );
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::r
+(
+    const volScalarField& nur,
+    const volScalarField& Omega,
+    const volScalarField& dTilda
+) const
+{
+    tmp<volScalarField> tr
+    (
+        min
+        (
+            nur
+           /(
+                max
+                (
+                    Omega,
+                    dimensionedScalar("SMALL", Omega.dimensions(), SMALL)
+                )
+                *sqr(kappa_*dTilda)
+            ),
+            scalar(10)
+        )
+    );
+    tr.ref().boundaryFieldRef() == 0.0;
+
+    return tr;
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::fw
+(
+    const volScalarField& Omega,
+    const volScalarField& dTilda
+) const
+{
+    const volScalarField r(this->r(nuTilda_, Omega, dTilda));
+    const volScalarField g(r + Cw2_*(pow6(r) - r));
+
+    return g*pow((1 + pow6(Cw3_))/(pow6(g) + pow6(Cw3_)), 1.0/6.0);
+}
+
+
 template<class BasicTurbulenceModel>
 tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::psi
 (
@@ -46,23 +159,26 @@ tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::psi
     const volScalarField& fv1
 ) const
 {
-    auto tpsi = tmp<volScalarField>::New
+    tmp<volScalarField> tpsi
     (
-        IOobject
+        new volScalarField
         (
-            IOobject::scopedName(type(), "psi"),
-            this->time().timeName(),
+            IOobject
+            (
+                type() + ":psi",
+                this->time().timeName(),
+                this->mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
             this->mesh(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        this->mesh(),
-        dimensionedScalar(dimless, Zero)
+            dimensionedScalar("one", dimless, 1)
+        )
     );
 
     if (lowReCorrection_)
     {
-        auto& psi = tpsi.ref();
+        volScalarField& psi = tpsi.ref();
 
         const volScalarField fv2(this->fv2(chi, fv1));
         const volScalarField ft2(this->ft2(chi));
@@ -73,8 +189,7 @@ tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::psi
                 min
                 (
                     scalar(100),
-                    (1 - this->Cb1_/(this->Cw1_*sqr(this->kappa_)*fwStar_)
-                   *(ft2 + (1 - ft2)*fv2))
+                    (1 - Cb1_/(Cw1_*sqr(kappa_)*fwStar_)*(ft2 + (1 - ft2)*fv2))
                    /max(SMALL, (fv1*max(scalar(1e-10), 1 - ft2)))
                 )
             );
@@ -84,57 +199,6 @@ tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::psi
 }
 
 
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::lengthScaleLES
-(
-    const volScalarField& chi,
-    const volScalarField& fv1
-) const
-{
-    return psi(chi, fv1)*CDES_*this->delta();
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::Stilda
-(
-    const volScalarField& chi,
-    const volScalarField& fv1,
-    const volTensorField& gradU,
-    const volScalarField& dTilda
-) const
-{
-    if (useSigma_)
-    {
-        const volScalarField& lRAS = this->y_;
-        const volScalarField fv2(this->fv2(chi, fv1));
-        const volScalarField lLES(this->lengthScaleLES(chi, fv1));
-        const volScalarField Omega(this->Omega(gradU));
-        const volScalarField Ssigma(this->Ssigma(gradU));
-        const volScalarField SsigmaDES
-        (
-            pos(dTilda - lRAS)*Omega + (scalar(1) - pos(dTilda - lRAS))*Ssigma
-        );
-
-        return
-            max
-            (
-                SsigmaDES + fv2*this->nuTilda_/sqr(this->kappa_*dTilda),
-                this->Cs_*SsigmaDES
-            );
-    }
-
-    return
-        SpalartAllmarasBase<DESModel<BasicTurbulenceModel>>::Stilda
-        (
-            chi,
-            fv1,
-            gradU,
-            dTilda
-        );
-}
-
-
 template<class BasicTurbulenceModel>
 tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::dTilda
 (
@@ -143,15 +207,33 @@ tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::dTilda
     const volTensorField& gradU
 ) const
 {
-    // Initialise with LES length scale
-    tmp<volScalarField> tdTilda = lengthScaleLES(chi, fv1);
-
-    // Take min vs wall distance
-    min(tdTilda.ref(), tdTilda(), this->y_);
+    tmp<volScalarField> tdTilda(psi(chi, fv1)*CDES_*this->delta());
+    min(tdTilda.ref().ref(), tdTilda(), y_);
     return tdTilda;
 }
 
 
+template<class BasicTurbulenceModel>
+void SpalartAllmarasDES<BasicTurbulenceModel>::correctNut
+(
+    const volScalarField& fv1
+)
+{
+    this->nut_ = nuTilda_*fv1;
+    this->nut_.correctBoundaryConditions();
+    fv::options::New(this->mesh_).correct(this->nut_);
+
+    BasicTurbulenceModel::correctNut();
+}
+
+
+template<class BasicTurbulenceModel>
+void SpalartAllmarasDES<BasicTurbulenceModel>::correctNut()
+{
+    correctNut(fv1(this->chi()));
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class BasicTurbulenceModel>
@@ -167,7 +249,7 @@ SpalartAllmarasDES<BasicTurbulenceModel>::SpalartAllmarasDES
     const word& type
 )
 :
-    SpalartAllmarasBase<DESModel<BasicTurbulenceModel>>
+    DESModel<BasicTurbulenceModel>
     (
         type,
         alpha,
@@ -179,13 +261,77 @@ SpalartAllmarasDES<BasicTurbulenceModel>::SpalartAllmarasDES
         propertiesName
     ),
 
-    useSigma_
+    sigmaNut_
     (
-        Switch::getOrAddToDict
+        dimensioned<scalar>::getOrAddToDict
         (
-            "useSigma",
+            "sigmaNut",
             this->coeffDict_,
-            false
+            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
         )
     ),
     CDES_
@@ -197,6 +343,15 @@ SpalartAllmarasDES<BasicTurbulenceModel>::SpalartAllmarasDES
             0.65
         )
     ),
+    ck_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "ck",
+            this->coeffDict_,
+            0.07
+        )
+    ),
     lowReCorrection_
     (
         Switch::getOrAddToDict
@@ -206,6 +361,24 @@ SpalartAllmarasDES<BasicTurbulenceModel>::SpalartAllmarasDES
             true
         )
     ),
+    Ct3_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Ct3",
+            this->coeffDict_,
+            1.2
+        )
+    ),
+    Ct4_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Ct4",
+            this->coeffDict_,
+            0.5
+        )
+    ),
     fwStar_
     (
         dimensioned<scalar>::getOrAddToDict
@@ -214,19 +387,25 @@ SpalartAllmarasDES<BasicTurbulenceModel>::SpalartAllmarasDES
             this->coeffDict_,
             0.424
         )
-    )
-{
-    // Note: Ctrans coeff is model specific; for S-A = 67.7
-    this->Ctrans_ =
-        dimensioned<scalar>::getOrAddToDict("Ctrans", this->coeffDict_, 67.7);
+    ),
 
+    nuTilda_
+    (
+        IOobject
+        (
+            "nuTilda",
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh_
+    ),
+
+    y_(wallDist::New(this->mesh_).y())
+{
     if (type == typeName)
     {
-        WarningInFunction
-            << "This model is not recommended and will be deprecated in future "
-            << "releases. Please consider using the DDES/IDDES versions instead"
-            << endl;
-
         this->printCoeffs(type);
     }
 }
@@ -237,11 +416,25 @@ SpalartAllmarasDES<BasicTurbulenceModel>::SpalartAllmarasDES
 template<class BasicTurbulenceModel>
 bool SpalartAllmarasDES<BasicTurbulenceModel>::read()
 {
-    if (SpalartAllmarasBase<DESModel<BasicTurbulenceModel>>::read())
+    if (DESModel<BasicTurbulenceModel>::read())
     {
-        useSigma_.readIfPresent("useSigma", this->coeffDict());
+        sigmaNut_.readIfPresent(this->coeffDict());
+        kappa_.readIfPresent(*this);
+
+        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());
+
         CDES_.readIfPresent(this->coeffDict());
+        ck_.readIfPresent(this->coeffDict());
+
         lowReCorrection_.readIfPresent("lowReCorrection", this->coeffDict());
+        Ct3_.readIfPresent(this->coeffDict());
+        Ct4_.readIfPresent(this->coeffDict());
         fwStar_.readIfPresent(this->coeffDict());
 
         return true;
@@ -252,24 +445,124 @@ bool SpalartAllmarasDES<BasicTurbulenceModel>::read()
 
 
 template<class BasicTurbulenceModel>
-tmp<volScalarField>
-SpalartAllmarasDES<BasicTurbulenceModel>::LESRegion() const
+tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::
+DnuTildaEff() const
+{
+    return tmp<volScalarField>
+    (
+        new volScalarField("DnuTildaEff", (nuTilda_ + this->nu())/sigmaNut_)
+    );
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::k() const
 {
     const volScalarField chi(this->chi());
     const volScalarField fv1(this->fv1(chi));
 
-    return tmp<volScalarField>::New
+    tmp<volScalarField> tdTilda
     (
-        IOobject
+        new volScalarField
         (
-            IOobject::scopedName("DES", "LESRegion"),
-            this->mesh_.time().timeName(),
+            IOobject
+            (
+                "dTilda",
+                this->mesh_.time().timeName(),
+                this->mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
             this->mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        neg(dTilda(chi, fv1, fvc::grad(this->U_)) - this->y_)
+            dimensionedScalar(dimLength, Zero),
+            zeroGradientFvPatchField<vector>::typeName
+        )
     );
+    volScalarField& dTildaF = tdTilda.ref();
+    dTildaF = dTilda(chi, fv1, fvc::grad(this->U_));
+    dTildaF.correctBoundaryConditions();
+
+    return sqr(this->nut()/ck_/dTildaF);
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::LESRegion() const
+{
+    const volScalarField chi(this->chi());
+    const volScalarField fv1(this->fv1(chi));
+
+    tmp<volScalarField> tLESRegion
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "DES::LESRegion",
+                this->mesh_.time().timeName(),
+                this->mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            neg(dTilda(chi, fv1, fvc::grad(this->U_)) - y_)
+        )
+    );
+
+    return tLESRegion;
+}
+
+
+template<class BasicTurbulenceModel>
+void SpalartAllmarasDES<BasicTurbulenceModel>::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_));
+
+    DESModel<BasicTurbulenceModel>::correct();
+
+    const volScalarField chi(this->chi());
+    const volScalarField fv1(this->fv1(chi));
+    const volScalarField ft2(this->ft2(chi));
+
+    tmp<volTensorField> tgradU = fvc::grad(U);
+    const volScalarField Omega(this->Omega(tgradU()));
+    const volScalarField dTilda(this->dTilda(chi, fv1, tgradU()));
+    const volScalarField Stilda(this->Stilda(chi, fv1, Omega, dTilda));
+
+    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();
 }
 
 

src/TurbulenceModels/turbulenceModels/DES/SpalartAllmarasDES/SpalartAllmarasDES.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2022 OpenCFD Ltd.
+    Copyright (C) 2015-2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,7 +32,7 @@ Group
 
 Description
     SpalartAllmarasDES DES turbulence model for incompressible and
-    compressible flows.
+    compressible flows
 
     Reference:
     \verbatim
@@ -44,12 +44,11 @@ Description
 
     Including the low Reynolds number correction documented in
     \verbatim
-        Spalart, P. R., Deck, S., Shur, M. L., Squires,
-        K. D., Strelets, M. K., & Travin, A. (2006).
-        A new version of detached-eddy simulation,
-        resistant to ambiguous grid densities.
-        Theoretical and computational fluid dynamics, 20(3), 181-195.
-        DOI:10.1007/s00162-006-0015-0
+        Spalart, P. R., Deck, S., Shur, M.L., Squires, K.D., Strelets, M.Kh,
+        Travin, A. (2006).
+        A new version of detached-eddy simulation, resistant to ambiguous grid
+        densities.
+        Theor. Comput. Fluid Dyn., 20, 181-195.
     \endverbatim
 
 Note
@@ -62,10 +61,9 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_SpalartAllmarasDES_H
-#define Foam_SpalartAllmarasDES_H
+#ifndef SpalartAllmarasDES_H
+#define SpalartAllmarasDES_H
 
-#include "SpalartAllmarasBase.H"
 #include "DESModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -82,7 +80,7 @@ namespace LESModels
 template<class BasicTurbulenceModel>
 class SpalartAllmarasDES
 :
-    public SpalartAllmarasBase<DESModel<BasicTurbulenceModel>>
+    public DESModel<BasicTurbulenceModel>
 {
     // Private Member Functions
 
@@ -95,47 +93,86 @@ class SpalartAllmarasDES
 
 protected:
 
-    // Protected Data
-
-        //- Switch to activate grey-area enhanced sigma-(D)DES
-        Switch useSigma_;
+    // Protected data
 
         // Model constants
 
-            //- DES coefficient
-            dimensionedScalar CDES_;
+            dimensionedScalar sigmaNut_;
+            dimensionedScalar kappa_;
 
-            //- Flag for low Reynolds number correction
-            Switch lowReCorrection_;
-            dimensionedScalar fwStar_;
+            dimensionedScalar Cb1_;
+            dimensionedScalar Cb2_;
+            dimensionedScalar Cw1_;
+            dimensionedScalar Cw2_;
+            dimensionedScalar Cw3_;
+            dimensionedScalar Cv1_;
+            dimensionedScalar Cs_;
+            dimensionedScalar CDES_;
+            dimensionedScalar ck_;
+
+
+            // Low Reynolds number correction
+
+                Switch lowReCorrection_;
+                dimensionedScalar Ct3_;
+                dimensionedScalar Ct4_;
+                dimensionedScalar fwStar_;
+
+
+        // Fields
+
+            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
 
-        //- Return the low Reynolds number correction function
-        virtual tmp<volScalarField> psi
+        tmp<volScalarField> chi() const;
+
+        tmp<volScalarField> fv1(const volScalarField& chi) const;
+
+        tmp<volScalarField> fv2
         (
             const volScalarField& chi,
             const volScalarField& fv1
         ) const;
 
-        //- Return the LES length scale
-        virtual tmp<volScalarField> lengthScaleLES
-        (
-            const volScalarField& chi,
-            const volScalarField& fv1
-        ) const;
+        tmp<volScalarField> ft2(const volScalarField& chi) const;
 
-        //- Return the production term
-        virtual tmp<volScalarField> Stilda
+        tmp<volScalarField> Omega(const volTensorField& gradU) const;
+
+        tmp<volScalarField> Stilda
         (
             const volScalarField& chi,
             const volScalarField& fv1,
-            const volTensorField& gradU,
+            const volScalarField& Omega,
             const volScalarField& dTilda
         ) const;
 
-        //- Return the length scale
+        tmp<volScalarField> r
+        (
+            const volScalarField& nur,
+            const volScalarField& Omega,
+            const volScalarField& dTilda
+        ) const;
+
+        tmp<volScalarField> fw
+        (
+            const volScalarField& Omega,
+            const volScalarField& dTilda
+        ) const;
+
+        tmp<volScalarField> psi
+        (
+            const volScalarField& chi,
+            const volScalarField& fv1
+        ) const;
+
+        //- Length scale
         virtual tmp<volScalarField> dTilda
         (
             const volScalarField& chi,
@@ -143,6 +180,9 @@ protected:
             const volTensorField& gradU
         ) const;
 
+        void correctNut(const volScalarField& fv1);
+        virtual void correctNut();
+
 
 public:
 
@@ -180,8 +220,22 @@ public:
         //- Re-read model coefficients if they have changed
         virtual bool read();
 
+        //- Return the effective diffusivity for nuTilda
+        tmp<volScalarField> DnuTildaEff() const;
+
+        //- Return SGS kinetic energy
+        virtual tmp<volScalarField> k() const;
+
+        tmp<volScalarField> nuTilda() const
+        {
+            return nuTilda_;
+        }
+
         //- Return the LES field indicator
         virtual tmp<volScalarField> LESRegion() const;
+
+        //- Correct nuTilda and related properties
+        virtual void correct();
 };
 
 

src/TurbulenceModels/turbulenceModels/DES/SpalartAllmarasIDDES/SpalartAllmarasIDDES.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2015-2022 OpenCFD Ltd.
+    Copyright (C) 2015-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -64,7 +64,7 @@ tmp<volScalarField> SpalartAllmarasIDDES<BasicTurbulenceModel>::ft
     const volScalarField& magGradU
 ) const
 {
-    return tanh(pow3(sqr(Ct_)*this->r(this->nut_, magGradU, this->y_)));
+    return tanh(pow3(sqr(Ct_)*rd(this->nut_, magGradU)));
 }
 
 
@@ -74,7 +74,36 @@ tmp<volScalarField> SpalartAllmarasIDDES<BasicTurbulenceModel>::fl
     const volScalarField& magGradU
 ) const
 {
-    return tanh(pow(sqr(Cl_)*this->r(this->nu(), magGradU, this->y_), 10));
+    return tanh(pow(sqr(Cl_)*rd(this->nu(), magGradU), 10));
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> SpalartAllmarasIDDES<BasicTurbulenceModel>::rd
+(
+    const volScalarField& nur,
+    const volScalarField& magGradU
+) const
+{
+    tmp<volScalarField> tr
+    (
+        min
+        (
+            nur
+           /(
+                max
+                (
+                    magGradU,
+                    dimensionedScalar("SMALL", magGradU.dimensions(), SMALL)
+                )
+               *sqr(this->kappa_*this->y_)
+            ),
+            scalar(10)
+        )
+    );
+    tr.ref().boundaryFieldRef() == 0.0;
+
+    return tr;
 }
 
 
@@ -84,7 +113,7 @@ tmp<volScalarField> SpalartAllmarasIDDES<BasicTurbulenceModel>::fdt
     const volScalarField& magGradU
 ) const
 {
-    return 1 - tanh(pow(Cdt1_*this->r(this->nut_, magGradU, this->y_), Cdt2_));
+    return 1 - tanh(pow(Cdt1_*rd(this->nut_, magGradU), Cdt2_));
 }
 
 
@@ -101,35 +130,31 @@ tmp<volScalarField> SpalartAllmarasIDDES<BasicTurbulenceModel>::dTilda
     const volScalarField magGradU(mag(gradU));
     const volScalarField psi(this->psi(chi, fv1));
 
-    const volScalarField& lRAS = this->y_;
+    const volScalarField& lRAS(this->y_);
     const volScalarField lLES(psi*this->CDES_*this->delta());
 
     const volScalarField alpha(this->alpha());
     const volScalarField expTerm(exp(sqr(alpha)));
 
     tmp<volScalarField> fB = min(2*pow(expTerm, -9.0), scalar(1));
+    tmp<volScalarField> fe1 =
+        2*(pos0(alpha)*pow(expTerm, -11.09) + neg(alpha)*pow(expTerm, -9.0));
+    tmp<volScalarField> fe2 = 1 - max(ft(magGradU), fl(magGradU));
+    tmp<volScalarField> fe = max(fe1 - 1, scalar(0))*psi*fe2;
+
     const volScalarField fdTilda(max(1 - fdt(magGradU), fB));
 
-    if (fe_)
-    {
-        tmp<volScalarField> fe1 =
-            2*(pos0(alpha)*pow(expTerm, -11.09) + neg(alpha)*pow(expTerm, -9.));
-        tmp<volScalarField> fe2 = 1 - max(ft(magGradU), fl(magGradU));
-        tmp<volScalarField> fe = max(fe1 - 1, scalar(0))*psi*fe2;
-
-        // Original formulation from Shur et al. paper (2008)
-        return max
-        (
-            fdTilda*(1 + fe)*lRAS + (1 - fdTilda)*lLES,
-            dimensionedScalar("SMALL", dimLength, SMALL)
-        );
-    }
-
-
     // Simplified formulation from Gritskevich et al. paper (2011) where fe = 0
+    // return max
+    // (
+    //     fdTilda*lRAS + (1 - fdTilda)*lLES,
+    //     dimensionedScalar("SMALL", dimLength, SMALL)
+    // );
+
+    // Original formulation from Shur et al. paper (2008)
     return max
     (
-        fdTilda*lRAS + (1 - fdTilda)*lLES,
+        fdTilda*(1 + fe)*lRAS + (1 - fdTilda)*lLES,
         dimensionedScalar("SMALL", dimLength, SMALL)
     );
 }
@@ -198,16 +223,6 @@ SpalartAllmarasIDDES<BasicTurbulenceModel>::SpalartAllmarasIDDES
             1.63
         )
     ),
-    fe_
-    (
-        Switch::getOrAddToDict
-        (
-            "fe",
-            this->coeffDict_,
-            true
-        )
-    ),
-
     IDDESDelta_(setDelta())
 {
     if (type == typeName)
@@ -236,17 +251,6 @@ bool SpalartAllmarasIDDES<BasicTurbulenceModel>::read()
 }
 
 
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmarasIDDES<BasicTurbulenceModel>::fd() const
-{
-    const volScalarField alpha(this->alpha());
-    const volScalarField expTerm(exp(sqr(alpha)));
-
-    tmp<volScalarField> fB = min(2*pow(expTerm, -9.0), scalar(1));
-    return max(1 - fdt(mag(fvc::grad(this->U_))), fB);
-}
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace LESModels

src/TurbulenceModels/turbulenceModels/DES/SpalartAllmarasIDDES/SpalartAllmarasIDDES.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2022 OpenCFD Ltd.
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,16 +31,15 @@ Group
     grpDESTurbulence
 
 Description
-    SpalartAllmaras IDDES turbulence model
-    for incompressible and compressible flows.
+    SpalartAllmaras IDDES turbulence model for incompressible and compressible
+    flows
 
     Reference:
     \verbatim
         Shur, M. L., Spalart, P. R., Strelets, M. K., & Travin, A. K. (2008).
-        A hybrid RANS-LES approach with delayed-DES
-        and wall-modelled LES capabilities.
-        International journal of heat and fluid flow, 29(6), 1638-1649.
-        DOI:10.1016/j.ijheatfluidflow.2008.07.001
+        A hybrid RANS-LES approach with delayed-DES and wall-modelled LES
+        capabilities.
+        International Journal of Heat and Fluid Flow, 29(6), 1638-1649.
     \endverbatim
 
 SourceFiles
@@ -48,8 +47,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_SpalartAllmarasIDDES_H
-#define Foam_SpalartAllmarasIDDES_H
+#ifndef SpalartAllmarasIDDES_H
+#define SpalartAllmarasIDDES_H
 
 #include "SpalartAllmarasDES.H"
 #include "IDDESDelta.H"
@@ -62,7 +61,7 @@ namespace LESModels
 {
 
 /*---------------------------------------------------------------------------*\
-                    Class SpalartAllmarasIDDES Declaration
+                        Class SpalartAllmarasIDDES Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class BasicTurbulenceModel>
@@ -76,11 +75,15 @@ class SpalartAllmarasIDDES
         const IDDESDelta& setDelta() const;
 
         tmp<volScalarField> alpha() const;
-
         tmp<volScalarField> ft(const volScalarField& magGradU) const;
-
         tmp<volScalarField> fl(const volScalarField& magGradU) const;
 
+        tmp<volScalarField> rd
+        (
+            const volScalarField& nur,
+            const volScalarField& magGradU
+        ) const;
+
         //- Delay function
         tmp<volScalarField> fdt(const volScalarField& magGradU) const;
 
@@ -93,7 +96,7 @@ class SpalartAllmarasIDDES
 
 protected:
 
-    // Protected Data
+    // Protected data
 
         // Model coefficients
 
@@ -101,16 +104,14 @@ protected:
             dimensionedScalar Cdt2_;
             dimensionedScalar Cl_;
             dimensionedScalar Ct_;
-            Switch fe_;
 
+        // Fields
 
-        //- IDDES delta
-        const IDDESDelta& IDDESDelta_;
-
+            const IDDESDelta& IDDESDelta_;
 
     // Protected Member Functions
 
-        //- Return the length scale
+        //- Length scale
         virtual tmp<volScalarField> dTilda
         (
             const volScalarField& chi,
@@ -154,9 +155,6 @@ public:
 
         //- Re-read model coefficients if they have changed
         virtual bool read();
-
-        //- Return the shielding function
-        virtual tmp<volScalarField> fd() const;
 };
 
 

src/TurbulenceModels/turbulenceModels/DES/kOmegaSSTDDES/kOmegaSSTDDES.C

@@ -6,8 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015 OpenFOAM Foundation
-    Copyright (C) 2022 Upstream CFD GmbH
-    Copyright (C) 2016-2022 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -38,44 +37,46 @@ namespace LESModels
 
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
+template<class BasicTurbulenceModel>
+tmp<volScalarField> kOmegaSSTDDES<BasicTurbulenceModel>::rd
+(
+    const volScalarField& magGradU
+) const
+{
+    tmp<volScalarField> tr
+    (
+        min
+        (
+            this->nuEff()
+           /(
+                max
+                (
+                    magGradU,
+                    dimensionedScalar("SMALL", magGradU.dimensions(), SMALL)
+                )
+                *sqr(this->kappa_*this->y_)
+            ),
+            scalar(10)
+        )
+    );
+    tr.ref().boundaryFieldRef() == 0.0;
+
+    return tr;
+}
+
+
 template<class BasicTurbulenceModel>
 tmp<volScalarField> kOmegaSSTDDES<BasicTurbulenceModel>::fd
 (
     const volScalarField& magGradU
 ) const
 {
-    return 1 - tanh(pow(Cd1_*this->r(this->nuEff(), magGradU), Cd2_));
+    return 1 - tanh(pow(Cd1_*rd(magGradU), Cd2_));
 }
 
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
-template<class BasicTurbulenceModel>
-tmp<volScalarField> kOmegaSSTDDES<BasicTurbulenceModel>::S2
-(
-    const volScalarField& F1,
-    const volTensorField& gradU
-) const
-{
-    tmp<volScalarField> tS2 =
-        kOmegaSSTBase<DESModel<BasicTurbulenceModel>>::S2(F1, gradU);
-
-    if (this->useSigma_)
-    {
-        volScalarField& S2 = tS2.ref();
-        const volScalarField CDES(this->CDES(F1));
-        const volScalarField Ssigma(this->Ssigma(gradU));
-
-        S2 -=
-            fd(mag(gradU))
-           *pos(this->lengthScaleRAS() - this->lengthScaleLES(CDES))
-           *(S2 - sqr(Ssigma));
-    }
-
-    return tS2;
-}
-
-
 template<class BasicTurbulenceModel>
 tmp<volScalarField> kOmegaSSTDDES<BasicTurbulenceModel>::dTilda
 (
@@ -83,46 +84,25 @@ tmp<volScalarField> kOmegaSSTDDES<BasicTurbulenceModel>::dTilda
     const volScalarField& CDES
 ) const
 {
-    const volScalarField lRAS(this->lengthScaleRAS());
-    const volScalarField lLES(this->lengthScaleLES(CDES));
-    const dimensionedScalar l0(dimLength, Zero);
+    const volScalarField& k = this->k_;
+    const volScalarField& omega = this->omega_;
+
+    const volScalarField lRAS(sqrt(k)/(this->betaStar_*omega));
+    const volScalarField lLES(CDES*this->delta());
 
     return max
     (
-        lRAS - fd(magGradU)*max(lRAS - lLES, l0),
+        lRAS
+      - fd(magGradU)
+       *max
+        (
+            lRAS - lLES,
+            dimensionedScalar(dimLength, Zero)
+        ),
         dimensionedScalar("small", dimLength, SMALL)
     );
 }
 
-template<class BasicTurbulenceModel>
-tmp<volScalarField::Internal> kOmegaSSTDDES<BasicTurbulenceModel>::GbyNu0
-(
-    const volTensorField& gradU,
-    const volScalarField& F1,
-    const volScalarField& S2
-) const
-{
-    if (this->useSigma_)
-    {
-        return S2();
-    }
-
-    return
-        kOmegaSSTBase<DESModel<BasicTurbulenceModel>>::GbyNu0(gradU, F1, S2);
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField::Internal> kOmegaSSTDDES<BasicTurbulenceModel>::GbyNu
-(
-    const volScalarField::Internal& GbyNu0,
-    const volScalarField::Internal& F2,
-    const volScalarField::Internal& S2
-) const
-{
-    return GbyNu0; // Unlimited
-}
-
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -153,19 +133,12 @@ kOmegaSSTDDES<BasicTurbulenceModel>::kOmegaSSTDDES
 
     Cd1_
     (
-        this->useSigma_
-          ? dimensioned<scalar>::getOrAddToDict
-            (
-                "Cd1Sigma",
-                this->coeffDict_,
-                22
-            )
-          : dimensioned<scalar>::getOrAddToDict
-            (
-                "Cd1",
-                this->coeffDict_,
-                20
-            )
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Cd1",
+            this->coeffDict_,
+            20
+        )
     ),
     Cd2_
     (
@@ -201,13 +174,6 @@ bool kOmegaSSTDDES<BasicTurbulenceModel>::read()
 }
 
 
-template<class BasicTurbulenceModel>
-tmp<volScalarField> kOmegaSSTDDES<BasicTurbulenceModel>::fd() const
-{
-    return fd(mag(fvc::grad(this->U_)));
-}
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace LESModels

src/TurbulenceModels/turbulenceModels/DES/kOmegaSSTDDES/kOmegaSSTDDES.H

@@ -6,8 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015 OpenFOAM Foundation
-    Copyright (C) 2022 Upstream CFD GmbH
-    Copyright (C) 2019, 2022 OpenCFD Ltd.
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,16 +31,14 @@ Group
     grpDESTurbulence
 
 Description
-    k-omega-SST DDES turbulence model for incompressible and compressible flows.
+    k-omega-SST DDES turbulence model for incompressible and compressible flows
 
     Reference:
     \verbatim
-        Gritskevich, M. S., Garbaruk, A. V.,
-        Schütze, J., & Menter, F. R. (2012).
-        Development of DDES and IDDES formulations for
-        the k-ω shear stress transport model.
-        Flow, turbulence and combustion, 88(3), 431-449.
-        DOI:10.1007/s10494-011-9378-4
+        Gritskevich, M.S., Garbaruk, A.V., Schuetze, J., Menter, F.R. (2011)
+        Development of DDES and IDDES Formulations for the k-omega
+        Shear Stress Transport Model, Flow, Turbulence and Combustion,
+        pp. 1-19
     \endverbatim
 
 SourceFiles
@@ -49,8 +46,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_kOmegaSSTDDES_H
-#define Foam_kOmegaSSTDDES_H
+#ifndef kOmegaSSTDDES_H
+#define kOmegaSSTDDES_H
 
 #include "kOmegaSSTDES.H"
 
@@ -72,9 +69,10 @@ class kOmegaSSTDDES
 {
     // Private Member Functions
 
-        //- Return the shielding function field
         tmp<volScalarField> fd(const volScalarField& magGradU) const;
 
+        tmp<volScalarField> rd(const volScalarField& magGradU) const;
+
         //- No copy construct
         kOmegaSSTDDES(const kOmegaSSTDDES&) = delete;
 
@@ -84,7 +82,7 @@ class kOmegaSSTDDES
 
 protected:
 
-    // Protected Data
+    // Protected data
 
         // Model coefficients
 
@@ -94,13 +92,6 @@ protected:
 
     // Protected Member Functions
 
-        //- Return square of strain rate
-        virtual tmp<volScalarField> S2
-        (
-            const volScalarField& F1,
-            const volTensorField& gradU
-        ) const;
-
         //- Length scale
         virtual tmp<volScalarField> dTilda
         (
@@ -108,21 +99,6 @@ protected:
             const volScalarField& CDES
         ) 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
-        (
-            const volScalarField::Internal& GbyNu0,
-            const volScalarField::Internal& F2,
-            const volScalarField::Internal& S2
-        ) const;
 
 public:
 
@@ -159,9 +135,6 @@ public:
 
         //- Re-read model coefficients if they have changed
         virtual bool read();
-
-        //- Return the shielding function
-        virtual tmp<volScalarField> fd() const;
 };
 
 

src/TurbulenceModels/turbulenceModels/DES/kOmegaSSTDES/kOmegaSSTDES.C

@@ -6,8 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015 OpenFOAM Foundation
-    Copyright (C) 2022 Upstream CFD GmbH
-    Copyright (C) 2016-2022 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -56,51 +55,6 @@ void kOmegaSSTDES<BasicTurbulenceModel>::correctNut()
 }
 
 
-template<class BasicTurbulenceModel>
-tmp<volScalarField> kOmegaSSTDES<BasicTurbulenceModel>::r
-(
-    const volScalarField& nur,
-    const volScalarField& magGradU
-) const
-{
-    const dimensionedScalar eps(magGradU.dimensions(), SMALL);
-
-    tmp<volScalarField> tr =
-        min(nur/(max(magGradU, eps)*sqr(this->kappa_*this->y_)), scalar(10));
-
-    tr.ref().boundaryFieldRef() == 0;
-
-    return tr;
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> kOmegaSSTDES<BasicTurbulenceModel>::S2
-(
-    const volScalarField& F1,
-    const volTensorField& gradU
-) const
-{
-    tmp<volScalarField> tS2 =
-        kOmegaSSTBase<DESModel<BasicTurbulenceModel>>::S2(F1, gradU);
-
-    if (this->useSigma_)
-    {
-        volScalarField& S2 = tS2.ref();
-        const volScalarField CDES(this->CDES(F1));
-        const volScalarField dTilda(this->dTilda(mag(gradU), CDES));
-        const volScalarField lengthScaleRAS(this->lengthScaleRAS());
-        const volScalarField Ssigma(this->Ssigma(gradU));
-
-        S2 =
-            pos(dTilda - lengthScaleRAS)*S2
-          + (scalar(1) - pos(dTilda - lengthScaleRAS))*sqr(Ssigma);
-    }
-
-    return tS2;
-}
-
-
 template<class BasicTurbulenceModel>
 tmp<volScalarField> kOmegaSSTDES<BasicTurbulenceModel>::dTilda
 (
@@ -108,7 +62,10 @@ tmp<volScalarField> kOmegaSSTDES<BasicTurbulenceModel>::dTilda
     const volScalarField& CDES
 ) const
 {
-    return min(lengthScaleLES(CDES), lengthScaleRAS());
+    const volScalarField& k = this->k_;
+    const volScalarField& omega = this->omega_;
+
+    return min(CDES*this->delta(), sqrt(k)/(this->betaStar_*omega));
 }
 
 
@@ -124,24 +81,6 @@ tmp<volScalarField::Internal> kOmegaSSTDES<BasicTurbulenceModel>::epsilonByk
 }
 
 
-template<class BasicTurbulenceModel>
-tmp<volScalarField::Internal> kOmegaSSTDES<BasicTurbulenceModel>::GbyNu0
-(
-    const volTensorField& gradU,
-    const volScalarField& F1,
-    const volScalarField& S2
-) const
-{
-    if (this->useSigma_)
-    {
-        return S2();
-    }
-
-    return
-        kOmegaSSTBase<DESModel<BasicTurbulenceModel>>::GbyNu0(gradU, F1, S2);
-}
-
-
 template<class BasicTurbulenceModel>
 tmp<volScalarField::Internal> kOmegaSSTDES<BasicTurbulenceModel>::GbyNu
 (
@@ -181,15 +120,6 @@ kOmegaSSTDES<BasicTurbulenceModel>::kOmegaSSTDES
         propertiesName
     ),
 
-    useSigma_
-    (
-        Switch::getOrAddToDict
-        (
-            "useSigma",
-            this->coeffDict_,
-            false
-        )
-    ),
     kappa_
     (
         dimensioned<scalar>::getOrAddToDict
@@ -218,17 +148,8 @@ kOmegaSSTDES<BasicTurbulenceModel>::kOmegaSSTDES
         )
     )
 {
-    // Note: Ctrans coeff is model specific; for k-w = 60
-    this->Ctrans_ =
-        dimensioned<scalar>::getOrAddToDict("Ctrans", this->coeffDict_, 60.0);
-
     if (type == typeName)
     {
-        WarningInFunction
-            << "This model is not recommended and will be deprecated in future "
-            << "releases. Please consider using the DDES/IDDES versions instead"
-            << endl;
-
         this->printCoeffs(type);
     }
 }
@@ -241,7 +162,6 @@ bool kOmegaSSTDES<BasicTurbulenceModel>::read()
 {
     if (kOmegaSSTBase<DESModel<BasicTurbulenceModel>>::read())
     {
-        useSigma_.readIfPresent("useSigma", this->coeffDict());
         kappa_.readIfPresent(this->coeffDict());
         CDESkom_.readIfPresent(this->coeffDict());
         CDESkeps_.readIfPresent(this->coeffDict());
@@ -253,28 +173,6 @@ bool kOmegaSSTDES<BasicTurbulenceModel>::read()
 }
 
 
-template<class BasicTurbulenceModel>
-Foam::tmp<Foam::volScalarField>
-kOmegaSSTDES<BasicTurbulenceModel>::lengthScaleRAS() const
-{
-    const volScalarField& k = this->k_;
-    const volScalarField& omega = this->omega_;
-
-    return sqrt(k)/(this->betaStar_*omega);
-}
-
-
-template<class BasicTurbulenceModel>
-Foam::tmp<Foam::volScalarField>
-kOmegaSSTDES<BasicTurbulenceModel>::lengthScaleLES
-(
-    const volScalarField& CDES
-) const
-{
-    return CDES*this->delta();
-}
-
-
 template<class BasicTurbulenceModel>
 tmp<volScalarField> kOmegaSSTDES<BasicTurbulenceModel>::LESRegion() const
 {
@@ -289,11 +187,31 @@ tmp<volScalarField> kOmegaSSTDES<BasicTurbulenceModel>::LESRegion() const
 
     const volScalarField F1(this->F1(CDkOmega));
 
-    return tmp<volScalarField>::New
+    tmp<volScalarField> tLESRegion
     (
-        IOobject::scopedName("DES", "LESRegion"),
-        neg(dTilda(mag(fvc::grad(U)), CDES(F1)) - lengthScaleRAS())
+        new volScalarField
+        (
+            IOobject
+            (
+                "DES::LESRegion",
+                this->mesh_.time().timeName(),
+                this->mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            neg
+            (
+                dTilda
+                (
+                    mag(fvc::grad(U)),
+                    F1*CDESkom_ + (1 - F1)*CDESkeps_
+                )
+              - sqrt(k)/(this->betaStar_*omega)
+            )
+        )
     );
+
+    return tLESRegion;
 }
 
 

src/TurbulenceModels/turbulenceModels/DES/kOmegaSSTDES/kOmegaSSTDES.H

@@ -31,13 +31,13 @@ Group
     grpDESTurbulence
 
 Description
-    k-omega-SST DES turbulence model for incompressible and compressible flows.
+    k-omega-SST DES turbulence model for incompressible and compressible flows
 
     Reference:
     \verbatim
-        Strelets, M. (2001).
-        Detached Eddy Simulation of Massively Separated Flows.
-        39th AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV.
+        Strelets, M. (2001)
+        Detached Eddy Simulation of Massively Separated Flows,
+        39th AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV
     \endverbatim
 
 Note
@@ -51,8 +51,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_kOmegaSSTDES_H
-#define Foam_kOmegaSSTDES_H
+#ifndef kOmegaSSTDES_H
+#define kOmegaSSTDES_H
 
 #include "DESModel.H"
 #include "kOmegaSSTBase.H"
@@ -65,7 +65,7 @@ namespace LESModels
 {
 
 /*---------------------------------------------------------------------------*\
-                        Class kOmegaSSTDES Declaration
+                        class kOmegaSSTDES Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class BasicTurbulenceModel>
@@ -84,16 +84,11 @@ class kOmegaSSTDES
 
 protected:
 
-    // Protected Data
-
-        //- Switch to activate grey-area enhanced sigma-(D)DES
-        Switch useSigma_;
+    // Protected data
 
         // Model coefficients
 
             dimensionedScalar kappa_;
-
-            //- DES coefficients
             dimensionedScalar CDESkom_;
             dimensionedScalar CDESkeps_;
 
@@ -109,20 +104,7 @@ protected:
         virtual void correctNut(const volScalarField& S2);
         virtual void correctNut();
 
-        tmp<volScalarField> r
-        (
-            const volScalarField& nur,
-            const volScalarField& magGradU
-        ) const;
-
-        //- Return square of strain rate
-        virtual tmp<volScalarField> S2
-        (
-            const volScalarField& F1,
-            const volTensorField& gradU
-        ) const;
-
-        //- Return length scale
+        //- Length scale
         virtual tmp<volScalarField> dTilda
         (
             const volScalarField& magGradU,
@@ -136,14 +118,6 @@ 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
         (
@@ -189,15 +163,6 @@ public:
         //- Re-read model coefficients if they have changed
         virtual bool read();
 
-        //- RAS length scale
-        virtual tmp<volScalarField> lengthScaleRAS() const;
-
-        //- LES length scale
-        virtual tmp<volScalarField> lengthScaleLES
-        (
-            const volScalarField& CDES
-        ) const;
-
         //- Return the LES field indicator
         virtual tmp<volScalarField> LESRegion() const;
 };

src/TurbulenceModels/turbulenceModels/DES/kOmegaSSTIDDES/kOmegaSSTIDDES.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015 OpenFOAM Foundation
-    Copyright (C) 2015-2022 OpenCFD Ltd.
+    Copyright (C) 2015-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -64,7 +64,7 @@ tmp<volScalarField> kOmegaSSTIDDES<BasicTurbulenceModel>::ft
     const volScalarField& magGradU
 ) const
 {
-    return tanh(pow3(sqr(Ct_)*this->r(this->nut_, magGradU)));
+    return tanh(pow3(sqr(Ct_)*rd(this->nut_, magGradU)));
 }
 
 
@@ -74,22 +74,51 @@ tmp<volScalarField> kOmegaSSTIDDES<BasicTurbulenceModel>::fl
     const volScalarField& magGradU
 ) const
 {
-    return tanh(pow(sqr(Cl_)*this->r(this->nu(), magGradU), 10));
+    return tanh(pow(sqr(Cl_)*rd(this->nu(), magGradU), 10));
 }
 
 
+template<class BasicTurbulenceModel>
+tmp<volScalarField> kOmegaSSTIDDES<BasicTurbulenceModel>::rd
+(
+    const volScalarField& nur,
+    const volScalarField& magGradU
+) const
+{
+    tmp<volScalarField> tr
+    (
+        min
+        (
+            nur
+           /(
+                max
+                (
+                    magGradU,
+                    dimensionedScalar("SMALL", magGradU.dimensions(), SMALL)
+                )
+                *sqr(this->kappa_*this->y_)
+            ),
+            scalar(10)
+        )
+    );
+    tr.ref().boundaryFieldRef() == 0.0;
+
+    return tr;
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
 template<class BasicTurbulenceModel>
 tmp<volScalarField> kOmegaSSTIDDES<BasicTurbulenceModel>::fdt
 (
     const volScalarField& magGradU
 ) const
 {
-    return 1 - tanh(pow(Cdt1_*this->r(this->nut_, magGradU), Cdt2_));
+    return 1 - tanh(pow(Cdt1_*rd(this->nut_, magGradU), Cdt2_));
 }
 
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
 template<class BasicTurbulenceModel>
 tmp<volScalarField> kOmegaSSTIDDES<BasicTurbulenceModel>::dTilda
 (
@@ -97,35 +126,35 @@ tmp<volScalarField> kOmegaSSTIDDES<BasicTurbulenceModel>::dTilda
     const volScalarField& CDES
 ) const
 {
-    const volScalarField lRAS(this->lengthScaleRAS());
-    const volScalarField lLES(this->lengthScaleLES(CDES));
+    const volScalarField& k = this->k_;
+    const volScalarField& omega = this->omega_;
+
+    const volScalarField lRAS(sqrt(k)/(this->betaStar_*omega));
+    const volScalarField lLES(CDES*this->delta());
 
     const volScalarField alpha(this->alpha());
     const volScalarField expTerm(exp(sqr(alpha)));
 
     tmp<volScalarField> fB = min(2*pow(expTerm, -9.0), scalar(1));
+    tmp<volScalarField> fe1 =
+        2*(pos0(alpha)*pow(expTerm, -11.09) + neg(alpha)*pow(expTerm, -9.0));
+    tmp<volScalarField> fe2 = 1 - max(ft(magGradU), fl(magGradU));
+    tmp<volScalarField> fe = max(fe1 - 1, scalar(0))*fe2;
+
     const volScalarField fdTilda(max(1 - fdt(magGradU), fB));
 
-    if (fe_)
-    {
-        tmp<volScalarField> fe1 =
-            2*(pos0(alpha)*pow(expTerm, -11.09) + neg(alpha)*pow(expTerm, -9.));
-        tmp<volScalarField> fe2 = 1 - max(ft(magGradU), fl(magGradU));
-        tmp<volScalarField> fe = max(fe1 - 1, scalar(0))*fe2;
-
-        // Original formulation from Shur et al. paper (2008)
-        return max
-        (
-            fdTilda*(1 + fe)*lRAS + (1 - fdTilda)*lLES,
-            dimensionedScalar(dimLength, SMALL)
-        );
-    }
-
     // Simplified formulation from Gritskevich et al. paper (2011) where fe = 0
+    // return max
+    // (
+    //     fdTilda*lRAS + (1 - fdTilda)*lLES,
+    //     dimensionedScalar("SMALL", dimLength, SMALL)
+    // );
+
+    // Original formulation from Shur et al. paper (2008)
     return max
     (
-        fdTilda*lRAS + (1 - fdTilda)*lLES,
-        dimensionedScalar(dimLength, SMALL)
+        fdTilda*(1 + fe)*lRAS + (1 - fdTilda)*lLES,
+        dimensionedScalar("SMALL", dimLength, SMALL)
     );
 }
 
@@ -193,16 +222,6 @@ kOmegaSSTIDDES<BasicTurbulenceModel>::kOmegaSSTIDDES
             1.87
         )
     ),
-    fe_
-    (
-        Switch::getOrAddToDict
-        (
-            "fe",
-            this->coeffDict_,
-            true
-        )
-    ),
-
     IDDESDelta_(setDelta())
 {
     if (type == typeName)
@@ -231,17 +250,6 @@ bool kOmegaSSTIDDES<BasicTurbulenceModel>::read()
 }
 
 
-template<class BasicTurbulenceModel>
-tmp<volScalarField> kOmegaSSTIDDES<BasicTurbulenceModel>::fd() const
-{
-    const volScalarField alpha(this->alpha());
-    const volScalarField expTerm(exp(sqr(alpha)));
-
-    tmp<volScalarField> fB = min(2*pow(expTerm, -9.0), scalar(1));
-    return max(1 - fdt(mag(fvc::grad(this->U_))), fB);
-}
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace LESModels

src/TurbulenceModels/turbulenceModels/DES/kOmegaSSTIDDES/kOmegaSSTIDDES.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015 OpenFOAM Foundation
-    Copyright (C) 2019-2022 OpenCFD Ltd.
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,17 +31,15 @@ Group
     grpDESTurbulence
 
 Description
-    k-omega-SST IDDES turbulence model for
-    incompressible and compressible flows.
+    k-omega-SST IDDES turbulence model for incompressible and compressible
+    flows
 
     Reference:
     \verbatim
-        Gritskevich, M. S., Garbaruk, A. V.,
-        Schütze, J., & Menter, F. R. (2012).
-        Development of DDES and IDDES formulations for
-        the k-ω shear stress transport model.
-        Flow, turbulence and combustion, 88(3), 431-449.
-        DOI:10.1007/s10494-011-9378-4
+        Gritskevich, M.S., Garbaruk, A.V., Schuetze, J., Menter, F.R. (2011)
+        Development of DDES and IDDES Formulations for the k-omega
+        Shear Stress Transport Model, Flow, Turbulence and Combustion,
+        pp. 1-19
     \endverbatim
 
 SourceFiles
@@ -49,8 +47,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_kOmegaSSTIDDES_H
-#define Foam_kOmegaSSTIDDES_H
+#ifndef kOmegaSSTIDDES_H
+#define kOmegaSSTIDDES_H
 
 #include "kOmegaSSTDES.H"
 
@@ -79,6 +77,12 @@ class kOmegaSSTIDDES
         tmp<volScalarField> ft(const volScalarField& magGradU) const;
         tmp<volScalarField> fl(const volScalarField& magGradU) const;
 
+        tmp<volScalarField> rd
+        (
+            const volScalarField& nur,
+            const volScalarField& magGradU
+        ) const;
+
         //- Delay function
         tmp<volScalarField> fdt(const volScalarField& magGradU) const;
 
@@ -91,7 +95,7 @@ class kOmegaSSTIDDES
 
 protected:
 
-    // Protected Data
+    // Protected data
 
         // Model coefficients
 
@@ -99,16 +103,14 @@ protected:
             dimensionedScalar Cdt2_;
             dimensionedScalar Cl_;
             dimensionedScalar Ct_;
-            Switch fe_;
 
+        // Fields
 
-        //- IDDES delta
-        const IDDESDelta& IDDESDelta_;
-
+            const IDDESDelta& IDDESDelta_;
 
     // Protected Member Functions
 
-        //- Return the length scale
+        //- Length scale
         virtual tmp<volScalarField> dTilda
         (
             const volScalarField& magGradU,
@@ -151,9 +153,6 @@ public:
 
         //- Re-read model coefficients if they have changed
         virtual bool read();
-
-        //- Return the shielding function
-        virtual tmp<volScalarField> fd() const;
 };
 
 

src/TurbulenceModels/turbulenceModels/LES/LESdeltas/DeltaOmegaTildeDelta/DeltaOmegaTildeDelta.C

@@ -1,196 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2022 Upstream CFD GmbH
-    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/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "DeltaOmegaTildeDelta.H"
-#include "fvcCurl.H"
-#include "maxDeltaxyz.H"
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace LESModels
-{
-    defineTypeNameAndDebug(DeltaOmegaTildeDelta, 0);
-    addToRunTimeSelectionTable(LESdelta, DeltaOmegaTildeDelta, dictionary);
-}
-}
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-void Foam::LESModels::DeltaOmegaTildeDelta::calcDelta()
-{
-    const fvMesh& mesh = turbulenceModel_.mesh();
-
-    const volVectorField& U0 = turbulenceModel_.U();
-    tmp<volVectorField> tvorticity = fvc::curl(U0);
-    const volVectorField& vorticity = tvorticity.cref();
-    const volVectorField nvecvort
-    (
-        vorticity
-      / (
-            max
-            (
-                mag(vorticity),
-                dimensionedScalar(dimless/dimTime, SMALL)
-            )
-        )
-    );
-    tvorticity.clear();
-
-    const cellList& cells = mesh.cells();
-    const vectorField& cellCentres = mesh.cellCentres();
-    const vectorField& faceCentres = mesh.faceCentres();
-
-    forAll(cells, celli)
-    {
-        const labelList& cFaces = cells[celli];
-        const point& cc = cellCentres[celli];
-        const vector& nv = nvecvort[celli];
-
-        scalar deltaMaxTmp = 0;
-
-        for (const label facei : cFaces)
-        {
-            const point& fc = faceCentres[facei];
-            const scalar tmp = 2.0*mag(nv ^ (fc - cc));
-
-            if (tmp > deltaMaxTmp)
-            {
-                deltaMaxTmp = tmp;
-            }
-        }
-
-        delta_[celli] = deltaCoeff_*deltaMaxTmp;
-    }
-
-    const label nD = mesh.nGeometricD();
-
-    if (nD == 2)
-    {
-        WarningInFunction
-            << "Case is 2D, LES is not strictly applicable" << nl
-            << endl;
-    }
-    else if (nD != 3)
-    {
-        FatalErrorInFunction
-            << "Case must be either 2D or 3D" << exit(FatalError);
-    }
-
-    // Handle coupled boundaries
-    delta_.correctBoundaryConditions();
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::LESModels::DeltaOmegaTildeDelta::DeltaOmegaTildeDelta
-(
-    const word& name,
-    const turbulenceModel& turbulence,
-    const dictionary& dict
-)
-:
-    LESdelta(name, turbulence),
-    hmaxPtr_(nullptr),
-    deltaCoeff_
-    (
-        dict.optionalSubDict(type() + "Coeffs").getOrDefault<scalar>
-        (
-            "deltaCoeff",
-            1.035
-        )
-    ),
-    requireUpdate_
-    (
-        dict.optionalSubDict(type() + "Coeffs").getOrDefault<bool>
-        (
-            "requireUpdate",
-            true
-        )
-    )
-{
-    if (dict.optionalSubDict(type() + "Coeffs").found("hmax"))
-    {
-        // User-defined hmax
-        hmaxPtr_ =
-            LESdelta::New
-            (
-                IOobject::groupName("hmax", turbulence.U().group()),
-                turbulence,
-                dict.optionalSubDict("hmaxCoeffs"),
-                "hmax"
-            );
-    }
-    else
-    {
-        Info<< "Employing " << maxDeltaxyz::typeName << " for hmax" << endl;
-
-        hmaxPtr_.reset
-        (
-            new maxDeltaxyz
-            (
-                IOobject::groupName("hmax", turbulence.U().group()),
-                turbulence,
-                dict.optionalSubDict("hmaxCoeffs")
-            )
-        );
-    }
-
-    calcDelta();
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::LESModels::DeltaOmegaTildeDelta::read(const dictionary& dict)
-{
-    const dictionary& coeffsDict = dict.optionalSubDict(type() + "Coeffs");
-
-    coeffsDict.readIfPresent<scalar>("deltaCoeff", deltaCoeff_);
-    coeffsDict.readIfPresent<bool>("requireUpdate", requireUpdate_);
-
-    calcDelta();
-}
-
-
-void Foam::LESModels::DeltaOmegaTildeDelta::correct()
-{
-    if (turbulenceModel_.mesh().changing() && requireUpdate_)
-    {
-        hmaxPtr_->correct();
-    }
-
-    calcDelta();
-}
-
-
-// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/LESdeltas/DeltaOmegaTildeDelta/DeltaOmegaTildeDelta.H

@@ -1,142 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2022 Upstream CFD GmbH
-    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::LESModels::DeltaOmegaTildeDelta
-
-Description
-    Delta formulation that accounts for the orientation of the vorticity
-    vector. In "2D-regions" (i.e. xy-plane), delta is of the order
-    max(delta_x,delta_y), so that the influence of delta_z is discarded.
-    This can help to accelerate the transition from RANS to LES in hybrid
-    RANS/LES simulations.
-
-    Reference:
-    \verbatim
-        Shur, M. L., Spalart, P. R., Strelets, M. K., & Travin, A. K. (2015).
-        An enhanced version of DES with rapid transition
-        from RANS to LES in separated flows.
-        Flow, turbulence and combustion, 95(4), 709-737.
-        DOI:10.1007/s10494-015-9618-0
-    \endverbatim
-
-SourceFiles
-    DeltaOmegaTildeDelta.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef LESModels_DeltaOmegaTildeDelta_H
-#define LESModels_DeltaOmegaTildeDelta_H
-
-#include "LESdelta.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace LESModels
-{
-
-/*---------------------------------------------------------------------------*\
-                    Class DeltaOmegaTildeDelta Declaration
-\*---------------------------------------------------------------------------*/
-
-class DeltaOmegaTildeDelta
-:
-    public LESdelta
-{
-    // Private Data
-
-        //- Run-time selectable delta for hmax
-        //  Defaults to the maxDeltaxyz model if not supplied
-        autoPtr<LESdelta> hmaxPtr_;
-
-        //- Model coefficient
-        scalar deltaCoeff_;
-
-        //- Flag to indicate whether hmax requires updating
-        bool requireUpdate_;
-
-
-    // Private Member Functions
-
-        //- Calculate the delta values
-        void calcDelta();
-
-        //- No copy construct
-        DeltaOmegaTildeDelta(const DeltaOmegaTildeDelta&) = delete;
-
-        //- No copy assignment
-        void operator=(const DeltaOmegaTildeDelta&) = delete;
-
-
-public:
-
-    //- Runtime type information
-    TypeName("DeltaOmegaTilde");
-
-
-    // Constructors
-
-        //- Construct from name, turbulenceModel and dictionary
-        DeltaOmegaTildeDelta
-        (
-            const word& name,
-            const turbulenceModel& turbulence,
-            const dictionary&
-        );
-
-
-    //- Destructor
-    virtual ~DeltaOmegaTildeDelta() = default;
-
-
-    // Member Functions
-
-        //- Read the LESdelta dictionary
-        virtual void read(const dictionary&);
-
-        //- Return the hmax delta field
-        const volScalarField& hmax() const
-        {
-            return hmaxPtr_();
-        }
-
-        // Correct values
-        void correct();
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace LESModels
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/LESdeltas/SLADelta/SLADelta.C

@@ -1,412 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2022 Upstream CFD GmbH
-    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/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "SLADelta.H"
-#include "wallDist.H"
-#include "maxDeltaxyz.H"
-#include "fvcGrad.H"
-#include "fvcCurl.H"
-#include "addToRunTimeSelectionTable.H"
-
-#include "oneField.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-template<class Type, class WeightType = oneField>
-tmp<scalarField> sumNeighbours
-(
-    const GeometricField<Type, fvPatchField, volMesh>& field,
-    GeometricField<Type, fvPatchField, volMesh>& result,
-    const WeightType& weights = oneField()
-)
-{
-    const fvMesh& mesh = field.mesh();
-
-    result == dimensioned<Type>(field.dimensions(), Zero);
-
-    auto tscaling = tmp<scalarField>::New(mesh.nCells(), Zero);
-    auto& scaling = tscaling.ref();
-
-    const auto& faceNeighbour = mesh.faceNeighbour();
-    const auto& faceOwner = mesh.faceOwner();
-
-    forAll(faceNeighbour, i)
-    {
-        const label own = faceOwner[i];
-        const label nbr = faceNeighbour[i];
-
-        result[own] += field[nbr];
-        result[nbr] += field[own];
-
-        scaling[own] = scaling[own] + weights[nbr];
-        scaling[nbr] = scaling[nbr] + weights[own];
-    }
-
-    const auto& pbm = mesh.boundaryMesh();
-
-    for (const polyPatch& pp : pbm)
-    {
-        const auto& pf = field.boundaryField()[pp.index()];
-        const labelUList& faceCells = pp.faceCells();
-
-        if (pf.coupled())
-        {
-            const scalarField nbrField(pf.patchNeighbourField());
-
-            forAll(faceCells, facei)
-            {
-                const label celli = faceCells[facei];
-                result[celli] += nbrField[facei];
-                scaling[celli] = scaling[celli] + weights[celli];
-            }
-        }
-    }
-
-    return tscaling;
-}
-
-} // End namespace Foam
-
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace LESModels
-{
-    defineTypeNameAndDebug(SLADelta, 0);
-    addToRunTimeSelectionTable(LESdelta, SLADelta, dictionary);
-}
-}
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-void Foam::LESModels::SLADelta::calcDelta()
-{
-    const fvMesh& mesh = turbulenceModel_.mesh();
-
-    tmp<volScalarField> tnut = turbulenceModel_.nut();
-    const volScalarField& nut = tnut.cref();
-
-    tmp<volScalarField> tnu = turbulenceModel_.nu();
-    const volScalarField& nu = tnu.cref();
-
-    // Calculate vortex tilting measure, VTM
-    const volVectorField& U0 = turbulenceModel_.U();
-
-    tmp<volVectorField> tvorticity = fvc::curl(U0);
-    const volVectorField& vorticity = tvorticity.cref();
-
-    tmp<volSymmTensorField> tS = symm(fvc::grad(U0));
-    const volSymmTensorField& S = tS.cref();
-
-    const dimensionedScalar nuMin(nu.dimensions(), SMALL);
-    const dimensionedScalar eps(dimless/pow3(dimTime), SMALL);
-
-    volScalarField vtm
-    (
-        max(scalar(1), 0.2*nu/max(nut, nuMin))
-       *sqrt(6.0)*mag((S & vorticity)^vorticity)
-       /max(magSqr(vorticity)*sqrt(3*magSqr(S) - sqr(tr(S))), eps)
-    );
-    vtm.correctBoundaryConditions();
-    tS.clear();
-
-    const dimensionedScalar vortMin(dimless/dimTime, SMALL);
-    const volVectorField nVecVort(vorticity/(max(mag(vorticity), vortMin)));
-    tvorticity.clear();
-
-
-    // Calculate averaged VTM
-    volScalarField vtmAve("vtmAve", vtm);
-    tmp<scalarField> tweights = sumNeighbours(vtm, vtmAve);
-
-    // Add cell centre values
-    vtmAve += vtm;
-
-    // Weights normalisation (add 1 for cell centres)
-    vtmAve.primitiveFieldRef() /= tweights + 1;
-
-
-    // Calculate DDES shielding function, fd
-    const volScalarField& y = wallDist::New(mesh).y();
-
-    const dimensionedScalar magGradUeps(dimless/dimTime, SMALL);
-    const dimensionedScalar nuEps(nu.dimensions(), ROOTSMALL);
-
-    tmp<volScalarField> tmagGradU = mag(fvc::grad(U0));
-
-    tmp<volScalarField> trd =
-        min
-        (
-            (nut + nu)/(max(tmagGradU, magGradUeps)*sqr(kappa_*y) + nuEps),
-            scalar(10)
-        );
-    tnut.clear();
-    tnu.clear();
-
-    const volScalarField fd(1.0 - tanh(pow(Cd1_*trd, Cd2_)));
-
-
-    // Assemble delta
-    const cellList& cells = mesh.cells();
-    const vectorField& cellCentres = mesh.cellCentres();
-    const vectorField& faceCentres = mesh.faceCentres();
-
-    forAll(cells, celli)
-    {
-        const labelList& cFaces = cells[celli];
-        const point& cc = cellCentres[celli];
-        const vector& nv = nVecVort[celli];
-
-        scalar deltaMaxTmp = 0;
-
-        for (const label facei : cFaces)
-        {
-            const point& fc = faceCentres[facei];
-            const scalar tmp = 2.0*mag(nv ^ (fc - cc));
-
-            if (tmp > deltaMaxTmp)
-            {
-                deltaMaxTmp = tmp;
-            }
-        }
-
-        scalar FKH =
-            max
-            (
-                FKHmin_,
-                min
-                (
-                    FKHmax_,
-                    FKHmin_
-                  + (FKHmax_ - FKHmin_)*(vtmAve[celli] - a1_)/(a2_ - a1_)
-                )
-            );
-
-        if ((shielding_ == true) && (fd[celli] < (1.0 - epsilon_)))
-        {
-            FKH = 1;
-        }
-
-        delta_[celli] = deltaCoeff_*deltaMaxTmp*FKH;
-    }
-
-    const label nD = mesh.nGeometricD();
-
-    if (nD == 2)
-    {
-        WarningInFunction
-            << "Case is 2D, LES is not strictly applicable" << nl
-            << endl;
-    }
-    else if (nD != 3)
-    {
-        FatalErrorInFunction
-            << "Case must be either 2D or 3D" << exit(FatalError);
-    }
-
-    // Handle coupled boundaries
-    delta_.correctBoundaryConditions();
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::LESModels::SLADelta::SLADelta
-(
-    const word& name,
-    const turbulenceModel& turbulence,
-    const dictionary& dict
-)
-:
-    LESdelta(name, turbulence),
-    hmaxPtr_(nullptr),
-    deltaCoeff_
-    (
-        dict.optionalSubDict(type() + "Coeffs").getOrDefault<scalar>
-        (
-            "deltaCoeff",
-            1.035
-        )
-    ),
-    requireUpdate_
-    (
-        dict.optionalSubDict(type() + "Coeffs").getOrDefault<bool>
-        (
-            "requireUpdate",
-            true
-        )
-    ),
-    shielding_
-    (
-        dict.optionalSubDict(type() + "Coeffs").getOrDefault<bool>
-        (
-            "shielding",
-            true
-        )
-    ),
-    FKHmin_
-    (
-        dict.optionalSubDict(type() + "Coeffs").getOrDefault<scalar>
-        (
-            "FKHmin",
-            0.1
-        )
-    ),
-    FKHmax_
-    (
-        dict.optionalSubDict(type() + "Coeffs").getOrDefault<scalar>
-        (
-            "FKHmax",
-            1.0
-        )
-    ),
-    a1_
-    (
-        dict.optionalSubDict(type() + "Coeffs").getOrDefault<scalar>
-        (
-            "a1",
-            0.15
-        )
-    ),
-    a2_
-    (
-        dict.optionalSubDict(type() + "Coeffs").getOrDefault<scalar>
-        (
-            "a2",
-            0.30
-        )
-    ),
-    epsilon_
-    (
-        dict.optionalSubDict(type() + "Coeffs").getOrDefault<scalar>
-        (
-            "epsilon",
-            0.01
-        )
-    ),
-    kappa_
-    (
-        dict.optionalSubDict(type() + "Coeffs").getOrDefault<scalar>
-        (
-            "kappa",
-            0.41
-        )
-    ),
-    Cd1_
-    (
-        dict.optionalSubDict(type() + "Coeffs").getOrDefault<scalar>
-        (
-            "Cd1",
-            20
-        )
-    ),
-    Cd2_
-    (
-        dict.optionalSubDict(type() + "Coeffs").getOrDefault<scalar>
-        (
-            "Cd2",
-            3
-        )
-    )
-{
-    if (dict.optionalSubDict(type() + "Coeffs").found("hmax"))
-    {
-        // User-defined hmax
-        hmaxPtr_ =
-            LESdelta::New
-            (
-                IOobject::groupName("hmax", turbulence.U().group()),
-                turbulence,
-                dict.optionalSubDict("hmaxCoeffs"),
-                "hmax"
-            );
-    }
-    else
-    {
-        Info<< "Employing " << maxDeltaxyz::typeName << " for hmax" << endl;
-
-        hmaxPtr_.reset
-        (
-            new maxDeltaxyz
-            (
-                IOobject::groupName("hmax", turbulence.U().group()),
-                turbulence,
-                dict.optionalSubDict("hmaxCoeffs")
-            )
-        );
-    }
-
-    if (mag(a2_ - a1_) < SMALL)
-    {
-        FatalIOErrorInFunction(dict)
-            << "Model coefficients a1 = " << a1_
-            << ", and a2 = " << a2_ << " cannot be equal."
-            << abort(FatalIOError);
-    }
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::LESModels::SLADelta::read(const dictionary& dict)
-{
-    const dictionary& coeffsDict(dict.optionalSubDict(type() + "Coeffs"));
-
-    coeffsDict.readIfPresent<scalar>("deltaCoeff", deltaCoeff_);
-    coeffsDict.readIfPresent<bool>("requireUpdate", requireUpdate_);
-    coeffsDict.readIfPresent<scalar>("FKHmin", FKHmin_);
-    coeffsDict.readIfPresent<scalar>("FKHmax", FKHmax_);
-    coeffsDict.readIfPresent<scalar>("a1", a1_);
-    coeffsDict.readIfPresent<scalar>("a2", a2_);
-    coeffsDict.readIfPresent<scalar>("epsilon", epsilon_);
-    coeffsDict.readIfPresent<scalar>("kappa", kappa_);
-    coeffsDict.readIfPresent<scalar>("Cd1", Cd1_);
-    coeffsDict.readIfPresent<scalar>("Cd2", Cd2_);
-
-    calcDelta();
-}
-
-
-void Foam::LESModels::SLADelta::correct()
-{
-    if (turbulenceModel_.mesh().changing() && requireUpdate_)
-    {
-        hmaxPtr_->correct();
-    }
-
-    calcDelta();
-}
-
-
-// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/LESdeltas/SLADelta/SLADelta.H

@@ -1,149 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2022 Upstream CFD GmbH
-    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::LESModels::SLADelta
-
-Description
-    Delta formulation that accounts for the orientation of the vorticity vector
-    and a flow-sensitised function to detect 2D flow regions. Delta value is
-    strongly reduced in these regions to accelerate the transition from RANS
-    to LES in hybrid RANS/LES simulations.
-
-    Reference:
-    \verbatim
-        Shur, M. L., Spalart, P. R., Strelets, M. K., & Travin, A. K. (2015).
-        An enhanced version of DES with rapid transition
-        from RANS to LES in separated flows.
-        Flow, turbulence and combustion, 95(4), 709-737.
-        DOI:10.1007/s10494-015-9618-0
-    \endverbatim
-
-SourceFiles
-    SLADelta.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef LESModels_SLADelta_H
-#define LESModels_SLADelta_H
-
-#include "LESdelta.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace LESModels
-{
-
-/*---------------------------------------------------------------------------*\
-                          Class SLADelta Declaration
-\*---------------------------------------------------------------------------*/
-
-class SLADelta
-:
-    public LESdelta
-{
-    // Private Data
-
-        //- Run-time selectable delta for hmax
-        //  Defaults to the maxDeltaxyz model if not supplied
-        autoPtr<LESdelta> hmaxPtr_;
-
-        // Model coefficients
-
-            scalar deltaCoeff_;
-            bool requireUpdate_;
-            bool shielding_;
-            scalar FKHmin_;
-            scalar FKHmax_;
-            scalar a1_;
-            scalar a2_;
-            scalar epsilon_;
-            scalar kappa_;
-            scalar Cd1_;
-            scalar Cd2_;
-
-
-    // Private Member Functions
-
-        // Calculate the delta values
-        void calcDelta();
-
-        //- No copy construct
-        SLADelta(const SLADelta&) = delete;
-
-        //- No copy assignment
-        void operator=(const SLADelta&) = delete;
-
-
-public:
-
-    //- Runtime type information
-    TypeName("SLADelta");
-
-
-    // Constructors
-
-        //- Construct from name, turbulenceModel and IOdictionary
-        SLADelta
-        (
-            const word& name,
-            const turbulenceModel& turbulence,
-            const dictionary&
-        );
-
-
-    //- Destructor
-    virtual ~SLADelta() = default;
-
-
-    // Member Functions
-
-        //- Read the LESdelta dictionary
-        void read(const dictionary&);
-
-        //- Return the hmax delta field
-        const volScalarField& hmax() const
-        {
-            return hmaxPtr_();
-        }
-
-        // Correct values
-        void correct();
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace LESModels
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/sigma/sigma.C

@@ -1,173 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2022 Upstream CFD GmbH
-    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/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "sigma.H"
-#include "DESModel.H"
-#include "fvOptions.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace LESModels
-{
-
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-template<class BasicTurbulenceModel>
-void sigma<BasicTurbulenceModel>::correctNut()
-{
-    this->nut_ =
-        sqr(this->delta())
-       *DESModel<BasicTurbulenceModel>::Ssigma(fvc::grad(this->U_), Csigma_);
-    this->nut_.correctBoundaryConditions();
-    fv::options::New(this->mesh_).correct(this->nut_);
-
-    BasicTurbulenceModel::correctNut();
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-template<class BasicTurbulenceModel> sigma<BasicTurbulenceModel>::sigma
-(
-    const alphaField& alpha,
-    const rhoField& rho,
-    const volVectorField& U,
-    const surfaceScalarField& alphaRhoPhi,
-    const surfaceScalarField& phi,
-    const transportModel& transport,
-    const word& propertiesName,
-    const word& type
-)
-:
-    LESeddyViscosity<BasicTurbulenceModel>
-    (
-        type,
-        alpha,
-        rho,
-        U,
-        alphaRhoPhi,
-        phi,
-        transport,
-        propertiesName
-    ),
-
-    Ck_
-    (
-        dimensioned<scalar>::lookupOrAddToDict
-        (
-            "Ck",
-            this->coeffDict_,
-            0.094
-        )
-    ),
-
-    Cw_
-    (
-        dimensioned<scalar>::lookupOrAddToDict
-        (
-            "Cw",
-            this->coeffDict_,
-            0.325
-        )
-    ),
-
-    Csigma_
-    (
-        dimensioned<scalar>::lookupOrAddToDict
-        (
-            "Csigma",
-            this->coeffDict_,
-            1.68
-        )
-    )
-{
-    if (type == typeName)
-    {
-        this->printCoeffs(type);
-    }
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class BasicTurbulenceModel>
-bool sigma<BasicTurbulenceModel>::read()
-{
-    if (LESeddyViscosity<BasicTurbulenceModel>::read())
-    {
-        Ck_.readIfPresent(this->coeffDict());
-        Cw_.readIfPresent(this->coeffDict());
-        Csigma_.readIfPresent(this->coeffDict());
-
-        return true;
-    }
-
-    return false;
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> sigma<BasicTurbulenceModel>::k() const
-{
-    return tmp<volScalarField>::New
-    (
-        IOobject::groupName("k", this->U_.group()),
-        (2.0*Ck_/this->Ce_)
-       *sqr(this->delta())
-       *magSqr(dev(symm(fvc::grad(this->U_))))
-    );
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> sigma<BasicTurbulenceModel>::epsilon() const
-{
-    return tmp<volScalarField>::New
-    (
-        IOobject::groupName("epsilon", this->U_.group()),
-        this->Ce_*k()*sqrt(k())/this->delta()
-    );
-}
-
-
-template<class BasicTurbulenceModel>
-void sigma<BasicTurbulenceModel>::correct()
-{
-    LESeddyViscosity<BasicTurbulenceModel>::correct();
-    correctNut();
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace LESModels
-} // End namespace Foam
-
-// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/sigma/sigma.H

@@ -1,173 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2022 Upstream CFD GmbH
-    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::LESModels::sigma
-
-Group
-    grpLESTurbulence
-
-Description
-    The sigma SGS model.
-
-    Reference:
-    \verbatim
-        Nicoud, F., Toda, H. B., Cabrit, O., Bose, S., & Lee, J. (2011).
-        Using singular values to build a subgrid-scale
-        model for large eddy simulations.
-        Physics of fluids, 23(8), 085106.
-        DOI:10.1063/1.3623274
-    \endverbatim
-
-    The default model coefficients correspond to the following:
-    \verbatim
-        sigmaCoeffs
-        {
-            ck          0.094;
-            Csigma      1.68;
-        }
-    \endverbatim
-
-Note
-    The default value of the Csigma constant implemented was calibrated for
-    OpenFOAM using decaying isotropic turbulence and is slightly higher than
-    the value suggested in the reference publication.
-
-SourceFiles
-    sigma.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef Foam_LESModels_sigma_H
-#define Foam_LESModels_sigma_H
-
-#include "LESModel.H"
-#include "LESeddyViscosity.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace LESModels
-{
-
-/*---------------------------------------------------------------------------*\
-                           Class sigma Declaration
-\*---------------------------------------------------------------------------*/
-
-template<class BasicTurbulenceModel>
-class sigma
-:
-    public LESeddyViscosity<BasicTurbulenceModel>
-{
-    // Private Member Functions
-
-        //- No copy construct
-        sigma(const sigma&) = delete;
-
-        //- No copy assignment
-        void operator=(const sigma&) = delete;
-
-
-protected:
-
-    // Protected Data
-
-        dimensionedScalar Ck_;
-        dimensionedScalar Cw_;
-        dimensionedScalar Csigma_;
-
-
-    // Protected Member Functions
-
-        //- Update the SGS eddy-viscosity
-        virtual void correctNut();
-
-
-public:
-
-    typedef typename BasicTurbulenceModel::alphaField alphaField;
-    typedef typename BasicTurbulenceModel::rhoField rhoField;
-    typedef typename BasicTurbulenceModel::transportModel transportModel;
-
-
-    //- Runtime type information
-    TypeName("sigma");
-
-
-    // Constructors
-
-        //- Construct from components
-        sigma
-        (
-            const alphaField& alpha,
-            const rhoField& rho,
-            const volVectorField& U,
-            const surfaceScalarField& alphaRhoPhi,
-            const surfaceScalarField& phi,
-            const transportModel& transport,
-            const word& propertiesName = turbulenceModel::propertiesName,
-            const word& type = typeName
-        );
-
-
-    //- Destructor
-    virtual ~sigma()
-    {}
-
-
-    // Member Functions
-
-        //- Read model coefficients if they have changed
-        virtual bool read();
-
-        //- Return SGS kinetic energy
-        virtual tmp<volScalarField> k() const;
-
-        //- Return SGS disipation rate
-        virtual tmp<volScalarField> epsilon() const;
-
-        //- Correct Eddy-Viscosity and related properties
-        virtual void correct();
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace LESModels
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-    #include "sigma.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/Make/files

@@ -10,8 +10,7 @@ $(LESdelta)/smoothDelta/smoothDelta.C
 $(LESdelta)/maxDeltaxyz/maxDeltaxyz.C
 $(LESdelta)/IDDESDelta/IDDESDelta.C
 $(LESdelta)/maxDeltaxyzCubeRootLESDelta/maxDeltaxyzCubeRootLESDelta.C
-$(LESdelta)/DeltaOmegaTildeDelta/DeltaOmegaTildeDelta.C
-$(LESdelta)/SLADelta/SLADelta.C
+
 
 LESfilters = LES/LESfilters
 

src/TurbulenceModels/turbulenceModels/RAS/SpalartAllmaras/SpalartAllmaras.C

@@ -27,6 +27,9 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "SpalartAllmaras.H"
+#include "fvOptions.H"
+#include "bound.H"
+#include "wallDist.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -38,14 +41,101 @@ namespace RASModels
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 template<class BasicTurbulenceModel>
-Foam::tmp<Foam::volScalarField> SpalartAllmaras<BasicTurbulenceModel>::dTilda
+tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::chi() const
+{
+    return nuTilda_/this->nu();
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::fv1
 (
-    const volScalarField& chi,
-    const volScalarField& fv1,
-    const volTensorField& gradU
+    const volScalarField& chi
 ) const
 {
-    return this->y_;
+    const volScalarField chi3(pow3(chi));
+
+    return chi3/(chi3 + pow3(Cv1_));
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField::Internal> SpalartAllmaras<BasicTurbulenceModel>::fv2
+(
+    const volScalarField::Internal& chi,
+    const volScalarField::Internal& fv1
+) const
+{
+    return scalar(1) - chi/(scalar(1) + chi*fv1);
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField::Internal> SpalartAllmaras<BasicTurbulenceModel>::Stilda()
+const
+{
+    const volScalarField chi(this->chi());
+
+    const volScalarField fv1(this->fv1(chi));
+
+    const volScalarField::Internal Omega
+    (
+        ::sqrt(scalar(2))*mag(skew(fvc::grad(this->U_)().v()))
+    );
+
+    return
+    (
+        max
+        (
+            Omega + fv2(chi(), fv1())*nuTilda_()/sqr(kappa_*y_),
+            Cs_*Omega
+        )
+    );
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField::Internal> SpalartAllmaras<BasicTurbulenceModel>::fw
+(
+    const volScalarField::Internal& Stilda
+) const
+{
+    const volScalarField::Internal r
+    (
+        min
+        (
+            nuTilda_
+           /(
+               max
+               (
+                   Stilda,
+                   dimensionedScalar(Stilda.dimensions(), SMALL)
+               )
+              *sqr(kappa_*y_)
+            ),
+            scalar(10)
+        )
+    );
+
+    const volScalarField::Internal g(r + Cw2_*(pow6(r) - r));
+
+    return
+        g*pow
+        (
+            (scalar(1) + pow6(Cw3_))/(pow6(g) + pow6(Cw3_)),
+            scalar(1)/scalar(6)
+        );
+}
+
+
+template<class BasicTurbulenceModel>
+void SpalartAllmaras<BasicTurbulenceModel>::correctNut()
+{
+    this->nut_ = nuTilda_*this->fv1(this->chi());
+    this->nut_.correctBoundaryConditions();
+    fv::options::New(this->mesh_).correct(this->nut_);
+
+    BasicTurbulenceModel::correctNut();
 }
 
 
@@ -64,7 +154,7 @@ SpalartAllmaras<BasicTurbulenceModel>::SpalartAllmaras
     const word& type
 )
 :
-    SpalartAllmarasBase<eddyViscosity<RASModel<BasicTurbulenceModel>>>
+    eddyViscosity<RASModel<BasicTurbulenceModel>>
     (
         type,
         alpha,
@@ -74,7 +164,97 @@ SpalartAllmaras<BasicTurbulenceModel>::SpalartAllmaras
         phi,
         transport,
         propertiesName
-    )
+    ),
+
+    sigmaNut_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "sigmaNut",
+            this->coeffDict_,
+            scalar(2)/scalar(3)
+        )
+    ),
+    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_) + (scalar(1) + 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
+        )
+    ),
+
+    nuTilda_
+    (
+        IOobject
+        (
+            "nuTilda",
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh_
+    ),
+
+    y_(wallDist::New(this->mesh_).y())
 {
     if (type == typeName)
     {
@@ -83,6 +263,153 @@ SpalartAllmaras<BasicTurbulenceModel>::SpalartAllmaras
 }
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+bool SpalartAllmaras<BasicTurbulenceModel>::read()
+{
+    if (eddyViscosity<RASModel<BasicTurbulenceModel>>::read())
+    {
+        sigmaNut_.readIfPresent(this->coeffDict());
+        kappa_.readIfPresent(this->coeffDict());
+
+        Cb1_.readIfPresent(this->coeffDict());
+        Cb2_.readIfPresent(this->coeffDict());
+        Cw1_ = Cb1_/sqr(kappa_) + (scalar(1) + Cb2_)/sigmaNut_;
+        Cw2_.readIfPresent(this->coeffDict());
+        Cw3_.readIfPresent(this->coeffDict());
+        Cv1_.readIfPresent(this->coeffDict());
+        Cs_.readIfPresent(this->coeffDict());
+
+        return true;
+    }
+
+    return false;
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::DnuTildaEff() const
+{
+    return tmp<volScalarField>::New
+    (
+        "DnuTildaEff",
+        (nuTilda_ + this->nu())/sigmaNut_
+    );
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::k() const
+{
+    // (B:Eq. 4.50)
+    const scalar Cmu = 0.09;
+
+    return tmp<volScalarField>::New
+    (
+        IOobject
+        (
+            IOobject::groupName("k", this->alphaRhoPhi_.group()),
+            this->runTime_.timeName(),
+            this->mesh_
+        ),
+        cbrt(this->fv1(this->chi()))
+        *nuTilda_
+        *::sqrt(scalar(2)/Cmu)
+        *mag(symm(fvc::grad(this->U_))),
+        this->nut_.boundaryField().types()
+    );
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::epsilon() const
+{
+    // (B:Eq. 4.50)
+    const scalar Cmu = 0.09;
+    const dimensionedScalar nutSMALL(sqr(dimLength)/dimTime, SMALL);
+
+    return tmp<volScalarField>::New
+    (
+        IOobject
+        (
+            IOobject::groupName("epsilon", this->alphaRhoPhi_.group()),
+            this->runTime_.timeName(),
+            this->mesh_
+        ),
+        pow(this->fv1(this->chi()), 0.5)
+        *pow(::sqrt(Cmu)*this->k(), 2)
+        /(nuTilda_ + this->nut_ + nutSMALL),
+        this->nut_.boundaryField().types()
+    );
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::omega() const
+{
+    // (P:p. 384)
+    const scalar betaStar = 0.09;
+    const dimensionedScalar k0(sqr(dimLength/dimTime), SMALL);
+
+    return tmp<volScalarField>::New
+    (
+        IOobject
+        (
+            IOobject::groupName("omega", this->alphaRhoPhi_.group()),
+            this->runTime_.timeName(),
+            this->mesh_
+        ),
+        this->epsilon()/(betaStar*(this->k() + k0)),
+        this->nut_.boundaryField().types()
+    );
+}
+
+
+template<class BasicTurbulenceModel>
+void SpalartAllmaras<BasicTurbulenceModel>::correct()
+{
+    if (!this->turbulence_)
+    {
+        return;
+    }
+
+    {
+        // Construct local convenience references
+        const alphaField& alpha = this->alpha_;
+        const rhoField& rho = this->rho_;
+        const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
+        fv::options& fvOptions(fv::options::New(this->mesh_));
+
+        eddyViscosity<RASModel<BasicTurbulenceModel>>::correct();
+
+        const volScalarField::Internal Stilda(this->Stilda());
+
+        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_()
+          - fvm::Sp(Cw1_*alpha()*rho()*fw(Stilda)*nuTilda_()/sqr(y_), 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();
+    }
+
+    // Update nut with latest available nuTilda
+    correctNut();
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace RASModels

src/TurbulenceModels/turbulenceModels/RAS/SpalartAllmaras/SpalartAllmaras.H

@@ -52,6 +52,18 @@ Description
             Spalart-Allmaras One-Equation Model without ft2 Term (SA-noft2).
             https://turbmodels.larc.nasa.gov/spalart.html#sanoft2
             (Retrieved:12-01-2021).
+
+        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
 
 Usage
@@ -84,6 +96,8 @@ Note
   - The model is implemented without the trip-term since the model has almost
     always been used in fully turbulent applications rather than those where
     laminar-turbulent transition occurs.
+  - It has been argued that the \c ft2 term is not needed in the absence of the
+    trip-term, hence \c ft2 term is also not implementated.
   - The \c Stilda generation term should never be allowed to be zero or negative
     to avoid potential numerical issues and unphysical results for complex
     flows. To this end, a limiter proposed by Spalart (R:Note-1(b)) is applied
@@ -98,12 +112,11 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_SpalartAllmaras_H
-#define Foam_SpalartAllmaras_H
+#ifndef SpalartAllmaras_H
+#define SpalartAllmaras_H
 
 #include "RASModel.H"
 #include "eddyViscosity.H"
-#include "SpalartAllmarasBase.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -119,7 +132,7 @@ namespace RASModels
 template<class BasicTurbulenceModel>
 class SpalartAllmaras
 :
-    public SpalartAllmarasBase<eddyViscosity<RASModel<BasicTurbulenceModel>>>
+    public eddyViscosity<RASModel<BasicTurbulenceModel>>
 {
     // Private Member Functions
 
@@ -132,16 +145,55 @@ class SpalartAllmaras
 
 protected:
 
+    // Protected Data
+
+        // Model coefficients
+
+            dimensionedScalar sigmaNut_;
+            dimensionedScalar kappa_;
+            dimensionedScalar Cb1_;
+            dimensionedScalar Cb2_;
+            dimensionedScalar Cw1_;
+            dimensionedScalar Cw2_;
+            dimensionedScalar Cw3_;
+            dimensionedScalar Cv1_;
+            dimensionedScalar Cs_;
+
+
+        // Fields
+
+            //- Modified kinematic viscosity [m2/s]
+            volScalarField nuTilda_;
+
+            //- Wall distance
+            //  Note: different to wall distance in parent RASModel
+            //  which is for near-wall cells only
+            const volScalarField::Internal& y_;
+
+
     // Protected Member Functions
 
-        //- Return the length scale
-        virtual tmp<volScalarField> dTilda
+        tmp<volScalarField> chi() const;
+
+        tmp<volScalarField> fv1(const volScalarField& chi) const;
+
+        tmp<volScalarField::Internal> fv2
         (
-            const volScalarField& chi,
-            const volScalarField& fv1,
-            const volTensorField& gradU
+            const volScalarField::Internal& chi,
+            const volScalarField::Internal& fv1
         ) const;
 
+        tmp<volScalarField::Internal> Stilda() const;
+
+        tmp<volScalarField::Internal> fw
+        (
+            const volScalarField::Internal& Stilda
+        ) const;
+
+
+        //- Update nut with the latest available nuTilda
+        virtual void correctNut();
+
 
 public:
 
@@ -172,6 +224,27 @@ public:
 
     //- Destructor
     virtual ~SpalartAllmaras() = 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;
+
+        //- Solve the turbulence equations and correct the turbulent viscosity
+        virtual void correct();
 };
 
 

src/dynamicMesh/motionSolvers/displacement/points0/points0MotionSolver.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2017 OpenFOAM Foundation
-    Copyright (C) 2015-2022 OpenCFD Ltd.
+    Copyright (C) 2015-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -84,7 +84,7 @@ Foam::points0MotionSolver::points0MotionSolver
 )
 :
     motionSolver(mesh, dict, type),
-    zoneMotion(coeffDict(), mesh),
+    zoneMotion(dict, mesh),
     points0_(points0IO(mesh))
 {
     if
@@ -128,7 +128,7 @@ Foam::points0MotionSolver::points0MotionSolver
 )
 :
     motionSolver(mesh, dict, type),
-    zoneMotion(coeffDict(), mesh),
+    zoneMotion(dict, mesh),
     points0_(points0)
 {
     if (points0_.size() != mesh.nPoints())

src/engine/include/StCorr.H

@@ -80,7 +80,7 @@
             );
             dimensionedScalar AkEst = gSum(mgb*mesh.V().field());
 
-            StCorr.value() = max(min((Ak/AkEst).value(), 10.0), 1.0);
+            StCorr.value() = Foam::max(Foam::min((Ak/AkEst).value(), 10.0), 1.0);
 
             Info<< "StCorr = " << StCorr.value() << endl;
         }