Merge commit 'OpenCFD/master' into HEAD

applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C

@@ -50,7 +50,7 @@ PDRkEpsilon::PDRkEpsilon
     const volScalarField& rho,
     const volVectorField& U,
     const surfaceScalarField& phi,
-    basicThermo& thermophysicalModel
+    const basicThermo& thermophysicalModel
 )
 :
     RASModel(typeName, rho, U, phi, thermophysicalModel),

applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.H

@@ -93,7 +93,7 @@ public:
             const volScalarField& rho,
             const volVectorField& U,
             const surfaceScalarField& phi,
-            basicThermo& thermophysicalModel
+            const basicThermo& thermophysicalModel
         );
 
 

applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.C

@@ -129,9 +129,20 @@ void Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::updateCoeffs()
     {
         Info<< "solidWallMixedTemperatureCoupledFvPatchScalarField::"
             << "updateCoeffs() :"
+            << " patch:" << patch().name()
             << " walltemperature "
-            << " min:" << gMin(*this)
+            << " min:"
-            << " max:" << gMax(*this)
+            <<  returnReduce
+                (
+                    (this->size() > 0 ? min(*this) : VGREAT),
+                    minOp<scalar>()
+                )
+            << " max:"
+            <<  returnReduce
+                (
+                    (this->size() > 0 ? max(*this) : -VGREAT),
+                    maxOp<scalar>()
+                )
             << " avg:" << gAverage(*this)
             << endl;
     }
@@ -163,7 +174,9 @@ void Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::updateCoeffs()
         label nTotSize = returnReduce(this->size(), sumOp<label>());
 
         Info<< "solidWallMixedTemperatureCoupledFvPatchScalarField::"
-            << "updateCoeffs() : Out of " << nTotSize
+            << "updateCoeffs() :"
+            << " patch:" << patch().name()
+            << " out of:" << nTotSize
             << " fixedBC:" << nFixed
             << " gradient:" << nTotSize-nFixed << endl;
     }
@@ -213,9 +226,22 @@ void Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::evaluate
 
         if (debug)
         {
-            Info<< "Setting master and slave to wall temperature "
+            Info<< "solidWallMixedTemperatureCoupledFvPatchScalarField::"
-                << " min:" << gMin(*this)
+                << "updateCoeffs() :"
-                << " max:" << gMax(*this)
+                << " patch:" << patch().name()
+                << " setting master and slave to wall temperature "
+                << " min:"
+                <<  returnReduce
+                    (
+                        (this->size() > 0 ? min(*this) : VGREAT),
+                        minOp<scalar>()
+                    )
+                << " max:"
+                <<  returnReduce
+                    (
+                        (this->size() > 0 ? max(*this) : -VGREAT),
+                        maxOp<scalar>()
+                    )
                 << " avg:" << gAverage(*this)
                 << endl;
         }

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

@@ -488,17 +488,18 @@ labelList getNonRegionCells(const labelList& cellRegion, const label regionI)
 }
 
 
-// Get per region-region interface the sizes.
+// Get per region-region interface the sizes. If sumParallel sums sizes.
-// If sumParallel does merge.
+// Returns interfaces as straight list for looping in identical order.
-EdgeMap<label> getInterfaceSizes
+void getInterfaceSizes
 (
     const polyMesh& mesh,
     const labelList& cellRegion,
-    const bool sumParallel
+    const bool sumParallel,
+
+    edgeList& interfaces,
+    EdgeMap<label>& interfaceSizes
 )
 {
-    EdgeMap<label> interfaceSizes;
-
     forAll(mesh.faceNeighbour(), faceI)
     {
         label ownRegion = cellRegion[mesh.faceOwner()[faceI]];
@@ -585,7 +586,12 @@ EdgeMap<label> getInterfaceSizes
         }
     }
 
-    return interfaceSizes;
+    // Make sure all processors have interfaces in same order
+    interfaces = interfaceSizes.toc();
+    if (sumParallel)
+    {
+        Pstream::scatter(interfaces);
+    }
 }
 
 
@@ -705,11 +711,7 @@ autoPtr<mapPolyMesh> createRegionMesh
 
         if (otherRegion != -1)
         {
-            edge interface
+            edge interface(regionI, otherRegion);
-            (
-                min(regionI, otherRegion),
-                max(regionI, otherRegion)
-            );
 
             // Find the patch.
             if (regionI < otherRegion)
@@ -848,6 +850,7 @@ void createAndWriteRegion
 
 
     const polyBoundaryMesh& newPatches = newMesh().boundaryMesh();
+    newPatches.checkParallelSync(true);
 
     // Delete empty patches
     // ~~~~~~~~~~~~~~~~~~~~
@@ -863,15 +866,14 @@ void createAndWriteRegion
     {
         const polyPatch& pp = newPatches[patchI];
 
-        if
+        if (!isA<processorPolyPatch>(pp))
-        (
+        {
-            !isA<processorPolyPatch>(pp)
+            if (returnReduce(pp.size(), sumOp<label>()) > 0)
-         && returnReduce(pp.size(), sumOp<label>()) > 0
-        )
             {
                 oldToNew[patchI] = newI++;
             }
         }
+    }
 
     // Same for processor patches (but need no reduction)
     forAll(newPatches, patchI)
@@ -983,10 +985,15 @@ void createAndWriteRegion
 }
 
 
+// Create for every region-region interface with non-zero size two patches.
+// First one is for minimumregion to maximumregion.
+// Note that patches get created in same order on all processors (if parallel)
+// since looping over synchronised 'interfaces'.
 EdgeMap<label> addRegionPatches
 (
     fvMesh& mesh,
     const regionSplit& cellRegion,
+    const edgeList& interfaces,
     const EdgeMap<label>& interfaceSizes,
     const wordList& regionNames
 )
@@ -998,15 +1005,12 @@ EdgeMap<label> addRegionPatches
 
     EdgeMap<label> interfaceToPatch(cellRegion.nRegions());
 
-    // Keep start of added patches for later.
+    forAll(interfaces, interI)
-    label minAddedPatchI = labelMax;
-
-    forAllConstIter(EdgeMap<label>, interfaceSizes, iter)
     {
-        if (iter() > 0)
+        const edge& e = interfaces[interI];
-        {
-            const edge& e = iter.key();
 
+        if (interfaceSizes[e] > 0)
+        {
             label patchI = addPatch
             (
                 mesh,
@@ -1025,12 +1029,9 @@ EdgeMap<label> addRegionPatches
                 << " " << mesh.boundaryMesh()[patchI].name()
                 << endl;
 
-            interfaceToPatch.insert(iter.key(), patchI);
+            interfaceToPatch.insert(e, patchI);
-
-            minAddedPatchI = min(minAddedPatchI, patchI);
         }
     }
-    //Info<< "minAddedPatchI:" << minAddedPatchI << endl;
     return interfaceToPatch;
 }
 
@@ -1348,24 +1349,26 @@ int main(int argc, char *argv[])
 
     // Sizes of interface between regions. From pair of regions to number of
     // faces.
-    EdgeMap<label> interfaceSizes
+    edgeList interfaces;
-    (
+    EdgeMap<label> interfaceSizes;
     getInterfaceSizes
     (
         mesh,
         cellRegion,
-            true       // sum in parallel?
+        true,      // sum in parallel?
-        )
+
+        interfaces,
+        interfaceSizes
     );
 
     Info<< "Region\tRegion\tFaces" << nl
         << "------\t------\t-----" << endl;
 
-    forAllConstIter(EdgeMap<label>, interfaceSizes, iter)
+    forAll(interfaces, interI)
     {
-        const edge& e = iter.key();
+        const edge& e = interfaces[interI];
 
-        Info<< e[0] << '\t' << e[1] << '\t' << iter() << nl;
+        Info<< e[0] << '\t' << e[1] << '\t' << interfaceSizes[e] << nl;
     }
     Info<< endl;
 
@@ -1511,6 +1514,7 @@ int main(int argc, char *argv[])
             (
                 mesh,
                 cellRegion,
+                interfaces,
                 interfaceSizes,
                 regionNames
             )

src/OSspecific/Unix/Unix.C

@@ -91,7 +91,9 @@ Foam::string Foam::getEnv(const word& envName)
     }
     else
     {
-        return string::null;
+        // Return null-constructed string rather than string::null
+        // to avoid cyclic dependencies in the construction of globals
+        return string();
     }
 }
 
@@ -277,7 +279,9 @@ Foam::fileName Foam::findEtcFile(const fileName& name, bool mandatory)
         ::exit(1);
     }
 
-    return fileName::null;
+    // Return null-constructed fileName rather than fileName::null
+    // to avoid cyclic dependencies in the construction of globals
+    return fileName();
 }
 
 

src/OpenFOAM/global/global.Cver

@@ -38,11 +38,6 @@ Description
 const char* const Foam::FOAMversion = "VERSION_STRING";
 const char* const Foam::FOAMbuild = "BUILD_STRING";
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-// Static initializers for string::null, word::null and fileName::null
-
-#include "stringsGlobals.C"
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 // Setup an error handler for the global new operator
 

src/OpenFOAM/include/OSspecific.H

@@ -62,7 +62,7 @@ pid_t pgid();
 bool env(const word&);
 
 //- Return environment variable of given name
-//  Return string::null if the environment is undefined
+//  Return string() if the environment is undefined
 string getEnv(const word& name);
 
 //- Set an environment variable
@@ -101,7 +101,7 @@ bool chDir(const fileName& dir);
 //  -# shipped settings:
 //    - $WM_PROJECT_DIR/etc/
 //
-//  @return the full path name or fileName::null if the name cannot be found
+//  @return the full path name or fileName() if the name cannot be found
 //  Optionally abort if the file cannot be found
 fileName findEtcFile(const fileName& name, bool mandatory=false);
 

src/OpenFOAM/primitives/strings/fileName/fileName.C

@@ -33,6 +33,7 @@ License
 
 const char* const Foam::fileName::typeName = "fileName";
 int Foam::fileName::debug(debug::debugSwitch(fileName::typeName, 0));
+const Foam::fileName Foam::fileName::null;
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 

src/OpenFOAM/primitives/strings/string/string.C

@@ -29,6 +29,9 @@ License
 
 /* * * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * * */
 
+const char* const Foam::string::typeName = "string";
+int Foam::string::debug(debug::debugSwitch(string::typeName, 0));
+const Foam::string Foam::string::null;
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 

src/OpenFOAM/primitives/strings/string/string.H

@@ -80,6 +80,8 @@ public:
 
     // Static data members
 
+        static const char* const typeName;
+        static int debug;
         static const string null;
 
 
@@ -116,8 +118,6 @@ public:
 
     // Member Functions
 
-    // Access
-
         //- Count and return the number of a given character in the string
         size_type count(const char) const;
 
@@ -130,9 +130,6 @@ public:
         template<class String>
         static inline bool meta(const string&, const char quote='\\');
 
-
-    // Edit
-
         //- Strip invalid characters from the given string
         template<class String>
         static inline bool stripInvalid(string&);
@@ -145,7 +142,6 @@ public:
         template<class String>
         static inline string quotemeta(const string&, const char quote='\\');
 
-
         //- Replace first occurence of sub-string oldStr with newStr
         //  starting at start
         string& replace

src/OpenFOAM/primitives/strings/stringsGlobals.C

@@ -38,8 +38,4 @@ Description
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-const Foam::string Foam::string::null;
-const Foam::word Foam::word::null;
-const Foam::fileName Foam::fileName::null;
-
 // ************************************************************************* //

src/OpenFOAM/primitives/strings/word/word.C

@@ -31,5 +31,6 @@ License
 
 const char* const Foam::word::typeName = "word";
 int Foam::word::debug(Foam::debug::debugSwitch(word::typeName, 0));
+const Foam::word Foam::word::null;
 
 // ************************************************************************* //

src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C

@@ -211,7 +211,6 @@ Foam::refinementSurfaces::refinementSurfaces
 
             minLevel_[globalRegionI] = iter();
             maxLevel_[globalRegionI] = regionMaxLevel[surfI][iter.key()];
-            perpendicularAngle_[globalRegionI] = regionAngle[surfI][iter.key()];
 
             // Check validity
             if
@@ -231,6 +230,13 @@ Foam::refinementSurfaces::refinementSurfaces
                     << exit(FatalError);
             }
         }
+        forAllConstIter(Map<scalar>, regionAngle[surfI], iter)
+        {
+            label globalRegionI = regionOffset_[surfI] + iter.key();
+
+            perpendicularAngle_[globalRegionI] = regionAngle[surfI][iter.key()];
+        }
+
 
         //// Optional patch names and patch types
         //forAllConstIter(HashTable<word>, regionPatchName[surfI], iter)
@@ -387,7 +393,6 @@ Foam::refinementSurfaces::refinementSurfaces
 
             minLevel_[globalRegionI] = iter();
             maxLevel_[globalRegionI] = regionMaxLevel[surfI][iter.key()];
-            perpendicularAngle_[globalRegionI] = regionAngle[surfI][iter.key()];
 
             // Check validity
             if
@@ -407,6 +412,12 @@ Foam::refinementSurfaces::refinementSurfaces
                     << exit(FatalError);
             }
         }
+        forAllConstIter(Map<scalar>, regionAngle[surfI], iter)
+        {
+            label globalRegionI = regionOffset_[surfI] + iter.key();
+
+            perpendicularAngle_[globalRegionI] = regionAngle[surfI][iter.key()];
+        }
     }
 }
 

src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.C

@@ -26,7 +26,6 @@ License
 
 #include "SpalartAllmaras.H"
 #include "addToRunTimeSelectionTable.H"
-#include "wallDist.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -54,7 +53,7 @@ tmp<volScalarField> SpalartAllmaras::fv1() const
 tmp<volScalarField> SpalartAllmaras::fv2() const
 {
     volScalarField chi = nuTilda_/nu();
-    return 1.0/pow3(scalar(1) + chi/Cv2_);
+    return 1/pow3(scalar(1) + chi/Cv2_);
 }
 
 
@@ -71,29 +70,30 @@ tmp<volScalarField> SpalartAllmaras::fv3() const
 }
 
 
-tmp<volScalarField> SpalartAllmaras::calcS(const volTensorField& gradU)
+tmp<volScalarField> SpalartAllmaras::S(const volTensorField& gradU) const
 {
-    return ::sqrt(2.0)*mag(skew(gradU));
+    return sqrt(2.0)*mag(skew(gradU));
 }
 
 
-tmp<volScalarField> SpalartAllmaras::calcSTilda(const volTensorField& gradU)
+tmp<volScalarField> SpalartAllmaras::STilda
+(
+    const volScalarField& S,
+    const volScalarField& dTilda
+) const
 {
-    return fv3()*calcS(gradU) + fv2()*nuTilda_/sqr(kappa_*dTilda_);
+    return fv3()*S + fv2()*nuTilda_/sqr(kappa_*dTilda);
 }
 
 
 tmp<volScalarField> SpalartAllmaras::r
 (
     const volScalarField& visc,
-    const volScalarField& S
+    const volScalarField& S,
+    const volScalarField& dTilda
 ) const
 {
-    tmp<volScalarField> tr
+    return min
-    (
-        new volScalarField
-        (
-            min
     (
         visc
        /(
@@ -102,7 +102,7 @@ tmp<volScalarField> SpalartAllmaras::r
                S,
                dimensionedScalar("SMALL", S.dimensions(), SMALL)
            )
-                   *sqr(kappa_*dTilda_)
+          *sqr(kappa_*dTilda)
          + dimensionedScalar
            (
                "ROOTVSMALL",
@@ -110,31 +110,28 @@ tmp<volScalarField> SpalartAllmaras::r
                ROOTVSMALL
            )
         ),
-                scalar(10.0)
+        scalar(10)
-            )
-        )
     );
-
-    return tr;
 }
 
 
-tmp<volScalarField> SpalartAllmaras::fw(const volScalarField& S) const
+tmp<volScalarField> SpalartAllmaras::fw
+(
+    const volScalarField& S,
+    const volScalarField& dTilda
+) const
 {
-    volScalarField r = this->r(nuTilda_, S);
+    volScalarField r = this->r(nuTilda_, S, dTilda);
 
     volScalarField g = r + Cw2_*(pow6(r) - r);
 
-    return g*pow((1.0 + pow6(Cw3_))/(pow6(g) + pow6(Cw3_)), 1.0/6.0);
+    return g*pow((1 + pow6(Cw3_))/(pow6(g) + pow6(Cw3_)), 1.0/6.0);
 }
 
 
-void SpalartAllmaras::dTildaUpdate(const volScalarField&)
+tmp<volScalarField> SpalartAllmaras::dTilda(const volScalarField&) const
 {
-    if (mesh_.changing())
+    return min(CDES_*delta(), y_);
-    {
-        dTilda_ = min(CDES_*delta(), wallDist(mesh_).y());
-    }
 }
 
 
@@ -242,6 +239,8 @@ SpalartAllmaras::SpalartAllmaras
         )
     ),
 
+    y_(mesh_),
+
     nuTilda_
     (
         IOobject
@@ -255,8 +254,6 @@ SpalartAllmaras::SpalartAllmaras
         mesh_
     ),
 
-    dTilda_(min(CDES_*delta(), wallDist(mesh_).y())),
-
     nuSgs_
     (
         IOobject
@@ -277,11 +274,15 @@ void SpalartAllmaras::correct(const tmp<volTensorField>& gradU)
 {
     LESModel::correct(gradU);
 
-    const volScalarField STilda = calcSTilda(gradU);
+    if (mesh_.changing())
+    {
+        y_.correct();
+        y_.boundaryField() = max(y_.boundaryField(), VSMALL);
+    }
 
-    const volScalarField S = calcS(gradU);
+    const volScalarField S = this->S(gradU);
-
+    const volScalarField dTilda = this->dTilda(S);
-    dTildaUpdate(S);
+    const volScalarField STilda = this->STilda(S, dTilda);
 
     fvScalarMatrix nuTildaEqn
     (
@@ -296,7 +297,7 @@ void SpalartAllmaras::correct(const tmp<volTensorField>& gradU)
       - alphaNut_*Cb2_*magSqr(fvc::grad(nuTilda_))
      ==
         Cb1_*STilda*nuTilda_
-      - fvm::Sp(Cw1_*fw(STilda)*nuTilda_/sqr(dTilda_), nuTilda_)
+      - fvm::Sp(Cw1_*fw(STilda, dTilda)*nuTilda_/sqr(dTilda), nuTilda_)
     );
 
     nuTildaEqn.relax();
@@ -310,9 +311,15 @@ void SpalartAllmaras::correct(const tmp<volTensorField>& gradU)
 }
 
 
+tmp<volScalarField> SpalartAllmaras::k() const
+{
+    return sqr(nuSgs()/ck_/dTilda(S(fvc::grad(U()))));
+}
+
+
 tmp<volScalarField> SpalartAllmaras::epsilon() const
 {
-    return 2.0*nuEff()*magSqr(symm(fvc::grad(U())));
+    return 2*nuEff()*magSqr(symm(fvc::grad(U())));
 }
 
 

src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.H

@@ -38,6 +38,7 @@ SourceFiles
 
 #include "LESModel.H"
 #include "volFields.H"
+#include "wallDist.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -86,30 +87,39 @@ protected:
 
         // Fields
 
+            wallDist y_;
             volScalarField nuTilda_;
-            volScalarField dTilda_;
             volScalarField nuSgs_;
 
 
     // Protected member functions
 
-        // Helper functions
-
         virtual tmp<volScalarField> fv1() const;
         virtual tmp<volScalarField> fv2() const;
         virtual tmp<volScalarField> fv3() const;
-            //-
+        virtual tmp<volScalarField> S(const volTensorField& gradU) const;
-            virtual tmp<volScalarField> calcS(const volTensorField& gradU);
+
-            virtual tmp<volScalarField> calcSTilda(const volTensorField& gradU);
+        virtual tmp<volScalarField> STilda
+        (
+            const volScalarField& S,
+            const volScalarField& dTilda
+        ) const;
+
         virtual tmp<volScalarField> r
         (
             const volScalarField& visc,
-                const volScalarField& S
+            const volScalarField& S,
+            const volScalarField& dTilda
         ) const;
-            virtual tmp<volScalarField> fw(const volScalarField& S) const;
 
-        //- Length scale calculation
+        virtual tmp<volScalarField> fw
-        virtual void dTildaUpdate(const volScalarField& S);
+        (
+            const volScalarField& S,
+            const volScalarField& dTilda
+        ) const;
+
+        //- Length scale
+        virtual tmp<volScalarField> dTilda(const volScalarField& S) const;
 
 
 public:
@@ -138,10 +148,7 @@ public:
     // Member Functions
 
         //- Return SGS kinetic energy
-        virtual tmp<volScalarField> k() const
+        virtual tmp<volScalarField> k() const;
-        {
-            return sqr(nuSgs()/ck_/dTilda_);
-        }
 
         //- Return sub-grid disipation rate
         virtual tmp<volScalarField> epsilon() const;

src/turbulenceModels/incompressible/LES/SpalartAllmarasDDES/SpalartAllmarasDDES.C

@@ -26,7 +26,6 @@ License
 
 #include "SpalartAllmarasDDES.H"
 #include "addToRunTimeSelectionTable.H"
-#include "wallDist.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -50,13 +49,7 @@ tmp<volScalarField> SpalartAllmarasDDES::rd
     const volScalarField& S
 ) const
 {
-    volScalarField d = wallDist(mesh_).y();
+    return min
-
-    tmp<volScalarField> trd
-    (
-        new volScalarField
-        (
-            min
     (
         visc
         /(
@@ -64,36 +57,33 @@ tmp<volScalarField> SpalartAllmarasDDES::rd
             (
                 S,
                 dimensionedScalar("SMALL", S.dimensions(), SMALL)
-                    )*sqr(kappa_*d)
+            )*sqr(kappa_*y_)
           + dimensionedScalar
             (
                 "ROOTVSMALL",
                 dimensionSet(0, 2 , -1, 0, 0),
                 ROOTVSMALL
             )
-                ), scalar(10.0)
+        ),
-            )
+        scalar(10)
-        )
     );
-
-    return trd;
 }
 
 
-tmp<volScalarField> SpalartAllmarasDDES::fd(const volScalarField& S)
+tmp<volScalarField> SpalartAllmarasDDES::fd(const volScalarField& S) const
 {
-    return 1.0 - tanh(pow3(8.0*rd(nuSgs_ + transport().nu(), S)));
+    return 1 - tanh(pow3(8*rd(nuEff(), S)));
 }
 
 
-void SpalartAllmarasDDES::dTildaUpdate(const volScalarField& S)
+tmp<volScalarField> SpalartAllmarasDDES::dTilda(const volScalarField& S) const
 {
-    dTilda_ =
+    return max
-        wallDist(mesh_).y()
-      - fd(S)*max
     (
-            dimensionedScalar("zero", dimLength, 0.0),
+        y_
-            wallDist(mesh_).y() - CDES_*delta()
+      - fd(S)
+       *max(y_ - CDES_*delta(), dimensionedScalar("zero", dimLength, 0)),
+        dimensionedScalar("small", dimLength, SMALL)
     );
 }
 

src/turbulenceModels/incompressible/LES/SpalartAllmarasDDES/SpalartAllmarasDDES.H

@@ -62,6 +62,14 @@ class SpalartAllmarasDDES
 {
     // Private member functions
 
+        tmp<volScalarField> fd(const volScalarField& S) const;
+
+        tmp<volScalarField> rd
+        (
+            const volScalarField& visc,
+            const volScalarField& S
+        ) const;
+
         // Disallow default bitwise copy construct and assignment
         SpalartAllmarasDDES(const SpalartAllmarasDDES&);
         SpalartAllmarasDDES& operator=(const SpalartAllmarasDDES&);
@@ -71,15 +79,8 @@ protected:
 
     // Protected member functions
 
-        tmp<volScalarField> fd(const volScalarField& S);
+        //- Length scale
-        tmp<volScalarField> rd
+        virtual tmp<volScalarField> dTilda(const volScalarField& S) const;
-        (
-            const volScalarField& visc,
-            const volScalarField& S
-        ) const;
-
-        //- Length scale calculation
-        virtual void dTildaUpdate(const volScalarField& S);
 
 
 public:

src/turbulenceModels/incompressible/LES/SpalartAllmarasIDDES/SpalartAllmarasIDDES.C

@@ -26,7 +26,6 @@ License
 
 #include "SpalartAllmarasIDDES.H"
 #include "addToRunTimeSelectionTable.H"
-#include "wallDist.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -46,22 +45,29 @@ addToRunTimeSelectionTable(LESModel, SpalartAllmarasIDDES, dictionary);
 
 tmp<volScalarField> SpalartAllmarasIDDES::alpha() const
 {
-    return
+    return max
-        0.25
+    (
-      - wallDist(mesh_).y()
+        0.25 - y_/dimensionedScalar("hMax", dimLength, max(cmptMax(delta()))),
-       /dimensionedScalar("hMax", dimLength, max(cmptMax(delta())));
+        scalar(-5)
+    );
 }
 
 
-tmp<volScalarField> SpalartAllmarasIDDES::ft(const volScalarField& S) const
+tmp<volScalarField> SpalartAllmarasIDDES::ft
+(
+    const volScalarField& S
+) const
 {
-    return tanh(pow3(sqr(ct_)*r(nuSgs_, S)));
+    return tanh(pow3(sqr(ct_)*rd(nuSgs_, S)));
 }
 
 
-tmp<volScalarField> SpalartAllmarasIDDES::fl(const volScalarField& S) const
+tmp<volScalarField> SpalartAllmarasIDDES::fl
+(
+    const volScalarField& S
+) const
 {
-    return tanh(pow(sqr(cl_)*r(transport().nu(), S), 10));
+    return tanh(pow(sqr(cl_)*rd(nu(), S), 10));
 }
 
 
@@ -71,13 +77,7 @@ tmp<volScalarField> SpalartAllmarasIDDES::rd
     const volScalarField& S
 ) const
 {
-    volScalarField d = wallDist(mesh_).y();
+    return min
-
-    tmp<volScalarField> trd
-    (
-        new volScalarField
-        (
-            min
     (
         visc
        /(
@@ -85,19 +85,16 @@ tmp<volScalarField> SpalartAllmarasIDDES::rd
            (
                S,
                dimensionedScalar("SMALL", S.dimensions(), SMALL)
-                   )*sqr(kappa_*d)
+           )*sqr(kappa_*y_)
          + dimensionedScalar
            (
                "ROOTVSMALL",
                dimensionSet(0, 2 , -1, 0, 0),
                ROOTVSMALL
            )
-               ), scalar(10.0)
+       ),
-            )
+       scalar(10)
-        )
     );
-
-    return trd;
 }
 
 
@@ -105,43 +102,38 @@ tmp<volScalarField> SpalartAllmarasIDDES::rd
 
 tmp<volScalarField> SpalartAllmarasIDDES::fd(const volScalarField& S) const
 {
-    return 1.0 - tanh(pow3(8.0*rd(nuSgs_+transport().nu(), S)));
+    return 1 - tanh(pow3(8*rd(nuEff(), S)));
 }
 
 
-void SpalartAllmarasIDDES::dTildaUpdate(const volScalarField& S)
+tmp<volScalarField> SpalartAllmarasIDDES::dTilda(const volScalarField& S) const
 {
     volScalarField alpha = this->alpha();
-
     volScalarField expTerm = exp(sqr(alpha));
 
     volScalarField fHill =
-        2.0*(pos(alpha)*pow(expTerm, -11.09) + neg(alpha)*pow(expTerm, -9.0));
+        2*(pos(alpha)*pow(expTerm, -11.09) + neg(alpha)*pow(expTerm, -9.0));
 
+    volScalarField fStep = min(2*pow(expTerm, -9.0), scalar(1));
+    volScalarField fHyb = max(1 - fd(S), fStep);
+    volScalarField fAmp = 1 - max(ft(S), fl(S));
+    volScalarField fRestore = max(fHill - 1, scalar(0))*fAmp;
 
-    volScalarField fStep = min(2.0*pow(expTerm, -9.0), scalar(1));
+    // IGNORING ft2 terms
-    volScalarField fHyb = max(1.0 - fd(S), fStep);
-
-    volScalarField fAmp = 1.0 - max(ft(S), fl(S));
-
-    volScalarField fRestore = max(fHill - 1.0, scalar(0))*fAmp;
-
-    // volScalarField ft2 = IGNORING ft2 terms
-
     volScalarField Psi = sqrt
     (
         min
         (
             scalar(100),
-            (1.0 - Cb1_/(Cw1_*sqr(kappa_)*fwStar_)*fv2())/max(SMALL, fv1())
+            (1 - Cb1_/(Cw1_*sqr(kappa_)*fwStar_)*fv2())/max(SMALL, fv1())
         )
     );
 
-    dTilda_ = max
+    return max
     (
-        dimensionedScalar("zero", dimLength, 0.0),
+        dimensionedScalar("SMALL", dimLength, SMALL),
-        fHyb*(1.0 + fRestore*Psi)*wallDist(mesh_).y()
+        fHyb*(1 + fRestore*Psi)*y_
-      + (1.0 - fHyb)*CDES_*Psi*delta()
+      + (1 - fHyb)*CDES_*Psi*delta()
     );
 }
 

src/turbulenceModels/incompressible/LES/SpalartAllmarasIDDES/SpalartAllmarasIDDES.H

@@ -69,12 +69,16 @@ class SpalartAllmarasIDDES
         tmp<volScalarField> alpha() const;
         tmp<volScalarField> ft(const volScalarField& S) const;
         tmp<volScalarField> fl(const volScalarField& S) const;
+
         tmp<volScalarField> rd
         (
             const volScalarField& visc,
             const volScalarField& S
         ) const;
 
+        //- Delay function
+        tmp<volScalarField> fd(const volScalarField& S) const;
+
         // Disallow default bitwise copy construct and assignment
         SpalartAllmarasIDDES(const SpalartAllmarasIDDES&);
         SpalartAllmarasIDDES& operator=(const SpalartAllmarasIDDES&);
@@ -84,11 +88,9 @@ protected:
 
     // Protected member functions
 
-        //- Delay function
+        //- Length scale
-        tmp<volScalarField> fd(const volScalarField& S) const;
+        virtual tmp<volScalarField> dTilda(const volScalarField& S) const;
 
-        //- Length scale calculation
-        virtual void dTildaUpdate(const volScalarField& S);
 
 public: