FEATURE: Pre-release bug fixes to new features. Author: Hrvoje Jasak. Merge: Dominik Christ.

applications/solvers/compressible/sonicDyMFoam/sonicDyMFoam.C

@@ -25,12 +25,15 @@ Application
     sonicDyMFoam
 
 Description
-    Transient solver for trans-sonic/supersonic for laminar or flow of a
+    Transient solver for trans-sonic/supersonic for laminar or turbulent
-    compressible gas with support for mesh motion and topological changes
+    flow of a compressible gas with support for mesh motion and
+    topological changes
 
     Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
     pseudo-transient simulations.
 
+    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
+
     Updated from sonicFoamAutoMotion by Hrvoje Jasak
 
 Author
@@ -82,9 +85,10 @@ int main(int argc, char *argv[])
         // Mesh motion update
         if (meshChanged)
         {
-            T.correctBoundaryConditions();
+            T = max(T, TMin);
-            p.correctBoundaryConditions();
+            p = max(p, pMin);
-            e.correctBoundaryConditions();
+            e = max(e, thermo.Cv()*TMin);
+
             thermo.correct();
             rho = thermo.rho();
 

src/coupledMatrix/coupledLduMatrix/coupledLduMatrix.C

@@ -195,7 +195,8 @@ void Foam::coupledLduMatrix::initMatrixInterfaces
     const lduInterfaceFieldPtrsListList& interfaces,
     const FieldField<Field, scalar>& x,
     FieldField<Field, scalar>& result,
-    const direction cmpt
+    const direction cmpt,
+    const bool switchToLhs
 ) const
 {
     const PtrList<lduMatrix>& matrices = *this;
@@ -205,7 +206,7 @@ void Foam::coupledLduMatrix::initMatrixInterfaces
     // interfaces.  The reason is that non-processor coupled
     // interfaces require a complex comms pattern involving more than
     // pairwise communications.
-    // Under normal circumstances this is achieved naturall, since
+    // Under normal circumstances this is achieved naturally, since
     // processor interfaces come last on the list and other coupled
     // interfaces execute complex comms at init() level.
     // For coupled matrices, the update loop needs to be split over
@@ -244,22 +245,66 @@ void Foam::coupledLduMatrix::initMatrixInterfaces
                             (
                                 Pstream::defaultCommsType()
                             ),
-                            false
+                            switchToLhs
                         );
                     }
                 }
             }
         }
+        else if (Pstream::defaultCommsType() == Pstream::scheduled)
+        {
+            // ERROR: Does not work with scheduled comms.
+            // To investigate.  HJ, 11/Jun/2015
+            forAll (matrices, rowI)
+            {
+                const lduSchedule& patchSchedule =
+                    matrices[rowI].patchSchedule();
+
+                // Loop over the "global" patches are on the list of
+                // interfaces but beyond the end of the schedule
+                // which only handles "normal" patches
+                for
+                (
+                    label interfaceI = patchSchedule.size()/2;
+                    interfaceI < interfaces.size();
+                    interfaceI++
+                )
+                {
+                    if (interfaces[rowI].set(interfaceI))
+                    {
+                        if
+                        (
+                            !isA<processorLduInterfaceField>
+                            (
+                                interfaces[rowI][interfaceI]
+                            )
+                        )
+                        {
+                            interfaces[rowI][interfaceI].
+                                initInterfaceMatrixUpdate
+                                (
+                                    x[rowI],
+                                    result[rowI],
+                                    matrices[rowI],
+                                    coupleCoeffs[rowI][interfaceI],
+                                    cmpt,
+                                    static_cast<const Pstream::commsTypes>
+                                    (
+                                        Pstream::defaultCommsType()
+                                    ),
+                                    switchToLhs
+                                );
+                        }
+                    }
+                }
+            }
+        }
         else
         {
-            matrices[rowI].initMatrixInterfaces
+            FatalErrorIn("void coupledLduMatrix::initMatrixInterfaces")
-            (
+                << "Unsuported communications type "
-                coupleCoeffs[rowI],
+                << Pstream::commsTypeNames[Pstream::defaultCommsType()]
-                interfaces[rowI],
+                << exit(FatalError);
-                x[rowI],
-                result[rowI],
-                cmpt
-            );
         }
     }
 
@@ -295,12 +340,67 @@ void Foam::coupledLduMatrix::initMatrixInterfaces
                             (
                                 Pstream::defaultCommsType()
                             ),
-                            false
+                            switchToLhs
                         );
                     }
                 }
             }
         }
+        else if (Pstream::defaultCommsType() == Pstream::scheduled)
+        {
+            // ERROR: Does not work with scheduled comms.
+            // To investigate.  HJ, 11/Jun/2015
+            forAll (matrices, rowI)
+            {
+                const lduSchedule& patchSchedule =
+                    matrices[rowI].patchSchedule();
+
+                // Loop over the "global" patches are on the list of
+                // interfaces but beyond the end of the schedule
+                // which only handles "normal" patches
+                for
+                (
+                    label interfaceI = patchSchedule.size()/2;
+                    interfaceI < interfaces.size();
+                    interfaceI++
+                )
+                {
+                    if (interfaces[rowI].set(interfaceI))
+                    {
+                        if
+                        (
+                            isA<processorLduInterfaceField>
+                            (
+                                interfaces[rowI][interfaceI]
+                            )
+                        )
+                        {
+                            interfaces[rowI][interfaceI].
+                                initInterfaceMatrixUpdate
+                                (
+                                    x[rowI],
+                                    result[rowI],
+                                    matrices[rowI],
+                                    coupleCoeffs[rowI][interfaceI],
+                                    cmpt,
+                                    static_cast<const Pstream::commsTypes>
+                                    (
+                                        Pstream::defaultCommsType()
+                                    ),
+                                    switchToLhs
+                                );
+                        }
+                    }
+                }
+            }
+        }
+        else
+        {
+            FatalErrorIn("void coupledLduMatrix::initMatrixInterfaces")
+                << "Unsuported communications type "
+                << Pstream::commsTypeNames[Pstream::defaultCommsType()]
+                << exit(FatalError);
+        }
     }
 }
 
@@ -311,7 +411,8 @@ void Foam::coupledLduMatrix::updateMatrixInterfaces
     const lduInterfaceFieldPtrsListList& interfaces,
     const FieldField<Field, scalar>& x,
     FieldField<Field, scalar>& result,
-    const direction cmpt
+    const direction cmpt,
+    const bool switchToLhs
 ) const
 {
     const PtrList<lduMatrix>& matrices = *this;
@@ -324,7 +425,8 @@ void Foam::coupledLduMatrix::updateMatrixInterfaces
             interfaces[rowI],
             x[rowI],
             result[rowI],
-            cmpt
+            cmpt,
+            switchToLhs
         );
     }
 }

src/coupledMatrix/coupledLduMatrix/coupledLduMatrix.H

@@ -128,7 +128,8 @@ public:
                 const lduInterfaceFieldPtrsListList& interfaces,
                 const FieldField<Field, scalar>& x,
                 FieldField<Field, scalar>& result,
-                const direction cmpt
+                const direction cmpt,
+                const bool switchToLhs = false
             ) const;
 
             //- Update coupled interfaces for matrix operations
@@ -138,7 +139,8 @@ public:
                 const lduInterfaceFieldPtrsListList& interfaces,
                 const FieldField<Field, scalar>& x,
                 FieldField<Field, scalar>& result,
-                const direction cmpt
+                const direction cmpt,
+                const bool switchToLhs = false
             ) const;
 
 };

src/coupledMatrix/coupledLduPrecon/GaussSeidelPrecon/coupledGaussSeidelPrecon.C

@@ -340,7 +340,8 @@ void Foam::coupledGaussSeidelPrecon::precondition
                 interfaces_,
                 x,
                 bPrime_,
-                cmpt
+                cmpt,
+                true         // switch to lhs
             );
 
             matrix_.updateMatrixInterfaces
@@ -349,7 +350,8 @@ void Foam::coupledGaussSeidelPrecon::precondition
                 interfaces_,
                 x,
                 bPrime_,
-                cmpt
+                cmpt,
+                true         // switch to lhs
             );
         }
 

src/decompositionMethods/decompositionMethods/patchConstrainedDecomp/patchConstrainedDecomp.C

@@ -103,6 +103,11 @@ Foam::labelList Foam::patchConstrainedDecomp::decompose
                 << abort(FatalError);
         }
 
+        Info<< "Putting patch named " << patchConstraints_[i].first()
+            << " index " << patchID
+            << " to processor " << procID
+            << endl;
+
         const labelList fc = mesh_.boundaryMesh()[patchID].faceCells();
 
         forAll (fc, fcI)

src/dynamicMesh/dynamicMesh/polyTopoChange/polyTopoChanger/polyTopoChangerChangeMesh.C

@@ -276,8 +276,8 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChanger::changeMesh
     const labelHashSet& removedFaces = ref.removedFaces();
     const labelHashSet& removedCells = ref.removedCells();
 
-    // Grab the untouched points.
+    // Grab the untouched live points
-    forAll (points, pointI)
+    for (label pointI = 0; pointI < nOldPoints; pointI++)
     {
         // Check if the point has been removed; if not add it to the list
         if (!removedPoints.found(pointI))
@@ -373,6 +373,32 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChanger::changeMesh
         debugPointCounter = nNewPoints;
     }
 
+    // Grab the untouched auxiliary points
+    for (label pointI = nOldPoints; pointI < points.size(); pointI++)
+    {
+        // Check if the point has been removed; if not add it to the list
+        if (!removedPoints.found(pointI))
+        {
+            // Grab a point
+            newPointsZeroVol[nNewPoints] = points[pointI];
+            newPointsMotion[nNewPoints] = points[pointI];
+
+            // Grab addressing
+            renumberPoints[pointI] = nNewPoints;
+            pointMap[nNewPoints] = pointI;
+
+            nNewPoints++;
+        }
+    }
+
+    if (debug)
+    {
+        Pout<< "Added retired points: untouched = "
+            << nNewPoints - debugPointCounter;
+
+        debugPointCounter = nNewPoints;
+    }
+
     // Grab auxiliary points
     forAll (mp, mpI)
     {
@@ -922,6 +948,29 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChanger::changeMesh
     }
 
     // Add freely-standing faces to the back of the list
+
+    // Freely-standing faces from the original mesh
+    // Insert untouched internal faces
+    for (label faceI = nOldFaces; faceI < faces.size(); faceI++)
+    {
+        if (!removedFaces.found(faceI))
+        {
+            newFaces[nNewFaces] = faces[faceI];
+            renumberFaces[faceI] = nNewFaces;
+            faceMap[nNewFaces]= faceI;
+            nNewFaces++;
+        }
+    }
+
+    if (debug)
+    {
+        Pout<< "Added zone-only faces: untouched = "
+            << nNewFaces - debugFaceCounter;
+
+            debugFaceCounter = nNewFaces;
+    }
+
+    // Freely-standing modified faces
     forAll (mf, mfI)
     {
         if (mf[mfI].onlyInZone())
@@ -935,12 +984,12 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChanger::changeMesh
 
     if (debug)
     {
-        Pout<< "Added zone-only faces: modified = "
+        Pout<< " modified = " << nNewFaces - debugFaceCounter;
-            << nNewFaces - debugFaceCounter;
 
             debugFaceCounter = nNewFaces;
     }
 
+    // Freely-standing added faces
     forAll (af, afI)
     {
         if (af[afI].onlyInZone())

src/finiteVolume/fields/fvPatchFields/constraint/mixingPlane/mixingPlaneFvPatchField.C

@@ -205,7 +205,10 @@ mixingPlaneFvPatchField<Type>::mixingPlaneFvPatchField
     phiName_("phi"),
     fluxMask_(),
     fluxWeights_(p.size(), 0)
-{}
+{
+    // Cannot read mixing type here: internal field may not be set
+    // HJ, 3/Jun/2015
+}
 
 
 template<class Type>
@@ -246,7 +249,7 @@ mixingPlaneFvPatchField<Type>::mixingPlaneFvPatchField
     }
 
     // Read mixing type
-    readMixingType();
+    this->readMixingType();
 }
 
 
@@ -387,6 +390,12 @@ tmp<Field<Type> > mixingPlaneFvPatchField<Type>::patchNeighbourField() const
     // Get shadow patch internalField field
     Field<Type> sField = this->shadowPatchField().patchInternalField();
 
+    // If mixing type is unknown, read it
+    if (mixing_ == mixingPlaneInterpolation::MIXING_UNKNOWN)
+    {
+        this->readMixingType();
+    }
+
     if (mixing_ == mixingPlaneInterpolation::AREA_AVERAGING)
     {
         // Area-weighted averaging
@@ -454,6 +463,12 @@ void mixingPlaneFvPatchField<Type>::initEvaluate
 
     const scalarField& w = this->patch().weights();
 
+    // If mixing type is unknown, read it
+    if (mixing_ == mixingPlaneInterpolation::MIXING_UNKNOWN)
+    {
+        this->readMixingType();
+    }
+
     if
     (
         mixing_ == mixingPlaneInterpolation::AREA_AVERAGING
@@ -505,6 +520,12 @@ tmp<Field<Type> > mixingPlaneFvPatchField<Type>::valueInternalCoeffs
     const tmp<scalarField>& w
 ) const
 {
+    // If mixing type is unknown, read it
+    if (mixing_ == mixingPlaneInterpolation::MIXING_UNKNOWN)
+    {
+        this->readMixingType();
+    }
+
     if (mixing_ == mixingPlaneInterpolation::AREA_AVERAGING)
     {
         return pTraits<Type>::one*w;
@@ -560,6 +581,12 @@ tmp<Field<Type> > mixingPlaneFvPatchField<Type>::valueBoundaryCoeffs
     const tmp<scalarField>& w
 ) const
 {
+    // If mixing type is unknown, read it
+    if (mixing_ == mixingPlaneInterpolation::MIXING_UNKNOWN)
+    {
+        this->readMixingType();
+    }
+
     if (mixing_ == mixingPlaneInterpolation::AREA_AVERAGING)
     {
         return pTraits<Type>::one*(1.0 - w);
@@ -606,6 +633,12 @@ tmp<Field<Type> > mixingPlaneFvPatchField<Type>::valueBoundaryCoeffs
 // tmp<Field<Type> >
 // mixingPlaneFvPatchField<Type>::gradientInternalCoeffs() const
 // {
+//     // If mixing type is unknown, read it
+//     if (mixing_ == mixingPlaneInterpolation::MIXING_UNKNOWN)
+//     {
+//         this->readMixingType();
+//     }
+//
 //     if
 //     (
 //         mixing_ == mixingPlaneInterpolation::AREA_AVERAGING
@@ -645,6 +678,12 @@ tmp<Field<Type> > mixingPlaneFvPatchField<Type>::valueBoundaryCoeffs
 // tmp<Field<Type> >
 // mixingPlaneFvPatchField<Type>::gradientBoundaryCoeffs() const
 // {
+//     // If mixing type is unknown, read it
+//     if (mixing_ == mixingPlaneInterpolation::MIXING_UNKNOWN)
+//     {
+//         this->readMixingType();
+//     }
+//
 //     if
 //     (
 //         mixing_ == mixingPlaneInterpolation::AREA_AVERAGING
@@ -691,8 +730,11 @@ void mixingPlaneFvPatchField<Type>::patchInterpolate
     // HJ, 13/Jun/2013
     const label patchI = this->patch().index();
 
-    // Read mixing type
+    // If mixing type is unknown, read it
-    readMixingType();
+    if (mixing_ == mixingPlaneInterpolation::MIXING_UNKNOWN)
+    {
+        this->readMixingType();
+    }
 
     // Use circumferential average of internal field when interpolating to
     // patch.  HJ and MB, 13/Jun/2013
@@ -742,8 +784,11 @@ void mixingPlaneFvPatchField<Type>::patchInterpolate
     // HJ, 13/Jun/2013
     const label patchI = this->patch().index();
 
-    // Read mixing type
+    // If mixing type is unknown, read it
-    readMixingType();
+    if (mixing_ == mixingPlaneInterpolation::MIXING_UNKNOWN)
+    {
+        this->readMixingType();
+    }
 
     // Use circumferential average of internal field when interpolating to
     // patch.  HJ and MB, 13/Jun/2013
@@ -791,8 +836,11 @@ void mixingPlaneFvPatchField<Type>::patchFlux
 {
     const label patchI = this->patch().index();
 
-    // Read mixing type
+    // If mixing type is unknown, read it
-    readMixingType();
+    if (mixing_ == mixingPlaneInterpolation::MIXING_UNKNOWN)
+    {
+        this->readMixingType();
+    }
 
     if
     (
@@ -867,6 +915,12 @@ void mixingPlaneFvPatchField<Type>::initInterfaceMatrixUpdate
     // Communication is allowed either before or after processor
     // patch comms.  HJ, 11/Jul/2011
 
+    // If mixing type is unknown, read it
+    if (mixing_ == mixingPlaneInterpolation::MIXING_UNKNOWN)
+    {
+        this->readMixingType();
+    }
+
     if
     (
         mixing_ == mixingPlaneInterpolation::AREA_AVERAGING

src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresGrad.C

@@ -94,7 +94,7 @@ leastSquaresGrad<Type>::grad
     const unallocLabelList& own = mesh.owner();
     const unallocLabelList& nei = mesh.neighbour();
 
-    forAll(own, facei)
+    forAll (own, facei)
     {
         register label ownFaceI = own[facei];
         register label neiFaceI = nei[facei];
@@ -106,7 +106,7 @@ leastSquaresGrad<Type>::grad
     }
 
     // Boundary faces
-    forAll(vsf.boundaryField(), patchi)
+    forAll (vsf.boundaryField(), patchi)
     {
         const fvsPatchVectorField& patchOwnLs = ownLs.boundaryField()[patchi];
 
@@ -118,7 +118,7 @@ leastSquaresGrad<Type>::grad
             Field<Type> neiVsf =
                 vsf.boundaryField()[patchi].patchNeighbourField();
 
-            forAll(neiVsf, patchFaceI)
+            forAll (neiVsf, patchFaceI)
             {
                 lsGrad[faceCells[patchFaceI]] +=
                     patchOwnLs[patchFaceI]
@@ -129,7 +129,7 @@ leastSquaresGrad<Type>::grad
         {
             const fvPatchField<Type>& patchVsf = vsf.boundaryField()[patchi];
 
-            forAll(patchVsf, patchFaceI)
+            forAll (patchVsf, patchFaceI)
             {
                 lsGrad[faceCells[patchFaceI]] +=
                      patchOwnLs[patchFaceI]

src/foam/containers/Lists/IndirectList/IndirectListI.H

@@ -34,7 +34,23 @@ inline Foam::IndirectList<T>::IndirectList
 :
     completeList_(const_cast<UList<T>&>(completeList)),
     addressing_(addr)
-{}
+{
+#   ifdef FULLDEBUG
+    if (completeList_.empty() && !addressing_.empty())
+    {
+        FatalErrorIn
+        (
+            "inline Foam::IndirectList<T>::IndirectList\n"
+            "(\n"
+            "    const UList<T>& completeList,\n"
+            "    const UList<label>& addr\n"
+            ")"
+        )   << "Incorrect definition of indirect list.  "
+            << "Complete list is empty and addressing is not"
+            << abort(FatalError);
+    }
+#   endif
+}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

src/foam/db/IOstreams/Pstreams/Pstream.C

@@ -269,7 +269,9 @@ const Foam::debug::optimisationSwitch
 Foam::Pstream::defaultCommsType
 (
     "commsType",
-    "blocking",
+    "nonBlocking",
+//     "scheduled",
+//     "blocking",
     "blocking, nonBlocking, scheduled"
 );
 

src/foam/global/controlSwitches/constantsSwitch.H

@@ -59,6 +59,15 @@ class constantsSwitch
 :
     public controlSwitches<scalar>
 {
+    // Private member functions
+
+        //- Disallow construct as copy
+        constantsSwitch(const constantsSwitch&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const constantsSwitch&);
+
+
 public:
 
     // Constructors

src/foam/global/controlSwitches/controlSwitches.C

@@ -39,20 +39,6 @@ namespace debug
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-template<class Type>
-Foam::debug::controlSwitches<Type>::controlSwitches()
-:
-    switchValue_(Type(0)),
-    switchDescription_("")
-{}
-
-template<class T>
-Foam::debug::controlSwitches<T>::controlSwitches(const T& switchValue)
-:
-    switchValue_(switchValue)
-{}
-
-
 template<class T>
 Foam::debug::controlSwitches<T>::controlSwitches
 (
@@ -101,19 +87,6 @@ Foam::debug::controlSwitches<T>::controlSwitches
 }
 
 
-template<class T>
-Foam::debug::controlSwitches<T>::controlSwitches
-(
-    const debug::controlSwitches<T>& csw
-)
-:
-    switchName_(csw.switchName_),
-    switchValue_(csw.switchValue_),
-    switchDescription_(csw.switchDescription_)
-{}
-
-
-
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 template<class T>
@@ -137,37 +110,7 @@ Foam::debug::controlSwitches<T>::~controlSwitches()
 }
 
 
-// * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
-
-template<class T>
-void Foam::debug::controlSwitches<T>::operator=
-(
-    const debug::controlSwitches<T>& rhs
-)
-{
-    // Check for assignment to self
-    if (this == &rhs)
-    {
-        std::cerr
-            << "Foam::debug::controlSwitches<T>::operator="
-            << "(const Foam::controlSwitches<T>&)"
-            << "--> FOAM FATAL ERROR: "
-            << "Attempted assignment to self"
-            << std::endl;
-
-        std::abort();
-
-        exit(-1);
-    }
-    else
-    {
-        switchValue_ = rhs.switchValue_;
-        switchDescription_ = rhs.switchDescription_;
-    }
-}
-
-
-// * * * * * * * * * * * * * * * *  IOStream operators * * * * * * * * * * * //
 
 template<class T>
 void printControlSwitches

src/foam/global/controlSwitches/controlSwitches.H

@@ -231,16 +231,19 @@ class controlSwitches
             switchValuesTable_;
 
 
+    // Private member functions
+
+        //- Disallow construct as copy
+        controlSwitches(const controlSwitches<T>&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const controlSwitches<T>&);
+
+
 public:
 
     // Constructors
 
-        //- Construct null
-        controlSwitches();
-
-        //- Construct from components
-        explicit controlSwitches(const T& data);
-
         controlSwitches
         (
             const std::string& switchName,
@@ -251,9 +254,6 @@ public:
             switchesValues
         );
 
-        //- Construct as copy
-        controlSwitches(const controlSwitches&);
-
 
     //- Destructor
     ~controlSwitches();
@@ -262,7 +262,7 @@ public:
     // Member Functions
 
         // Check
-        bool boolean_test() const
+        bool test() const
         {
             // Perform Boolean logic here
             return switchValue_ != 0;
@@ -274,9 +274,6 @@ public:
         //- Assignement operator
         void operator=(const T&);
 
-        //- Assignement operator
-        void operator=(const controlSwitches&);
-
         //- & operator
         const T operator&(const T&);
 

src/foam/global/controlSwitches/debugSwitch.H

@@ -58,6 +58,15 @@ class debugSwitch
 :
     public controlSwitches<int>
 {
+    // Private member functions
+
+        //- Disallow construct as copy
+        debugSwitch(const debugSwitch&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const debugSwitch&);
+
+
 public:
 
     // Constructors

src/foam/global/controlSwitches/infoSwitch.H

@@ -54,8 +54,18 @@ ListInfoControlSwitches;
 extern ListInfoControlSwitches* infoSwitchValues_;
 
 class infoSwitch
-    : public controlSwitches<int>
+:
+    public controlSwitches<int>
 {
+    // Private member functions
+
+        //- Disallow construct as copy
+        infoSwitch(const infoSwitch&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const infoSwitch&);
+
+
 public:
 
     // Constructors

src/foam/global/controlSwitches/optimisationSwitch.H

@@ -59,6 +59,15 @@ class optimisationSwitch
 :
     public controlSwitches<int>
 {
+    // Private member functions
+
+        //- Disallow construct as copy
+        optimisationSwitch(const optimisationSwitch&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const optimisationSwitch&);
+
+
 public:
 
     // Constructors

src/foam/global/controlSwitches/safeBool.H

@@ -78,7 +78,7 @@ public:
 
     operator boolType() const
     {
-        return (static_cast<const T*>(this))->boolean_test()
+        return (static_cast<const T*>(this))->test()
             ? &safeBoolBase::thisTypeDoesNotSupportComparisons : 0;
     }
 

src/foam/global/controlSwitches/tolerancesSwitch.H

@@ -58,6 +58,15 @@ class tolerancesSwitch
 :
     public controlSwitches<scalar>
 {
+    // Private member functions
+
+        //- Disallow construct as copy
+        tolerancesSwitch(const tolerancesSwitch&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const tolerancesSwitch&);
+
+
 public:
 
     // Constructors

src/foam/matrices/blockLduMatrix/BlockAmg/coarseBlockAmgLevel.C

@@ -222,7 +222,7 @@ void Foam::coarseBlockAmgLevel<Type>::solve
     }
 
     // Switch of debug in top-level direct solve
-    debug::debugSwitch oldDebug = BlockLduMatrix<Type>::debug;
+    label oldDebug = BlockLduMatrix<Type>::debug();
 
     if (BlockLduMatrix<Type>::debug >= 4)
     {

src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrix.C

@@ -351,35 +351,54 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const BlockLduMatrix<Type>& ldum)
     else
     {
         // Dummy write for consistency
-        os.writeKeyword("diagonal")  << typename BlockLduMatrix<Type>::TypeCoeffField
+        os.writeKeyword("diagonal")
-            (ldum.lduAddr().size()) << token::END_STATEMENT << nl;
+            << typename BlockLduMatrix<Type>::TypeCoeffField
+               (
+                   ldum.lduAddr().size()
+               )
+            << token::END_STATEMENT << nl;
     }
 
     if (ldum.upperPtr_)
     {
-        os.writeKeyword("upper")  << *ldum.upperPtr_ <<  token::END_STATEMENT << nl;
+        os.writeKeyword("upper")
+            << *ldum.upperPtr_
+            << token::END_STATEMENT << nl;
     }
     else
     {
         // Dummy write for consistency
-        os.writeKeyword("upper")  << typename BlockLduMatrix<Type>::TypeCoeffField
+        os.writeKeyword("upper")
-            (ldum.lduAddr().lowerAddr().size()) <<  token::END_STATEMENT << nl;
+            << typename BlockLduMatrix<Type>::TypeCoeffField
+               (
+                   ldum.lduAddr().lowerAddr().size()
+               )
+            <<  token::END_STATEMENT << nl;
     }
 
     if (ldum.lowerPtr_)
     {
-        os.writeKeyword("lower")  << *ldum.lowerPtr_ <<  token::END_STATEMENT << nl;
+        os.writeKeyword("lower")
+            << *ldum.lowerPtr_ <<  token::END_STATEMENT << nl;
     }
     else
     {
         // Dummy write for consistency
-        os.writeKeyword("lower")  << typename BlockLduMatrix<Type>::TypeCoeffField
+        os.writeKeyword("lower")
-            (ldum.lduAddr().lowerAddr().size()) <<  token::END_STATEMENT << nl;
+            << typename BlockLduMatrix<Type>::TypeCoeffField
+               (
+                   ldum.lduAddr().lowerAddr().size()
+               )
+            <<  token::END_STATEMENT << nl;
     }
 
-    os.writeKeyword("coupleUpper") << ldum.coupleUpper_ << token::END_STATEMENT << endl;
+    os.writeKeyword("coupleUpper")
+        << ldum.coupleUpper_
+        << token::END_STATEMENT << endl;
 
-    os.writeKeyword("coupleLower") << ldum.coupleLower_ << token::END_STATEMENT << endl;
+    os.writeKeyword("coupleLower")
+        << ldum.coupleLower_
+        << token::END_STATEMENT << endl;
 
     os.check("Ostream& operator<<(Ostream&, const BlockLduMatrix<Type>&");
 

src/foam/matrices/blockLduMatrix/BlockLduMatrix/blockLduMatrices.C

@@ -38,11 +38,12 @@ namespace Foam
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-defineNamedTemplateTypeNameAndDebug(blockScalarMatrix, 0);
+// Note: Matrix debug level defaults to 1.  HJ, 4/Jun/2015
-defineNamedTemplateTypeNameAndDebug(blockVectorMatrix, 0);
+defineNamedTemplateTypeNameAndDebug(blockScalarMatrix, 1);
-defineNamedTemplateTypeNameAndDebug(blockSphericalTensorMatrix, 0);
+defineNamedTemplateTypeNameAndDebug(blockVectorMatrix, 1);
-defineNamedTemplateTypeNameAndDebug(blockSymmTensorMatrix, 0);
+defineNamedTemplateTypeNameAndDebug(blockSphericalTensorMatrix, 1);
-defineNamedTemplateTypeNameAndDebug(blockTensorMatrix, 0);
+defineNamedTemplateTypeNameAndDebug(blockSymmTensorMatrix, 1);
+defineNamedTemplateTypeNameAndDebug(blockTensorMatrix, 1);
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/foam/matrices/lduMatrix/lduMatrix/lduMatrixUpdateMatrixInterfaces.C

@@ -54,7 +54,10 @@ void Foam::lduMatrix::initMatrixInterfaces
                     *this,
                     coupleCoeffs[interfaceI],
                     cmpt,
-                    static_cast<Pstream::commsTypes>(Pstream::defaultCommsType()),
+                    static_cast<Pstream::commsTypes>
+                    (
+                        Pstream::defaultCommsType()
+                    ),
                     switchToLhs
                 );
             }
@@ -68,8 +71,8 @@ void Foam::lduMatrix::initMatrixInterfaces
         // beyond the end of the schedule which only handles "normal" patches
         for
         (
-            label interfaceI=patchSchedule.size()/2;
+            label interfaceI = patchSchedule.size()/2;
-            interfaceI<interfaces.size();
+            interfaceI < interfaces.size();
             interfaceI++
         )
         {
@@ -185,8 +188,8 @@ void Foam::lduMatrix::updateMatrixInterfaces
         // beyond the end of the schedule which only handles "normal" patches
         for
         (
-            label interfaceI=patchSchedule.size()/2;
+            label interfaceI = patchSchedule.size()/2;
-            interfaceI<interfaces.size();
+            interfaceI < interfaces.size();
             interfaceI++
         )
         {

src/foam/meshes/polyMesh/polyMeshInitMesh.C

@@ -223,9 +223,11 @@ void Foam::polyMesh::initMesh()
 
     string meshInfo =
         "nPoints: " + Foam::name(nPoints())
-      + " nCells: " + Foam::name(this->nCells())
+      + " nAllPoints: " + Foam::name(allPoints().size())
+      + " nInternalFaces: " + Foam::name(nInternalFaces())
       + " nFaces: " + Foam::name(nFaces())
-      + " nInternalFaces: " + Foam::name(nInternalFaces());
+      + " nAllFaces: " + Foam::name(allFaces().size())
+      + " nCells: " + Foam::name(this->nCells());
 
     owner_.note() = meshInfo;
     neighbour_.note() = meshInfo;

src/foam/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.C

@@ -902,7 +902,11 @@ void Foam::ggiPolyPatch::calcTransforms() const
 
     if (debug > 1 && master())
     {
-        if (patchToPatch().uncoveredMasterFaces().size() > 0)
+        if
+        (
+            !empty()
+         && patchToPatch().uncoveredMasterFaces().size() > 0
+        )
         {
             // Write uncovered master faces
             Info<< "Writing uncovered master faces for patch "
@@ -925,7 +929,11 @@ void Foam::ggiPolyPatch::calcTransforms() const
             );
         }
 
-        if (patchToPatch().uncoveredSlaveFaces().size() > 0)
+        if
+        (
+            !shadow().empty()
+         && patchToPatch().uncoveredSlaveFaces().size() > 0
+        )
         {
             // Write uncovered master faces
             Info<< "Writing uncovered shadow faces for patch "
@@ -935,7 +943,11 @@ void Foam::ggiPolyPatch::calcTransforms() const
 
             fileName fvPath(mesh.time().path()/"VTK");
             mkDir(fvPath);
-
+            Pout<< "shadow().localFaces(): " << shadow().localFaces().size()
+                << " patchToPatch().uncoveredSlaveFaces().size(): "
+                << patchToPatch().uncoveredSlaveFaces().size()
+                << " shadow().localPoints(): " << shadow().localPoints().size()
+                << endl;
             indirectPrimitivePatch::writeVTK
             (
                 fvPath/fileName("uncoveredGgiFaces" + shadowName()),
@@ -953,11 +965,17 @@ void Foam::ggiPolyPatch::calcTransforms() const
         {
             if
             (
-                patchToPatch().uncoveredMasterFaces().size() > 0
+                (
-             || patchToPatch().uncoveredSlaveFaces().size() > 0
+                    patchToPatch().uncoveredMasterFaces().size() > 0
+                    && !empty()
+                )
+             || (
+                    !shadow().empty()
+                 && patchToPatch().uncoveredSlaveFaces().size() > 0
+                )
             )
             {
-                FatalErrorIn("label ggiPolyPatch::shadowIndex() const")
+                FatalErrorIn("label ggiPolyPatch::calcTransforms() const")
                     << "ggi patch " << name() << " with shadow "
                     << shadowName() << " has "
                     << patchToPatch().uncoveredMasterFaces().size()
@@ -969,14 +987,14 @@ void Foam::ggiPolyPatch::calcTransforms() const
         }
         else
         {
-            InfoIn("label ggiPolyPatch::shadowIndex() const")
+            InfoIn("label ggiPolyPatch::calcTransforms() const")
                 << "ggi patch " << name() << " with shadow "
                 << shadowName() << " has "
                 << patchToPatch().uncoveredMasterFaces().size()
                 << " uncovered master faces and "
                 << patchToPatch().uncoveredSlaveFaces().size()
                 << " uncovered slave faces.  Bridging is switched on. "
-                << abort(FatalError);
+                << endl;
         }
     }
 }

src/foam/meshes/polyMesh/polyPatches/derived/cohesive/cohesivePolyPatch.C

@@ -119,7 +119,7 @@ bool cohesivePolyPatch::order
     labelList& rotation
 ) const
 {
-    // Grab faceMap from polyTopoChanger
+    // Grab faceMap
     SortableList<label> topoFaceMap(faceMap);
 
     faceMap.setSize(pp.size());

src/foam/primitives/strings/word/wordI.H

@@ -25,8 +25,6 @@ License
 
 #include <cctype>
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 inline void Foam::word::stripInvalid()

src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatch/immersedBoundaryFvPatch.C

@@ -1127,7 +1127,7 @@ void Foam::immersedBoundaryFvPatch::makeIbCellCells() const
     // Note: the algorithm is originally written with inward-facing normals
     // and subsequently changed: IB surface normals point outwards
     // HJ, 21/May/2012
-    const vectorField& ibn = ibNormals();
+//     const vectorField& ibn = ibNormals();
 
     forAll (cellCells, cellI)
     {
@@ -1152,8 +1152,8 @@ void Foam::immersedBoundaryFvPatch::makeIbCellCells() const
             // Collect the cells within rM of the fitting cell
             if (r <= rM[cellI])
             {
-                scalar angleLimit =
+//                 scalar angleLimit =
-                    -Foam::cos(angleFactor_()*mathematicalConstant::pi/180);
+//                     -Foam::cos(angleFactor_()*mathematicalConstant::pi/180);
 
                 vector dir = (C[curCell] - ibp[cellI]);
                 dir /= mag(dir) + SMALL;

src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatchField/immersedBoundaryFvPatchField.C

@@ -39,7 +39,8 @@ template<class Type>
 const Foam::debug::optimisationSwitch
 immersedBoundaryFvPatchField<Type>::nBcIter_
 (
-    debug::optimisationSwitch("immersedBoundaryNBCIter", 5)
+    "immersedBoundaryNBCIter",
+    5
 );
 
 
@@ -251,7 +252,7 @@ immersedBoundaryFvPatchField<Type>::imposeDirichletCondition() const
 
         reduce(maxError, maxOp<scalar>());
     }
-    while (maxError > bcTolerance_() && counter < nBcIter_());
+    while (maxError > bcTolerance_ && counter < nBcIter_);
 
     if (counter == nBcIter_() && debug)
     {

src/lduSolvers/amg/coarseAmgLevel.C

@@ -209,7 +209,7 @@ void Foam::coarseAmgLevel::solve
     }
 
     // Switch of debug in top-level direct solve
-    debug::debugSwitch oldDebug = lduMatrix::debug;
+    label oldDebug = lduMatrix::debug();
 
     if (matrixPtr_->matrix().symmetric())
     {

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.C

@@ -158,9 +158,9 @@ Foam::autoHexMeshDriver::autoHexMeshDriver
     {
         meshRefinement::debug = debug_;
         autoHexMeshDriver::debug = debug_;
-        autoRefineDriver::debug = debug;
+        autoRefineDriver::debug = debug_;
-        autoSnapDriver::debug = debug;
+        autoSnapDriver::debug = debug_;
-        autoLayerDriver::debug = debug;
+        autoLayerDriver::debug = debug_;
     }
 
     refinementParameters refineParams(dict, 1);

src/thermophysicalModels/basic/psiThermo/ePsiThermo/ePsiThermo.C

@@ -162,7 +162,7 @@ void Foam::ePsiThermo<MixtureType>::correct()
         Info<< "entering ePsiThermo<MixtureType>::correct()" << endl;
     }
 
-    // force the saving of the old-time values
+    // Force the saving of the old-time values
     this->psi_.oldTime();
 
     calculate();

src/turbulenceModels/compressible/RAS/kOmegaSST/kOmegaSST.C

@@ -85,6 +85,31 @@ tmp<volScalarField> kOmegaSST::F2() const
 }
 
 
+tmp<volScalarField> kOmegaSST::F3() const
+{
+    tmp<volScalarField> arg3 = min
+    (
+        150*(mu()/rho_)/(omega_*sqr(y_)),
+        scalar(10)
+    );
+
+    return 1 - tanh(pow4(arg3));
+}
+
+
+tmp<volScalarField> kOmegaSST::F23() const
+{
+    tmp<volScalarField> f23(F2());
+
+    if (F3_)
+    {
+        f23() *= F3();
+    }
+
+    return f23;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 kOmegaSST::kOmegaSST
@@ -103,7 +128,7 @@ kOmegaSST::kOmegaSST
         (
             "alphaK1",
             coeffDict_,
-            0.85034
+            0.85
         )
     ),
     alphaK2_
@@ -130,7 +155,7 @@ kOmegaSST::kOmegaSST
         (
             "alphaOmega2",
             coeffDict_,
-            0.85616
+            0.856
         )
     ),
     Prt_
@@ -148,7 +173,7 @@ kOmegaSST::kOmegaSST
         (
             "gamma1",
             coeffDict_,
-            0.5532
+            5.0/9.0
         )
     ),
     gamma2_
@@ -157,7 +182,7 @@ kOmegaSST::kOmegaSST
         (
             "gamma2",
             coeffDict_,
-            0.4403
+            0.44
         )
     ),
     beta1_
@@ -196,6 +221,15 @@ kOmegaSST::kOmegaSST
             0.31
         )
     ),
+    b1_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "b1",
+            coeffDict_,
+            1.0
+        )
+    ),
     c1_
     (
         dimensioned<scalar>::lookupOrAddToDict
@@ -205,6 +239,15 @@ kOmegaSST::kOmegaSST
             10.0
         )
     ),
+    F3_
+    (
+        Switch::lookupOrAddToDict
+        (
+            "F3",
+            coeffDict_,
+            false
+        )
+    ),
 
     y_(mesh_),
 
@@ -257,6 +300,7 @@ kOmegaSST::kOmegaSST
         autoCreateAlphat("alphat", mesh_)
     )
 {
+    bound(k_, k0_);
     bound(omega_, omega0_);
 
     mut_ =
@@ -265,7 +309,7 @@ kOmegaSST::kOmegaSST
         max
         (
             a1_*omega_,
-            F2()*sqrt(2*magSqr(symm(fvc::grad(U_))))
+            b1_*F23()*sqrt(2.0)*mag(symm(fvc::grad(U_)))
         )
     );
     mut_ = min(mut_, muRatio()*mu());
@@ -345,7 +389,9 @@ bool kOmegaSST::read()
         beta2_.readIfPresent(coeffDict());
         betaStar_.readIfPresent(coeffDict());
         a1_.readIfPresent(coeffDict());
+        b1_.readIfPresent(coeffDict());
         c1_.readIfPresent(coeffDict());
+        F3_.readIfPresent("F3", coeffDict());
 
         return true;
     }
@@ -371,7 +417,7 @@ void kOmegaSST::correct()
         // Re-calculate viscosity
         mut_ =
             a1_*rho_*k_
-           /max(a1_*omega_, F2()*sqrt(2*magSqr(symm(fvc::grad(U_)))));
+           /max(a1_*omega_, F2()*sqrt(2.0)*mag(symm(fvc::grad(U_))));
         mut_ = min(mut_, muRatio()*mu());
         mut_.correctBoundaryConditions();
 
@@ -384,7 +430,7 @@ void kOmegaSST::correct()
 
     RASModel::correct();
 
-    volScalarField divU = fvc::div(phi_/fvc::interpolate(rho_));
+    volScalarField divU(fvc::div(phi_/fvc::interpolate(rho_)));
 
     if (mesh_.changing())
     {
@@ -397,19 +443,21 @@ void kOmegaSST::correct()
     }
 
     tmp<volTensorField> tgradU = fvc::grad(U_);
-    volScalarField S2 = magSqr(symm(tgradU()));
+    volScalarField S2(2*magSqr(symm(tgradU())));
-    volScalarField GbyMu = (tgradU() && dev(twoSymm(tgradU())));
+    volScalarField GbyMu((tgradU() && dev(twoSymm(tgradU()))));
     volScalarField G("RASModel::G", mut_*GbyMu);
     tgradU.clear();
 
     // Update omega and G at the wall
     omega_.boundaryField().updateCoeffs();
 
-    volScalarField CDkOmega =
+    volScalarField CDkOmega
-        (2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_;
+    (
+        (2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_
+    );
 
-    volScalarField F1 = this->F1(CDkOmega);
+    volScalarField F1(this->F1(CDkOmega));
-    volScalarField rhoGammaF1 = rho_*gamma(F1);
+    volScalarField rhoGammaF1(rho_*gamma(F1));
 
     // Turbulent frequency equation
     tmp<fvScalarMatrix> omegaEqn
@@ -418,7 +466,12 @@ void kOmegaSST::correct()
       + fvm::div(phi_, omega_)
       - fvm::laplacian(DomegaEff(F1), omega_)
      ==
-        rhoGammaF1*GbyMu
+        rhoGammaF1
+       *min
+        (
+            GbyMu,
+            (c1_/a1_)*betaStar_*omega_*max(a1_*omega_, b1_*F23()*sqrt(S2))
+        )
       - fvm::SuSp((2.0/3.0)*rhoGammaF1*divU, omega_)
       - fvm::Sp(rho_*beta(F1)*omega_, omega_)
       - fvm::SuSp
@@ -455,7 +508,8 @@ void kOmegaSST::correct()
 
 
     // Re-calculate viscosity
-    mut_ = a1_*rho_*k_/max(a1_*omega_, F2()*sqrt(2*S2));
+    // Fixed sqrt(2) error.  HJ, 10/Jun/2015
+    mut_ = a1_*rho_*k_/max(a1_*omega_, b1_*F23()*sqrt(S2));
     mut_ = min(mut_, muRatio()*mu());
     mut_.correctBoundaryConditions();
 

src/turbulenceModels/compressible/RAS/kOmegaSST/kOmegaSST.H

@@ -27,12 +27,35 @@ Class
 Description
     Implementation of the k-omega-SST turbulence model for compressible flows.
 
-    Turbulence model described in:
+    Turbulence model described in
     @verbatim
-        Menter, F., Esch, T.
+        Menter, F., Esch, T.,
-        "Elements of Industrial Heat Transfer Prediction"
+        "Elements of Industrial Heat Transfer Prediction",
         16th Brazilian Congress of Mechanical Engineering (COBEM),
-        Nov. 2001
+        Nov. 2001.
+    @endverbatim
+
+    with updated coefficients from
+    @verbatim
+        Menter, F. R., Kuntz, M., and Langtry, R.,
+        "Ten Years of Industrial Experience with the SST Turbulence Model",
+        Turbulence, Heat and Mass Transfer 4, 2003,
+        pp. 625 - 632.
+    @endverbatim
+
+    but with the consistent production terms from the 2001 paper as form in the
+    2003 paper is a typo, see
+    @verbatim
+        http://turbmodels.larc.nasa.gov/sst.html
+    @endverbatim
+
+    and the addition of the optional F3 term for rough walls from
+    @verbatim
+        Hellsten, A.
+        "Some Improvements in Menter’s k-omega-SST turbulence model"
+        29th AIAA Fluid Dynamics Conference,
+        AIAA-98-2554,
+        June 1998.
     @endverbatim
 
     Note that this implementation is written in terms of alpha diffusion
@@ -44,29 +67,24 @@ Description
     Also note that the error in the last term of equation (2) relating to
     sigma has been corrected.
 
-    Wall-functions are applied in this implementation by using equations (14)
-    to specify the near-wall omega as appropriate.
-
-    The blending functions (15) and (16) are not currently used because of the
-    uncertainty in their origin, range of applicability and that is y+ becomes
-    sufficiently small blending u_tau in this manner clearly becomes nonsense.
-
     The default model coefficients correspond to the following:
     @verbatim
         kOmegaSSTCoeffs
         {
-            alphaK1     0.85034;
+            alphaK1     0.85;
             alphaK2     1.0;
             alphaOmega1 0.5;
-            alphaOmega2 0.85616;
+            alphaOmega2 0.856;
-            Prt         1.0;    // only for compressible
             beta1       0.075;
             beta2       0.0828;
             betaStar    0.09;
-            gamma1      0.5532;
+            gamma1      5/9;
-            gamma2      0.4403;
+            gamma2      0.44;
             a1          0.31;
+            b1          1.0;
             c1          10.0;
+            F3          no;
+            Prt         1.0;
         }
     @endverbatim
 
@@ -122,8 +140,11 @@ class kOmegaSST
             dimensionedScalar betaStar_;
 
             dimensionedScalar a1_;
+            dimensionedScalar b1_;
             dimensionedScalar c1_;
 
+            Switch F3_;
+
 
         //- Wall distance
         //  Note: different to wall distance in parent RASModel
@@ -141,6 +162,8 @@ class kOmegaSST
 
         tmp<volScalarField> F1(const volScalarField& CDkOmega) const;
         tmp<volScalarField> F2() const;
+        tmp<volScalarField> F3() const;
+        tmp<volScalarField> F23() const;
 
         tmp<volScalarField> blend
         (

src/turbulenceModels/incompressible/RAS/coupledKEpsilon/coupledKEpsilon.C

@@ -302,6 +302,16 @@ void coupledKEpsilon::correct()
         keEqn.insertEquation(0, kEqn);
     }
 
+    // Add coupling term: C1*Cmu*(symm(grad(U))) k but with wall function
+    // corrections: must be calculated from G.  HJ, 27/Apr/2015
+
+    // Add coupling term: epsilon source depends on k
+    // k, e sink terms cannot be changed because of boundedness
+    keEqn.insertEquationCoupling
+    (
+        1, 0, -C1_*G*epsilon_/sqr(k_)
+    );
+
     // Update source coupling: coupling terms eliminated from source
     keEqn.updateSourceCoupling();
 

src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C

@@ -86,6 +86,31 @@ tmp<volScalarField> kOmegaSST::F2() const
 }
 
 
+tmp<volScalarField> kOmegaSST::F3() const
+{
+    tmp<volScalarField> arg3 = min
+    (
+        150*nu()/(omega_*sqr(y_)),
+        scalar(10)
+    );
+
+    return 1 - tanh(pow4(arg3));
+}
+
+
+tmp<volScalarField> kOmegaSST::F23() const
+{
+    tmp<volScalarField> f23(F2());
+
+    if (F3_)
+    {
+        f23() *= F3();
+    }
+
+    return f23;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 kOmegaSST::kOmegaSST
@@ -103,7 +128,7 @@ kOmegaSST::kOmegaSST
         (
             "alphaK1",
             coeffDict_,
-            0.85034
+            0.85
         )
     ),
     alphaK2_
@@ -130,7 +155,7 @@ kOmegaSST::kOmegaSST
         (
             "alphaOmega2",
             coeffDict_,
-            0.85616
+            0.856
         )
     ),
     gamma1_
@@ -139,7 +164,7 @@ kOmegaSST::kOmegaSST
         (
             "gamma1",
             coeffDict_,
-            0.5532
+            5.0/9.0
         )
     ),
     gamma2_
@@ -148,7 +173,7 @@ kOmegaSST::kOmegaSST
         (
             "gamma2",
             coeffDict_,
-            0.4403
+            0.44
         )
     ),
     beta1_
@@ -187,6 +212,15 @@ kOmegaSST::kOmegaSST
             0.31
         )
     ),
+    b1_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "b1",
+            coeffDict_,
+            1.0
+        )
+    ),
     c1_
     (
         dimensioned<scalar>::lookupOrAddToDict
@@ -196,6 +230,15 @@ kOmegaSST::kOmegaSST
             10.0
         )
     ),
+    F3_
+    (
+        Switch::lookupOrAddToDict
+        (
+            "F3",
+            coeffDict_,
+            false
+        )
+    ),
 
     y_(mesh_),
 
@@ -236,10 +279,18 @@ kOmegaSST::kOmegaSST
         autoCreateNut("nut", mesh_, U_.db())
     )
 {
+    bound(k_, k0_);
     bound(omega_, omega0_);
 
-    nut_ = a1_*k_/max(a1_*omega_, F2()*sqrt(2.0)*mag(symm(fvc::grad(U_))));
+    nut_ =
-    nut_ = min(nut_, nuRatio()*nu());
+    (
+        a1_*k_/
+        max
+        (
+            a1_*omega_,
+            b1_*F23()*sqrt(2.0)*mag(symm(fvc::grad(U_)))
+        )
+    );
     nut_.correctBoundaryConditions();
 
     printCoeffs();
@@ -313,7 +364,9 @@ bool kOmegaSST::read()
         beta2_.readIfPresent(coeffDict());
         betaStar_.readIfPresent(coeffDict());
         a1_.readIfPresent(coeffDict());
+        b1_.readIfPresent(coeffDict());
         c1_.readIfPresent(coeffDict());
+        F3_.readIfPresent("F3", coeffDict());
 
         return true;
     }
@@ -346,16 +399,18 @@ void kOmegaSST::correct()
         y_.correct();
     }
 
-    volScalarField S2 = magSqr(symm(fvc::grad(U_)));
+    const volScalarField S2(2*magSqr(symm(fvc::grad(U_))));
     volScalarField G("RASModel::G", nut_*2*S2);
 
     // Update omega and G at the wall
     omega_.boundaryField().updateCoeffs();
 
-    volScalarField CDkOmega =
+    const volScalarField CDkOmega
-        (2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_;
+    (
+        (2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_
+    );
 
-    volScalarField F1 = this->F1(CDkOmega);
+    const volScalarField F1(this->F1(CDkOmega));
 
     // Turbulent frequency equation
     tmp<fvScalarMatrix> omegaEqn
@@ -365,7 +420,8 @@ void kOmegaSST::correct()
       + fvm::SuSp(-fvc::div(phi_), omega_)
       - fvm::laplacian(DomegaEff(F1), omega_)
      ==
-        gamma(F1)*2*S2
+        gamma(F1)
+       *min(S2, (c1_/a1_)*betaStar_*omega_*max(a1_*omega_, b1_*F23()*sqrt(S2)))
       - fvm::Sp(beta(F1)*omega_, omega_)
       - fvm::SuSp
         (
@@ -401,7 +457,8 @@ void kOmegaSST::correct()
 
 
     // Re-calculate viscosity
-    nut_ = a1_*k_/max(a1_*omega_, F2()*sqrt(2*S2));
+    // Fixed sqrt(2) error.  HJ, 10/Jun/2015
+    nut_ = a1_*k_/max(a1_*omega_, b1_*F23()*sqrt(S2));
     nut_ = min(nut_, nuRatio()*nu());
     nut_.correctBoundaryConditions();
 }

src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.H

@@ -33,9 +33,32 @@ Description
         Menter, F., Esch, T.
         "Elements of Industrial Heat Transfer Prediction"
         16th Brazilian Congress of Mechanical Engineering (COBEM),
-        Nov. 2001
+        Nov. 2001.
     @endverbatim
 
+    with updated coefficients from
+    @verbatim
+        Menter, F. R., Kuntz, M., and Langtry, R.,
+        "Ten Years of Industrial Experience with the SST Turbulence Model",
+        Turbulence, Heat and Mass Transfer 4, 2003,
+        pp. 625 - 632.
+    @endverbatim
+
+    but with the consistent production terms from the 2001 paper as form in the
+    2003 paper is a typo, see
+    @verbatim
+        http://turbmodels.larc.nasa.gov/sst.html
+    @endverbatim
+
+    and the addition of the optional F3 term for rough walls from
+    \verbatim
+        Hellsten, A.
+        "Some Improvements in Menter’s k-omega-SST turbulence model"
+        29th AIAA Fluid Dynamics Conference,
+        AIAA-98-2554,
+        June 1998.
+    \endverbatim
+
     Note that this implementation is written in terms of alpha diffusion
     coefficients rather than the more traditional sigma (alpha = 1/sigma) so
     that the blending can be applied to all coefficuients in a consistent
@@ -45,28 +68,23 @@ Description
     Also note that the error in the last term of equation (2) relating to
     sigma has been corrected.
 
-    Wall-functions are applied in this implementation by using equations (14)
-    to specify the near-wall omega as appropriate.
-
-    The blending functions (15) and (16) are not currently used because of the
-    uncertainty in their origin, range of applicability and that is y+ becomes
-    sufficiently small blending u_tau in this manner clearly becomes nonsense.
-
     The default model coefficients correspond to the following:
     @verbatim
         kOmegaSSTCoeffs
         {
-            alphaK1     0.85034;
+            alphaK1     0.85;
             alphaK2     1.0;
             alphaOmega1 0.5;
-            alphaOmega2 0.85616;
+            alphaOmega2 0.856;
             beta1       0.075;
             beta2       0.0828;
             betaStar    0.09;
-            gamma1      0.5532;
+            gamma1      5/9;
-            gamma2      0.4403;
+            gamma2      0.44;
             a1          0.31;
+            b1          1.0;
             c1          10.0;
+            F3          no;
         }
     @endverbatim
 
@@ -116,8 +134,12 @@ class kOmegaSST
             dimensionedScalar betaStar_;
 
             dimensionedScalar a1_;
+            dimensionedScalar b1_;
             dimensionedScalar c1_;
 
+            Switch F3_;
+
+
         //- Wall distance field
         //  Note: different to wall distance in parent RASModel
         wallDist y_;
@@ -133,6 +155,8 @@ class kOmegaSST
 
         tmp<volScalarField> F1(const volScalarField& CDkOmega) const;
         tmp<volScalarField> F2() const;
+        tmp<volScalarField> F3() const;
+        tmp<volScalarField> F23() const;
 
         tmp<volScalarField> blend
         (
@@ -195,9 +219,8 @@ public:
 
 
     //- Destructor
-
+    virtual ~kOmegaSST()
-        virtual ~kOmegaSST()
+    {}
-        {}
 
 
     // Member Functions

tutorials/finiteArea/surfactantFoam/sphereTransport/sphereTransport/system/controlDict

@@ -44,4 +44,3 @@ timePrecision       6;
 runTimeModifiable   yes;
 
 // ************************************************************************* //
-

tutorials/incompressible/icoFoam/cavity/constant/polyMesh/boundary

@@ -1,9 +1,9 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | foam-extend: Open Source CFD                    |
-|  \\    /   O peration     | Version:     3.2                                |
+|  \\    /   O peration     | Version:     3.1                                |
 |   \\  /    A nd           | Web:         http://www.foam-extend.org         |
-|    \\/     M anipulation  |                                                 |
+|    \\/     M anipulation  | For copyright notice see file Copyright         |
 \*---------------------------------------------------------------------------*/
 FoamFile
 {

tutorials/incompressible/pUCoupledFoam/cavity/constant/polyMesh/boundary

@@ -1,9 +1,9 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | foam-extend: Open Source CFD                    |
-|  \\    /   O peration     | Version:     3.2                                |
+|  \\    /   O peration     | Version:     3.1                                |
 |   \\  /    A nd           | Web:         http://www.foam-extend.org         |
-|    \\/     M anipulation  |                                                 |
+|    \\/     M anipulation  | For copyright notice see file Copyright         |
 \*---------------------------------------------------------------------------*/
 FoamFile
 {

tutorials/incompressible/simpleFoam/pitzDaily/0/epsilon

@@ -3,7 +3,7 @@
 | \\      /  F ield         | foam-extend: Open Source CFD                    |
 |  \\    /   O peration     | Version:     3.2                                |
 |   \\  /    A nd           | Web:         http://www.foam-extend.org         |
-|    \\/     M anipulation  |                                                 |
+|    \\/     M anipulation  | For copyright notice see file Copyright         |
 \*---------------------------------------------------------------------------*/
 FoamFile
 {

tutorials/incompressible/simpleFoam/pitzDaily/0/k

@@ -3,7 +3,7 @@
 | \\      /  F ield         | foam-extend: Open Source CFD                    |
 |  \\    /   O peration     | Version:     3.2                                |
 |   \\  /    A nd           | Web:         http://www.foam-extend.org         |
-|    \\/     M anipulation  |                                                 |
+|    \\/     M anipulation  | For copyright notice see file Copyright         |
 \*---------------------------------------------------------------------------*/
 FoamFile
 {
@@ -46,4 +46,5 @@ boundaryField
     }
 }
 
+
 // ************************************************************************* //

tutorials/incompressible/simpleFoam/pitzDaily/0/nut

@@ -0,0 +1,57 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | foam-extend: Open Source CFD                    |
+|  \\    /   O peration     | Version:     3.1                                |
+|   \\  /    A nd           | Web:         http://www.foam-extend.org         |
+|    \\/     M anipulation  | For copyright notice see file Copyright         |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      nut;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    outlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    upperWall
+    {
+        type            nutWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 0;
+    }
+    lowerWall
+    {
+        type            nutWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 0;
+    }
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/incompressible/simpleFoam/pitzDaily/constant/polyMesh/boundary

@@ -0,0 +1,52 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | foam-extend: Open Source CFD                    |
+|  \\    /   O peration     | Version:     3.1                                |
+|   \\  /    A nd           | Web:         http://www.foam-extend.org         |
+|    \\/     M anipulation  | For copyright notice see file Copyright         |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       polyBoundaryMesh;
+    location    "constant/polyMesh";
+    object      boundary;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+5
+(
+    inlet
+    {
+        type            patch;
+        nFaces          30;
+        startFace       24170;
+    }
+    outlet
+    {
+        type            patch;
+        nFaces          57;
+        startFace       24200;
+    }
+    upperWall
+    {
+        type            wall;
+        nFaces          223;
+        startFace       24257;
+    }
+    lowerWall
+    {
+        type            wall;
+        nFaces          250;
+        startFace       24480;
+    }
+    frontAndBack
+    {
+        type            empty;
+        nFaces          24450;
+        startFace       24730;
+    }
+)
+
+// ************************************************************************* //