From fceb3b566e8240f2b5957fc1e7721c7e4a55db40 Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Wed, 3 Jun 2015 09:15:02 +0100
Subject: [PATCH 01/24] Formatting

---
 .../gradSchemes/leastSquaresGrad/leastSquaresGrad.C       | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresGrad.C b/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresGrad.C
index 928187873..7f3bdad0f 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresGrad.C
+++ b/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]

From 09bd0e1c70b2d5a722d172c6abf87f394375eb5c Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Wed, 3 Jun 2015 09:17:28 +0100
Subject: [PATCH 02/24] Formatting

---
 src/foam/primitives/strings/word/wordI.H | 2 --
 1 file changed, 2 deletions(-)

diff --git a/src/foam/primitives/strings/word/wordI.H b/src/foam/primitives/strings/word/wordI.H
index af61de9a0..c2dd8bab7 100644
--- a/src/foam/primitives/strings/word/wordI.H
+++ b/src/foam/primitives/strings/word/wordI.H
@@ -25,8 +25,6 @@ License
 
 #include <cctype>
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 inline void Foam::word::stripInvalid()

From 9b576b99164a8a2702a24ba413b8b9e7c7afd5b9 Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Wed, 3 Jun 2015 09:18:00 +0100
Subject: [PATCH 03/24] Formatting

---
 .../polyMesh/polyPatches/derived/cohesive/cohesivePolyPatch.C   | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/foam/meshes/polyMesh/polyPatches/derived/cohesive/cohesivePolyPatch.C b/src/foam/meshes/polyMesh/polyPatches/derived/cohesive/cohesivePolyPatch.C
index d07887e8a..e9ed148e2 100644
--- a/src/foam/meshes/polyMesh/polyPatches/derived/cohesive/cohesivePolyPatch.C
+++ b/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());

From 15e13ead7f378b16da7dcbccc784fa1008cc0474 Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Wed, 3 Jun 2015 09:18:16 +0100
Subject: [PATCH 04/24] Better mesh info data

---
 src/foam/meshes/polyMesh/polyMeshInitMesh.C | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/foam/meshes/polyMesh/polyMeshInitMesh.C b/src/foam/meshes/polyMesh/polyMeshInitMesh.C
index 46478dc78..2a4564ec9 100644
--- a/src/foam/meshes/polyMesh/polyMeshInitMesh.C
+++ b/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;

From d5e847baf5f63e9a5ca3dc2b2fdda558107e9179 Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Wed, 3 Jun 2015 09:19:27 +0100
Subject: [PATCH 05/24] Corrected handling of mixing type

---
 .../mixingPlane/mixingPlaneFvPatchField.C     | 70 ++++++++++++++++---
 1 file changed, 62 insertions(+), 8 deletions(-)

diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/mixingPlane/mixingPlaneFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/constraint/mixingPlane/mixingPlaneFvPatchField.C
index 59196778f..be9da33d9 100644
--- a/src/finiteVolume/fields/fvPatchFields/constraint/mixingPlane/mixingPlaneFvPatchField.C
+++ b/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
-    readMixingType();
+    // If mixing type is unknown, read it
+    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
-    readMixingType();
+    // If mixing type is unknown, read it
+    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
-    readMixingType();
+    // If mixing type is unknown, read it
+    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

From 37492e45f4db62545566a20f064f30f9e0e66251 Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Thu, 4 Jun 2015 17:53:51 +0100
Subject: [PATCH 06/24] Changed default debug level of block matrix

---
 .../blockLduMatrix/BlockLduMatrix/blockLduMatrices.C  | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/blockLduMatrices.C b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/blockLduMatrices.C
index 6824f8c62..84e367a80 100644
--- a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/blockLduMatrices.C
+++ b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/blockLduMatrices.C
@@ -38,11 +38,12 @@ namespace Foam
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-defineNamedTemplateTypeNameAndDebug(blockScalarMatrix, 0);
-defineNamedTemplateTypeNameAndDebug(blockVectorMatrix, 0);
-defineNamedTemplateTypeNameAndDebug(blockSphericalTensorMatrix, 0);
-defineNamedTemplateTypeNameAndDebug(blockSymmTensorMatrix, 0);
-defineNamedTemplateTypeNameAndDebug(blockTensorMatrix, 0);
+// Note: Matrix debug level defaults to 1.  HJ, 4/Jun/2015
+defineNamedTemplateTypeNameAndDebug(blockScalarMatrix, 1);
+defineNamedTemplateTypeNameAndDebug(blockVectorMatrix, 1);
+defineNamedTemplateTypeNameAndDebug(blockSphericalTensorMatrix, 1);
+defineNamedTemplateTypeNameAndDebug(blockSymmTensorMatrix, 1);
+defineNamedTemplateTypeNameAndDebug(blockTensorMatrix, 1);
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

From bf0f746ed555d0019dbc624c058f365fb34a2d84 Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Thu, 4 Jun 2015 17:54:17 +0100
Subject: [PATCH 07/24] Fixed debug treatment - cannot copy debug switch

---
 src/foam/matrices/blockLduMatrix/BlockAmg/coarseBlockAmgLevel.C | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/foam/matrices/blockLduMatrix/BlockAmg/coarseBlockAmgLevel.C b/src/foam/matrices/blockLduMatrix/BlockAmg/coarseBlockAmgLevel.C
index 197ba84ed..6e9dd3b40 100644
--- a/src/foam/matrices/blockLduMatrix/BlockAmg/coarseBlockAmgLevel.C
+++ b/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)
     {

From 48e45fb1fb535e95baeca8c994f31ac4a872d241 Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Thu, 4 Jun 2015 18:02:34 +0100
Subject: [PATCH 08/24] Formatting

---
 .../basic/psiThermo/ePsiThermo/ePsiThermo.C                     | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/thermophysicalModels/basic/psiThermo/ePsiThermo/ePsiThermo.C b/src/thermophysicalModels/basic/psiThermo/ePsiThermo/ePsiThermo.C
index 07bc449cf..0a9ab5438 100644
--- a/src/thermophysicalModels/basic/psiThermo/ePsiThermo/ePsiThermo.C
+++ b/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();

From 50dc19899e392fcd4603e21d32cffbe1af2bf523 Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Thu, 4 Jun 2015 18:03:34 +0100
Subject: [PATCH 09/24] Alternative form of k-e coupling

---
 .../RAS/coupledKEpsilon/coupledKEpsilon.C              | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/src/turbulenceModels/incompressible/RAS/coupledKEpsilon/coupledKEpsilon.C b/src/turbulenceModels/incompressible/RAS/coupledKEpsilon/coupledKEpsilon.C
index c5842af24..613ed4ea5 100644
--- a/src/turbulenceModels/incompressible/RAS/coupledKEpsilon/coupledKEpsilon.C
+++ b/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();
 

From 7f92fbcae4614711cdaa6eb323334c9a9ddddf22 Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Thu, 4 Jun 2015 18:04:33 +0100
Subject: [PATCH 10/24] Prevent copying of switches: this causes double
 deletion error

---
 .../global/controlSwitches/constantsSwitch.H  |  9 +++
 .../global/controlSwitches/controlSwitches.C  | 59 +------------------
 .../global/controlSwitches/controlSwitches.H  | 23 ++++----
 src/foam/global/controlSwitches/debugSwitch.H |  9 +++
 src/foam/global/controlSwitches/infoSwitch.H  | 12 +++-
 .../controlSwitches/optimisationSwitch.H      |  9 +++
 src/foam/global/controlSwitches/safeBool.H    |  2 +-
 .../global/controlSwitches/tolerancesSwitch.H |  9 +++
 .../BlockLduMatrix/BlockLduMatrix.C           | 39 ++++++++----
 src/lduSolvers/amg/coarseAmgLevel.C           |  2 +-
 10 files changed, 89 insertions(+), 84 deletions(-)

diff --git a/src/foam/global/controlSwitches/constantsSwitch.H b/src/foam/global/controlSwitches/constantsSwitch.H
index 01fafa035..dd6bb7f03 100644
--- a/src/foam/global/controlSwitches/constantsSwitch.H
+++ b/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
diff --git a/src/foam/global/controlSwitches/controlSwitches.C b/src/foam/global/controlSwitches/controlSwitches.C
index f107c8030..de28ac95c 100644
--- a/src/foam/global/controlSwitches/controlSwitches.C
+++ b/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  * * * * * * * * * * * * * * //
-
-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 * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
 
 template<class T>
 void printControlSwitches
diff --git a/src/foam/global/controlSwitches/controlSwitches.H b/src/foam/global/controlSwitches/controlSwitches.H
index 06e05bd5e..e6f144fa2 100644
--- a/src/foam/global/controlSwitches/controlSwitches.H
+++ b/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&);
 
diff --git a/src/foam/global/controlSwitches/debugSwitch.H b/src/foam/global/controlSwitches/debugSwitch.H
index 671091afe..22c560b7a 100644
--- a/src/foam/global/controlSwitches/debugSwitch.H
+++ b/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
diff --git a/src/foam/global/controlSwitches/infoSwitch.H b/src/foam/global/controlSwitches/infoSwitch.H
index 1f9200cac..5ae398d93 100644
--- a/src/foam/global/controlSwitches/infoSwitch.H
+++ b/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
diff --git a/src/foam/global/controlSwitches/optimisationSwitch.H b/src/foam/global/controlSwitches/optimisationSwitch.H
index afe9dea39..a66ee7498 100644
--- a/src/foam/global/controlSwitches/optimisationSwitch.H
+++ b/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
diff --git a/src/foam/global/controlSwitches/safeBool.H b/src/foam/global/controlSwitches/safeBool.H
index 39e00a80e..1dd9eade0 100644
--- a/src/foam/global/controlSwitches/safeBool.H
+++ b/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;
     }
 
diff --git a/src/foam/global/controlSwitches/tolerancesSwitch.H b/src/foam/global/controlSwitches/tolerancesSwitch.H
index 001f9b061..5f5f22bed 100644
--- a/src/foam/global/controlSwitches/tolerancesSwitch.H
+++ b/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
diff --git a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrix.C b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrix.C
index dad0bcbd3..86cf42339 100644
--- a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrix.C
+++ b/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
-            (ldum.lduAddr().size()) << token::END_STATEMENT << nl;
+        os.writeKeyword("diagonal")
+            << 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
-            (ldum.lduAddr().lowerAddr().size()) <<  token::END_STATEMENT << nl;
+        os.writeKeyword("upper")
+            << 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
-            (ldum.lduAddr().lowerAddr().size()) <<  token::END_STATEMENT << nl;
+        os.writeKeyword("lower")
+            << 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>&");
 
diff --git a/src/lduSolvers/amg/coarseAmgLevel.C b/src/lduSolvers/amg/coarseAmgLevel.C
index 3e026cbe5..c4a9deda5 100644
--- a/src/lduSolvers/amg/coarseAmgLevel.C
+++ b/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())
     {

From 4ccbdb4a75138d4e8d1f84d6c4a28c4cb49cc31b Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Wed, 10 Jun 2015 11:42:39 +0100
Subject: [PATCH 11/24] Updated version of k-omega SST model

---
 .../compressible/RAS/kOmegaSST/kOmegaSST.C    | 84 +++++++++++++++----
 .../compressible/RAS/kOmegaSST/kOmegaSST.H    | 55 ++++++++----
 .../incompressible/RAS/kOmegaSST/kOmegaSST.C  | 81 +++++++++++++++---
 .../incompressible/RAS/kOmegaSST/kOmegaSST.H  | 53 ++++++++----
 4 files changed, 215 insertions(+), 58 deletions(-)

diff --git a/src/turbulenceModels/compressible/RAS/kOmegaSST/kOmegaSST.C b/src/turbulenceModels/compressible/RAS/kOmegaSST/kOmegaSST.C
index c274fd623..8149b4082 100644
--- a/src/turbulenceModels/compressible/RAS/kOmegaSST/kOmegaSST.C
+++ b/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 GbyMu = (tgradU() && dev(twoSymm(tgradU())));
+    volScalarField S2(2*magSqr(symm(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 =
-        (2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_;
+    volScalarField CDkOmega
+    (
+        (2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_
+    );
 
-    volScalarField F1 = this->F1(CDkOmega);
-    volScalarField rhoGammaF1 = rho_*gamma(F1);
+    volScalarField F1(this->F1(CDkOmega));
+    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();
 
diff --git a/src/turbulenceModels/compressible/RAS/kOmegaSST/kOmegaSST.H b/src/turbulenceModels/compressible/RAS/kOmegaSST/kOmegaSST.H
index 7a72030f7..13595376f 100644
--- a/src/turbulenceModels/compressible/RAS/kOmegaSST/kOmegaSST.H
+++ b/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.
-        "Elements of Industrial Heat Transfer Prediction"
+        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
@@ -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;
-            Prt         1.0;    // only for compressible
+            alphaOmega2 0.856;
             beta1       0.075;
             beta2       0.0828;
             betaStar    0.09;
-            gamma1      0.5532;
-            gamma2      0.4403;
+            gamma1      5/9;
+            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
         (
diff --git a/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C b/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C
index ee4f91da0..842479757 100644
--- a/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C
+++ b/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_ = min(nut_, nuRatio()*nu());
+    nut_ =
+    (
+        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 =
-        (2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_;
+    const volScalarField CDkOmega
+    (
+        (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();
 }
diff --git a/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.H b/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.H
index 4da8488e4..5382b5f98 100644
--- a/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.H
+++ b/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;
-            gamma2      0.4403;
+            gamma1      5/9;
+            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

From 867da248b24c120112a3da2488947f732306f8dc Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Wed, 10 Jun 2015 11:42:56 +0100
Subject: [PATCH 12/24] Better messaging

---
 .../patchConstrainedDecomp/patchConstrainedDecomp.C          | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/src/decompositionMethods/decompositionMethods/patchConstrainedDecomp/patchConstrainedDecomp.C b/src/decompositionMethods/decompositionMethods/patchConstrainedDecomp/patchConstrainedDecomp.C
index 2c997a2cb..5250f9fa3 100644
--- a/src/decompositionMethods/decompositionMethods/patchConstrainedDecomp/patchConstrainedDecomp.C
+++ b/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)

From a8bc92b12ae64de6d2c17db28cd83f173ad8c249 Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Wed, 10 Jun 2015 11:43:21 +0100
Subject: [PATCH 13/24] Changed default Pstream comms type to scheduled

---
 src/foam/db/IOstreams/Pstreams/Pstream.C | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/foam/db/IOstreams/Pstreams/Pstream.C b/src/foam/db/IOstreams/Pstreams/Pstream.C
index 9089eb347..f5045df3d 100644
--- a/src/foam/db/IOstreams/Pstreams/Pstream.C
+++ b/src/foam/db/IOstreams/Pstreams/Pstream.C
@@ -269,7 +269,8 @@ const Foam::debug::optimisationSwitch
 Foam::Pstream::defaultCommsType
 (
     "commsType",
-    "blocking",
+    "scheduled",
+//     "blocking",
     "blocking, nonBlocking, scheduled"
 );
 

From b6b47d10176fd16eb248adfdff393de4c52bacbf Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Wed, 10 Jun 2015 11:43:55 +0100
Subject: [PATCH 14/24] Porting bug fix

---
 .../immersedBoundaryFvPatch/immersedBoundaryFvPatch.C       | 6 +++---
 .../immersedBoundaryFvPatchField.C                          | 5 +++--
 2 files changed, 6 insertions(+), 5 deletions(-)

diff --git a/src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatch/immersedBoundaryFvPatch.C b/src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatch/immersedBoundaryFvPatch.C
index df0d5aee1..1a19cc8c2 100644
--- a/src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatch/immersedBoundaryFvPatch.C
+++ b/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 =
-                    -Foam::cos(angleFactor_()*mathematicalConstant::pi/180);
+//                 scalar angleLimit =
+//                     -Foam::cos(angleFactor_()*mathematicalConstant::pi/180);
 
                 vector dir = (C[curCell] - ibp[cellI]);
                 dir /= mag(dir) + SMALL;
diff --git a/src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatchField/immersedBoundaryFvPatchField.C b/src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatchField/immersedBoundaryFvPatchField.C
index 432b4ccc5..33778c5dc 100644
--- a/src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatchField/immersedBoundaryFvPatchField.C
+++ b/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)
     {

From 8470a791b436f21a33eb5558aff2ec1833000dd1 Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Wed, 10 Jun 2015 11:45:11 +0100
Subject: [PATCH 15/24] Copying correct debug

---
 .../autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.C       | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.C
index dfbeb55af..c92bf01fb 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.C
@@ -158,9 +158,9 @@ Foam::autoHexMeshDriver::autoHexMeshDriver
     {
         meshRefinement::debug = debug_;
         autoHexMeshDriver::debug = debug_;
-        autoRefineDriver::debug = debug;
-        autoSnapDriver::debug = debug;
-        autoLayerDriver::debug = debug;
+        autoRefineDriver::debug = debug_;
+        autoSnapDriver::debug = debug_;
+        autoLayerDriver::debug = debug_;
     }
 
     refinementParameters refineParams(dict, 1);

From 9aef7695ce9703881bfde7367132679a475980c2 Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Wed, 10 Jun 2015 11:45:31 +0100
Subject: [PATCH 16/24] Bug fix: handling solo faces in topo change

---
 .../polyTopoChangerChangeMesh.C               | 57 +++++++++++++++++--
 1 file changed, 53 insertions(+), 4 deletions(-)

diff --git a/src/dynamicMesh/dynamicMesh/polyTopoChange/polyTopoChanger/polyTopoChangerChangeMesh.C b/src/dynamicMesh/dynamicMesh/polyTopoChange/polyTopoChanger/polyTopoChangerChangeMesh.C
index 56015ebbe..f82469a91 100644
--- a/src/dynamicMesh/dynamicMesh/polyTopoChange/polyTopoChanger/polyTopoChangerChangeMesh.C
+++ b/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.
-    forAll (points, pointI)
+    // Grab the untouched live points
+    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 = "
-            << nNewFaces - debugFaceCounter;
+        Pout<< " modified = " << nNewFaces - debugFaceCounter;
 
             debugFaceCounter = nNewFaces;
     }
 
+    // Freely-standing added faces
     forAll (af, afI)
     {
         if (af[afI].onlyInZone())

From e1fe806a11711cbc3e9cd4213627158a4c01b32e Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Wed, 10 Jun 2015 11:47:03 +0100
Subject: [PATCH 17/24] Clean-up

---
 .../sphereTransport/sphereTransport/system/controlDict        | 1 -
 .../incompressible/icoFoam/cavity/constant/polyMesh/boundary  | 4 ++--
 .../pUCoupledFoam/cavity/constant/polyMesh/boundary           | 4 ++--
 tutorials/incompressible/simpleFoam/pitzDaily/0/epsilon       | 2 +-
 tutorials/incompressible/simpleFoam/pitzDaily/0/k             | 3 ++-
 5 files changed, 7 insertions(+), 7 deletions(-)

diff --git a/tutorials/finiteArea/surfactantFoam/sphereTransport/sphereTransport/system/controlDict b/tutorials/finiteArea/surfactantFoam/sphereTransport/sphereTransport/system/controlDict
index 8c4c0d612..8dd1c7103 100644
--- a/tutorials/finiteArea/surfactantFoam/sphereTransport/sphereTransport/system/controlDict
+++ b/tutorials/finiteArea/surfactantFoam/sphereTransport/sphereTransport/system/controlDict
@@ -44,4 +44,3 @@ timePrecision       6;
 runTimeModifiable   yes;
 
 // ************************************************************************* //
-
diff --git a/tutorials/incompressible/icoFoam/cavity/constant/polyMesh/boundary b/tutorials/incompressible/icoFoam/cavity/constant/polyMesh/boundary
index d290d2f9f..5fe20704d 100644
--- a/tutorials/incompressible/icoFoam/cavity/constant/polyMesh/boundary
+++ b/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
 {
diff --git a/tutorials/incompressible/pUCoupledFoam/cavity/constant/polyMesh/boundary b/tutorials/incompressible/pUCoupledFoam/cavity/constant/polyMesh/boundary
index d0e19cb03..b8ec883be 100644
--- a/tutorials/incompressible/pUCoupledFoam/cavity/constant/polyMesh/boundary
+++ b/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
 {
diff --git a/tutorials/incompressible/simpleFoam/pitzDaily/0/epsilon b/tutorials/incompressible/simpleFoam/pitzDaily/0/epsilon
index b9d0572a9..580f5ddde 100644
--- a/tutorials/incompressible/simpleFoam/pitzDaily/0/epsilon
+++ b/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
 {
diff --git a/tutorials/incompressible/simpleFoam/pitzDaily/0/k b/tutorials/incompressible/simpleFoam/pitzDaily/0/k
index 717e6e04b..55d5121f0 100644
--- a/tutorials/incompressible/simpleFoam/pitzDaily/0/k
+++ b/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
     }
 }
 
+
 // ************************************************************************* //

From 6eb03ca7ce852c1c02e382c830243fe417502b67 Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Wed, 10 Jun 2015 11:47:30 +0100
Subject: [PATCH 18/24] Clean-up

---
 .../incompressible/simpleFoam/pitzDaily/0/nut | 57 +++++++++++++++++++
 .../pitzDaily/constant/polyMesh/boundary      | 52 +++++++++++++++++
 2 files changed, 109 insertions(+)
 create mode 100644 tutorials/incompressible/simpleFoam/pitzDaily/0/nut
 create mode 100644 tutorials/incompressible/simpleFoam/pitzDaily/constant/polyMesh/boundary

diff --git a/tutorials/incompressible/simpleFoam/pitzDaily/0/nut b/tutorials/incompressible/simpleFoam/pitzDaily/0/nut
new file mode 100644
index 000000000..7574f3676
--- /dev/null
+++ b/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;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/simpleFoam/pitzDaily/constant/polyMesh/boundary b/tutorials/incompressible/simpleFoam/pitzDaily/constant/polyMesh/boundary
new file mode 100644
index 000000000..dc903581e
--- /dev/null
+++ b/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;
+    }
+)
+
+// ************************************************************************* //

From badd84700a8144f6aa529f93cdfcffa7ce4b6f99 Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Thu, 11 Jun 2015 15:55:26 +0100
Subject: [PATCH 19/24] Debug check

---
 .../Lists/IndirectList/IndirectListI.H         | 18 +++++++++++++++++-
 1 file changed, 17 insertions(+), 1 deletion(-)

diff --git a/src/foam/containers/Lists/IndirectList/IndirectListI.H b/src/foam/containers/Lists/IndirectList/IndirectListI.H
index 05908bbd4..ae39d8acb 100644
--- a/src/foam/containers/Lists/IndirectList/IndirectListI.H
+++ b/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  * * * * * * * * * * * * * //

From 0d34546d97c2676c9710c624cabb84a6c8418045 Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Thu, 11 Jun 2015 15:55:51 +0100
Subject: [PATCH 20/24] Debug handling for parallel cases, where ggi patch may
 be empty

---
 .../polyPatches/constraint/ggi/ggiPolyPatch.C | 34 ++++++++++++++-----
 1 file changed, 26 insertions(+), 8 deletions(-)

diff --git a/src/foam/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.C b/src/foam/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.C
index 919aef2bb..e1b009318 100644
--- a/src/foam/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.C
+++ b/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;
         }
     }
 }

From 5c3c481e85eed71657f33ea01c862f64238d7590 Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Thu, 11 Jun 2015 15:56:06 +0100
Subject: [PATCH 21/24] Field bounding after topo change

---
 .../compressible/sonicDyMFoam/sonicDyMFoam.C       | 14 +++++++++-----
 1 file changed, 9 insertions(+), 5 deletions(-)

diff --git a/applications/solvers/compressible/sonicDyMFoam/sonicDyMFoam.C b/applications/solvers/compressible/sonicDyMFoam/sonicDyMFoam.C
index 524e4168e..4fb6e2572 100644
--- a/applications/solvers/compressible/sonicDyMFoam/sonicDyMFoam.C
+++ b/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
-    compressible gas with support for mesh motion and topological changes
+    Transient solver for trans-sonic/supersonic for laminar or turbulent
+    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();
-            p.correctBoundaryConditions();
-            e.correctBoundaryConditions();
+            T = max(T, TMin);
+            p = max(p, pMin);
+            e = max(e, thermo.Cv()*TMin);
+
             thermo.correct();
             rho = thermo.rho();
 

From 6f4acbb9c52d75b8c88134ced84f3e457d4791d9 Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Thu, 11 Jun 2015 19:24:41 +0100
Subject: [PATCH 22/24] Formatting

---
 .../lduMatrix/lduMatrixUpdateMatrixInterfaces.C     | 13 ++++++++-----
 1 file changed, 8 insertions(+), 5 deletions(-)

diff --git a/src/foam/matrices/lduMatrix/lduMatrix/lduMatrixUpdateMatrixInterfaces.C b/src/foam/matrices/lduMatrix/lduMatrix/lduMatrixUpdateMatrixInterfaces.C
index 26c7bb731..90afa8b84 100644
--- a/src/foam/matrices/lduMatrix/lduMatrix/lduMatrixUpdateMatrixInterfaces.C
+++ b/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;
-            interfaceI<interfaces.size();
+            label interfaceI = patchSchedule.size()/2;
+            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;
-            interfaceI<interfaces.size();
+            label interfaceI = patchSchedule.size()/2;
+            interfaceI < interfaces.size();
             interfaceI++
         )
         {

From c6ffe26776cb12be8a1071593955e9d31efcc657 Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Thu, 11 Jun 2015 19:28:56 +0100
Subject: [PATCH 23/24] Updated default comms to nonBlocking

---
 src/foam/db/IOstreams/Pstreams/Pstream.C | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/foam/db/IOstreams/Pstreams/Pstream.C b/src/foam/db/IOstreams/Pstreams/Pstream.C
index f5045df3d..8f21cc1a9 100644
--- a/src/foam/db/IOstreams/Pstreams/Pstream.C
+++ b/src/foam/db/IOstreams/Pstreams/Pstream.C
@@ -269,7 +269,8 @@ const Foam::debug::optimisationSwitch
 Foam::Pstream::defaultCommsType
 (
     "commsType",
-    "scheduled",
+    "nonBlocking",
+//     "scheduled",
 //     "blocking",
     "blocking, nonBlocking, scheduled"
 );

From 44ee06943a31a406d300373cfb6f32d344e69342 Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Thu, 11 Jun 2015 19:29:44 +0100
Subject: [PATCH 24/24] Update for non-blocking comms.  NOTE: Scheduled comms
 are broken for coupled matrix

---
 .../coupledLduMatrix/coupledLduMatrix.C       | 130 ++++++++++++++++--
 .../coupledLduMatrix/coupledLduMatrix.H       |   6 +-
 .../coupledGaussSeidelPrecon.C                |   6 +-
 3 files changed, 124 insertions(+), 18 deletions(-)

diff --git a/src/coupledMatrix/coupledLduMatrix/coupledLduMatrix.C b/src/coupledMatrix/coupledLduMatrix/coupledLduMatrix.C
index 93c39a017..b8e7fc68a 100644
--- a/src/coupledMatrix/coupledLduMatrix/coupledLduMatrix.C
+++ b/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
-            (
-                coupleCoeffs[rowI],
-                interfaces[rowI],
-                x[rowI],
-                result[rowI],
-                cmpt
-            );
+            FatalErrorIn("void coupledLduMatrix::initMatrixInterfaces")
+                << "Unsuported communications type "
+                << Pstream::commsTypeNames[Pstream::defaultCommsType()]
+                << exit(FatalError);
         }
     }
 
@@ -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
         );
     }
 }
diff --git a/src/coupledMatrix/coupledLduMatrix/coupledLduMatrix.H b/src/coupledMatrix/coupledLduMatrix/coupledLduMatrix.H
index 67d5a1d22..8d762faf9 100644
--- a/src/coupledMatrix/coupledLduMatrix/coupledLduMatrix.H
+++ b/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;
 
 };
diff --git a/src/coupledMatrix/coupledLduPrecon/GaussSeidelPrecon/coupledGaussSeidelPrecon.C b/src/coupledMatrix/coupledLduPrecon/GaussSeidelPrecon/coupledGaussSeidelPrecon.C
index b7f46cdb6..361518b89 100644
--- a/src/coupledMatrix/coupledLduPrecon/GaussSeidelPrecon/coupledGaussSeidelPrecon.C
+++ b/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
             );
         }