From f983e0b42e623d2e72e3de76e492d7947d3f2936 Mon Sep 17 00:00:00 2001
From: henry <Henry Weller h.weller@opencfd.co.uk>
Date: Mon, 11 May 2009 21:19:26 +0100
Subject: [PATCH 1/7] Correct handling of MANPATH so that the system settings
 are not overridden.

---
 etc/apps/cint/bashrc | 2 +-
 etc/bashrc           | 2 +-
 etc/cshrc            | 2 +-
 3 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/etc/apps/cint/bashrc b/etc/apps/cint/bashrc
index 577923968a..3ecfd85517 100644
--- a/etc/apps/cint/bashrc
+++ b/etc/apps/cint/bashrc
@@ -41,7 +41,7 @@ do
         export CINTSYSDIR
         export PATH=$PATH:$CINTSYSDIR/bin
         export MANPATH=$MANPATH:$CINTSYSDIR/doc
-        export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:$CINTSYSDIR/lib
+        export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$CINTSYSDIR/lib
         break
     fi
 done
diff --git a/etc/bashrc b/etc/bashrc
index 26dc8c2dfc..17b84d853d 100644
--- a/etc/bashrc
+++ b/etc/bashrc
@@ -256,7 +256,7 @@ cleanEnv=`$cleanProg "$PATH"` && PATH="$cleanEnv"
 cleanEnv=`$cleanProg "$LD_LIBRARY_PATH"` && LD_LIBRARY_PATH="$cleanEnv"
 
 #- Clean MANPATH
-cleanEnv=`$cleanProg "$MANPATH"` && MANPATH="$cleanEnv"
+cleanEnv=`$cleanProg "$MANPATH"` && MANPATH="$cleanEnv:"
 
 export PATH LD_LIBRARY_PATH MANPATH
 
diff --git a/etc/cshrc b/etc/cshrc
index f4147422f9..0073521a33 100644
--- a/etc/cshrc
+++ b/etc/cshrc
@@ -257,7 +257,7 @@ endif
 setenv LD_LIBRARY_PATH `$cleanProg "$LD_LIBRARY_PATH"`
 
 #- Clean MANPATH
-setenv MANPATH `$cleanProg "$MANPATH"`
+setenv MANPATH `$cleanProg "$MANPATH"`:
 
 #- Clean LD_PRELOAD
 if ( $?LD_PRELOAD ) then

From f2d00b32307a1c086062fe0b5a1c42c294be7e05 Mon Sep 17 00:00:00 2001
From: henry <Henry Weller h.weller@opencfd.co.uk>
Date: Mon, 11 May 2009 21:19:55 +0100
Subject: [PATCH 2/7] Improve reference handling.

---
 .../cfdTools/general/findRefCell/findRefCell.C    | 13 ++++++++++++-
 .../cfdTools/general/findRefCell/findRefCell.H    | 15 +++++++++++----
 src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C   | 12 ++++--------
 src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H   |  2 +-
 4 files changed, 28 insertions(+), 14 deletions(-)

diff --git a/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C b/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C
index d815b809a3..ed450ea488 100644
--- a/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C
+++ b/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C
@@ -34,7 +34,7 @@ void Foam::setRefCell
     const dictionary& dict,
     label& refCelli,
     scalar& refValue,
-    bool forceReference
+    const bool forceReference
 )
 {
     if (field.needReference() || forceReference)
@@ -119,4 +119,15 @@ void Foam::setRefCell
 }
 
 
+Foam::scalar Foam::getRefCellValue
+(
+    const volScalarField& field,
+    const label refCelli
+)
+{
+    scalar refCellValue = (refCelli >= 0 ? field[refCelli] : 0.0);
+    return returnReduce<label>(refCellValue, sumOp<scalar>());
+}
+
+
 // ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.H b/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.H
index 9f09a9efea..4017d3c3eb 100644
--- a/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.H
+++ b/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.H
@@ -49,11 +49,18 @@ namespace Foam
 //  but which is not on a cyclic, symmetry or processor patch.
 void setRefCell
 (
-    const volScalarField&,
-    const dictionary&,
-    label& refCellI,
+    const volScalarField& field,
+    const dictionary& dict,
+    label& refCelli,
     scalar& refValue,
-    bool forceReference = false
+    const bool forceReference = false
+);
+
+//- Return the current value of field in the reference cell
+scalar getRefCellValue
+(
+    const volScalarField& field,
+    const label refCelli
 );
 
 }
diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
index 78fede73f8..5ad5dca8cd 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
@@ -472,22 +472,18 @@ void Foam::fvMatrix<Type>::setValues
 }
 
 
-// Set reference level for solution
 template<class Type>
 void Foam::fvMatrix<Type>::setReference
 (
-    const label cell,
+    const label celli,
     const Type& value,
     const bool forceReference
 )
 {
-    if (psi_.needReference() || forceReference)
+    if (celli >= 0 && (psi_.needReference() || forceReference))
     {
-        if (cell >= 0)
-        {
-            source()[cell] += diag()[cell]*value;
-            diag()[cell] += diag()[cell];
-        }
+        source()[celli] += diag()[celli]*value;
+        diag()[celli] += diag()[celli];
     }
 }
 
diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H
index 8ae502cb57..4fbc2430d7 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H
@@ -326,7 +326,7 @@ public:
             //- Set reference level for solution
             void setReference
             (
-                const label cell,
+                const label celli,
                 const Type& value,
                 const bool forceReference = false
             );

From 0ee817a060d196ffd10fd5f1d91f5133121ff6fa Mon Sep 17 00:00:00 2001
From: henry <Henry Weller h.weller@opencfd.co.uk>
Date: Mon, 11 May 2009 23:19:56 +0100
Subject: [PATCH 3/7] Also added a separate specification of the static
 pressure level but the location for the specification is the same as that for
 the pd level.

---
 .../multiphase/interDyMFoam/Make/options        |  6 ++----
 .../multiphase/interDyMFoam/createFields.H      | 17 +++--------------
 .../multiphase/interDyMFoam/interDyMFoam.C      |  3 +--
 .../solvers/multiphase/interFoam/createFields.H | 17 +++++++++++++++++
 .../solvers/multiphase/interFoam/interFoam.C    | 10 ++++++++++
 .../multiphaseInterFoam/createFields.H          | 17 +++++++++++++++++
 .../multiphaseInterFoam/multiphaseInterFoam.C   | 10 ++++++++++
 .../twoLiquidMixingFoam/createFields.H          | 17 +++++++++++++++++
 .../twoLiquidMixingFoam/twoLiquidMixingFoam.C   | 10 ++++++++++
 .../ras/damBreakWithObstacle/system/fvSolution  |  9 +++------
 .../ras/sloshingTank2D/system/fvSolution        | 15 +++------------
 .../ras/sloshingTank2D3DoF/system/fvSolution    | 15 +++------------
 .../ras/sloshingTank3D/system/fvSolution        | 15 +++------------
 .../ras/sloshingTank3D3DoF/system/fvSolution    | 15 +++------------
 .../ras/sloshingTank3D6DoF/system/fvSolution    | 15 +++------------
 15 files changed, 105 insertions(+), 86 deletions(-)

diff --git a/applications/solvers/multiphase/interDyMFoam/Make/options b/applications/solvers/multiphase/interDyMFoam/Make/options
index 74498326a9..b9aa770576 100644
--- a/applications/solvers/multiphase/interDyMFoam/Make/options
+++ b/applications/solvers/multiphase/interDyMFoam/Make/options
@@ -7,8 +7,7 @@ EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude
 
 EXE_LIBS = \
     -linterfaceProperties \
@@ -18,5 +17,4 @@ EXE_LIBS = \
     -lfiniteVolume \
     -ldynamicMesh \
     -lmeshTools \
-    -ldynamicFvMesh \
-    -lsampling
+    -ldynamicFvMesh
diff --git a/applications/solvers/multiphase/interDyMFoam/createFields.H b/applications/solvers/multiphase/interDyMFoam/createFields.H
index 0f1efe717f..bcceb9d748 100644
--- a/applications/solvers/multiphase/interDyMFoam/createFields.H
+++ b/applications/solvers/multiphase/interDyMFoam/createFields.H
@@ -116,34 +116,23 @@
         pd + rho*(g & mesh.C())
     );
 
-    autoPtr<probes> pRefProbe;
     label pdRefCell = 0;
     scalar pdRefValue = 0.0;
+    setRefCell(pd, mesh.solutionDict().subDict("PISO"), pdRefCell, pdRefValue);
+
     scalar pRefValue = 0.0;
 
     if (pd.needReference())
     {
-        pRefProbe.set
-        (
-            new probes
-            (
-                "pRefProbe",
-                mesh,
-                mesh.solutionDict().subDict("PISO").subDict("pRefProbe")
-            )
-        );
-
         pRefValue = readScalar
         (
             mesh.solutionDict().subDict("PISO").lookup("pRefValue")
         );
 
-        pdRefCell = pRefProbe->cells()[0];
-
         p += dimensionedScalar
         (
             "p",
             p.dimensions(),
-            pRefValue - pRefProbe->sample<scalar>("p")()[0]
+            pRefValue - getRefCellValue(p, pdRefCell)
         );
     }
diff --git a/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C b/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C
index 0386e6625c..046503a8c4 100644
--- a/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C
+++ b/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C
@@ -40,7 +40,6 @@ Description
 #include "interfaceProperties.H"
 #include "twoPhaseMixture.H"
 #include "turbulenceModel.H"
-#include "probes.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -123,7 +122,7 @@ int main(int argc, char *argv[])
             (
                 "p",
                 p.dimensions(),
-                pRefValue - pRefProbe->sample<scalar>("p")()[0]
+                pRefValue - getRefCellValue(p, pdRefCell)
             );
         }
 
diff --git a/applications/solvers/multiphase/interFoam/createFields.H b/applications/solvers/multiphase/interFoam/createFields.H
index 53dd01672a..7063774752 100644
--- a/applications/solvers/multiphase/interFoam/createFields.H
+++ b/applications/solvers/multiphase/interFoam/createFields.H
@@ -106,6 +106,23 @@
     scalar pdRefValue = 0.0;
     setRefCell(pd, mesh.solutionDict().subDict("PISO"), pdRefCell, pdRefValue);
 
+    scalar pRefValue = 0.0;
+
+    if (pd.needReference())
+    {
+        pRefValue = readScalar
+        (
+            mesh.solutionDict().subDict("PISO").lookup("pRefValue")
+        );
+
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pdRefCell)
+        );
+    }
+
 
     // Construct interface from alpha1 distribution
     interfaceProperties interface(alpha1, U, twoPhaseProperties);
diff --git a/applications/solvers/multiphase/interFoam/interFoam.C b/applications/solvers/multiphase/interFoam/interFoam.C
index 7d2048037a..4dd86bd8f7 100644
--- a/applications/solvers/multiphase/interFoam/interFoam.C
+++ b/applications/solvers/multiphase/interFoam/interFoam.C
@@ -91,6 +91,16 @@ int main(int argc, char *argv[])
 
         p = pd + rho*gh;
 
+        if (pd.needReference())
+        {
+            p += dimensionedScalar
+            (
+                "p",
+                p.dimensions(),
+                pRefValue - getRefCellValue(p, pdRefCell)
+            );
+        }
+
         turbulence->correct();
 
         runTime.write();
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/createFields.H b/applications/solvers/multiphase/multiphaseInterFoam/createFields.H
index 9119631d0c..07650fd4dd 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/createFields.H
+++ b/applications/solvers/multiphase/multiphaseInterFoam/createFields.H
@@ -68,6 +68,23 @@
     scalar pdRefValue = 0.0;
     setRefCell(pd, mesh.solutionDict().subDict("PISO"), pdRefCell, pdRefValue);
 
+    scalar pRefValue = 0.0;
+
+    if (pd.needReference())
+    {
+        pRefValue = readScalar
+        (
+            mesh.solutionDict().subDict("PISO").lookup("pRefValue")
+        );
+
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pdRefCell)
+        );
+    }
+
 
     // Construct incompressible turbulence model
     autoPtr<incompressible::turbulenceModel> turbulence
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C
index d0626cbf05..c5fdb5f3c2 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C
+++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C
@@ -83,6 +83,16 @@ int main(int argc, char *argv[])
 
         p = pd + rho*gh;
 
+        if (pd.needReference())
+        {
+            p += dimensionedScalar
+            (
+                "p",
+                p.dimensions(),
+                pRefValue - getRefCellValue(p, pdRefCell)
+            );
+        }
+
         turbulence->correct();
 
         runTime.write();
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H b/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H
index bd7c3e2c9f..a38135dc19 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H
@@ -95,6 +95,23 @@
     scalar pdRefValue = 0.0;
     setRefCell(pd, mesh.solutionDict().subDict("PISO"), pdRefCell, pdRefValue);
 
+    scalar pRefValue = 0.0;
+
+    if (pd.needReference())
+    {
+        pRefValue = readScalar
+        (
+            mesh.solutionDict().subDict("PISO").lookup("pRefValue")
+        );
+
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pdRefCell)
+        );
+    }
+
 
     // Construct incompressible turbulence model
     autoPtr<incompressible::turbulenceModel> turbulence
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C b/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C
index be8eea9061..1caf455708 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C
@@ -75,6 +75,16 @@ int main(int argc, char *argv[])
 
         p = pd + rho*gh;
 
+        if (pd.needReference())
+        {
+            p += dimensionedScalar
+            (
+                "p",
+                p.dimensions(),
+                pRefValue - getRefCellValue(p, pdRefCell)
+            );
+        }
+
         turbulence->correct();
 
         runTime.write();
diff --git a/tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/system/fvSolution b/tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/system/fvSolution
index 36d56bb947..2ca306273d 100644
--- a/tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/system/fvSolution
+++ b/tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/system/fvSolution
@@ -121,13 +121,10 @@ PISO
     nAlphaCorr      1;
     nAlphaSubCycles 3;
     cAlpha          1;
-    pdRefCell       0;
+
+    pdRefPoint      (0.51 0.51 0.51);
     pdRefValue      0;
-    pRefProbe
-    {
-        fields          ( pd );
-        probeLocations  ( ( 0.51 0.51 0.51 ) );
-    }
+    pRefValue       0;
 }
 
 
diff --git a/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D/system/fvSolution b/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D/system/fvSolution
index 2688ec1f98..c9dccc239e 100644
--- a/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D/system/fvSolution
+++ b/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D/system/fvSolution
@@ -98,19 +98,10 @@ PISO
     nAlphaSubCycles 3;
     cAlpha          1.5;
     correctPhi      no;
-    pRefProbe
-    {
-        fields          
-        (
-            p
-        );
-        probeLocations
-        (
-            ( 0 0 0.15 )
-        );
-    }
 
-    pRefValue       100000;
+    pdRefPoint      (0 0 0.15);
+    pdRefValue      0;
+    pRefValue       1e5;
 }
 
 relaxationFactors
diff --git a/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D3DoF/system/fvSolution b/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D3DoF/system/fvSolution
index ce22769b82..c9dccc239e 100644
--- a/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D3DoF/system/fvSolution
+++ b/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D3DoF/system/fvSolution
@@ -98,19 +98,10 @@ PISO
     nAlphaSubCycles 3;
     cAlpha          1.5;
     correctPhi      no;
-    pRefProbe
-    {
-        fields
-        (
-            p
-        );
-        probeLocations
-        (
-            ( 0 0 0.15 )
-        );
-    }
 
-    pRefValue       100000;
+    pdRefPoint      (0 0 0.15);
+    pdRefValue      0;
+    pRefValue       1e5;
 }
 
 relaxationFactors
diff --git a/tutorials/multiphase/interDyMFoam/ras/sloshingTank3D/system/fvSolution b/tutorials/multiphase/interDyMFoam/ras/sloshingTank3D/system/fvSolution
index ce22769b82..c9dccc239e 100644
--- a/tutorials/multiphase/interDyMFoam/ras/sloshingTank3D/system/fvSolution
+++ b/tutorials/multiphase/interDyMFoam/ras/sloshingTank3D/system/fvSolution
@@ -98,19 +98,10 @@ PISO
     nAlphaSubCycles 3;
     cAlpha          1.5;
     correctPhi      no;
-    pRefProbe
-    {
-        fields
-        (
-            p
-        );
-        probeLocations
-        (
-            ( 0 0 0.15 )
-        );
-    }
 
-    pRefValue       100000;
+    pdRefPoint      (0 0 0.15);
+    pdRefValue      0;
+    pRefValue       1e5;
 }
 
 relaxationFactors
diff --git a/tutorials/multiphase/interDyMFoam/ras/sloshingTank3D3DoF/system/fvSolution b/tutorials/multiphase/interDyMFoam/ras/sloshingTank3D3DoF/system/fvSolution
index ce22769b82..c9dccc239e 100644
--- a/tutorials/multiphase/interDyMFoam/ras/sloshingTank3D3DoF/system/fvSolution
+++ b/tutorials/multiphase/interDyMFoam/ras/sloshingTank3D3DoF/system/fvSolution
@@ -98,19 +98,10 @@ PISO
     nAlphaSubCycles 3;
     cAlpha          1.5;
     correctPhi      no;
-    pRefProbe
-    {
-        fields
-        (
-            p
-        );
-        probeLocations
-        (
-            ( 0 0 0.15 )
-        );
-    }
 
-    pRefValue       100000;
+    pdRefPoint      (0 0 0.15);
+    pdRefValue      0;
+    pRefValue       1e5;
 }
 
 relaxationFactors
diff --git a/tutorials/multiphase/interDyMFoam/ras/sloshingTank3D6DoF/system/fvSolution b/tutorials/multiphase/interDyMFoam/ras/sloshingTank3D6DoF/system/fvSolution
index ce22769b82..c9dccc239e 100644
--- a/tutorials/multiphase/interDyMFoam/ras/sloshingTank3D6DoF/system/fvSolution
+++ b/tutorials/multiphase/interDyMFoam/ras/sloshingTank3D6DoF/system/fvSolution
@@ -98,19 +98,10 @@ PISO
     nAlphaSubCycles 3;
     cAlpha          1.5;
     correctPhi      no;
-    pRefProbe
-    {
-        fields
-        (
-            p
-        );
-        probeLocations
-        (
-            ( 0 0 0.15 )
-        );
-    }
 
-    pRefValue       100000;
+    pdRefPoint      (0 0 0.15);
+    pdRefValue      0;
+    pRefValue       1e5;
 }
 
 relaxationFactors

From 58aa677a7445c7519ca9a6ec6ea00ec70f4c7f9d Mon Sep 17 00:00:00 2001
From: andy <a.heather@opencfd.co.uk>
Date: Tue, 12 May 2009 18:31:13 +0100
Subject: [PATCH 4/7] bugfix: removed reference to tmp field

---
 .../alphaSgsJayatillekeWallFunctionFvPatchScalarField.C         | 2 +-
 .../alphaSgsWallFunctionFvPatchScalarField.C                    | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/turbulenceModels/compressible/LES/derivedFvPatchFields/wallFunctions/alphaSgsWallFunctions/alphaSgsJayatillekeWallFunction/alphaSgsJayatillekeWallFunctionFvPatchScalarField.C b/src/turbulenceModels/compressible/LES/derivedFvPatchFields/wallFunctions/alphaSgsWallFunctions/alphaSgsJayatillekeWallFunction/alphaSgsJayatillekeWallFunctionFvPatchScalarField.C
index 10838a6dac..b868d8baa5 100644
--- a/src/turbulenceModels/compressible/LES/derivedFvPatchFields/wallFunctions/alphaSgsWallFunctions/alphaSgsJayatillekeWallFunction/alphaSgsJayatillekeWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/compressible/LES/derivedFvPatchFields/wallFunctions/alphaSgsWallFunctions/alphaSgsJayatillekeWallFunction/alphaSgsJayatillekeWallFunctionFvPatchScalarField.C
@@ -191,7 +191,7 @@ void alphaSgsJayatillekeWallFunctionFvPatchScalarField::evaluate
     const label patchI = patch().index();
 
     const scalarField& muw = lesModel.mu().boundaryField()[patchI];
-    const scalarField& muSgsw = lesModel.muSgs()().boundaryField()[patchI];
+    const scalarField muSgsw = lesModel.muSgs()().boundaryField()[patchI];
 
     const scalarField& alphaw = lesModel.alpha().boundaryField()[patchI];
     scalarField& alphaSgsw = *this;
diff --git a/src/turbulenceModels/compressible/LES/derivedFvPatchFields/wallFunctions/alphaSgsWallFunctions/alphaSgsWallFunction/alphaSgsWallFunctionFvPatchScalarField.C b/src/turbulenceModels/compressible/LES/derivedFvPatchFields/wallFunctions/alphaSgsWallFunctions/alphaSgsWallFunction/alphaSgsWallFunctionFvPatchScalarField.C
index 3c697c5c17..5d59e7f6d7 100644
--- a/src/turbulenceModels/compressible/LES/derivedFvPatchFields/wallFunctions/alphaSgsWallFunctions/alphaSgsWallFunction/alphaSgsWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/compressible/LES/derivedFvPatchFields/wallFunctions/alphaSgsWallFunctions/alphaSgsWallFunction/alphaSgsWallFunctionFvPatchScalarField.C
@@ -137,7 +137,7 @@ void alphaSgsWallFunctionFvPatchScalarField::evaluate
     const scalar Prt = lesModel.Prt().value();
 
     // Get the turbulent viscosity at the wall
-    const scalarField& muSgsw =
+    const scalarField muSgsw =
         lesModel.muSgs()().boundaryField()[patch().index()];
 
     operator==(muSgsw/Prt);

From 0d2fd6b2524a3c4dbd32862fafe8a80db00d410c Mon Sep 17 00:00:00 2001
From: henry <Henry Weller h.weller@opencfd.co.uk>
Date: Wed, 13 May 2009 00:02:36 +0100
Subject: [PATCH 5/7] Added stabilisation to the nut update.

---
 .../nutRoughWallFunctionFvPatchScalarField.C  | 56 +++++++++++--------
 1 file changed, 32 insertions(+), 24 deletions(-)

diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutRoughWallFunction/nutRoughWallFunctionFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutRoughWallFunction/nutRoughWallFunctionFvPatchScalarField.C
index 138e783333..42e537295e 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutRoughWallFunction/nutRoughWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutRoughWallFunction/nutRoughWallFunctionFvPatchScalarField.C
@@ -49,21 +49,20 @@ scalar nutRoughWallFunctionFvPatchScalarField::fnRough
     const scalar kappa
 ) const
 {
-    // Set deltaB based on non-dimensional roughness height
-    scalar deltaB = 0.0;
+    // Return fn based on non-dimensional roughness height
+
     if (KsPlus < 90.0)
     {
-        deltaB =
-            1.0/kappa
-            *log((KsPlus - 2.25)/87.75 + Cs*KsPlus)
-            *sin(0.4258*(log(KsPlus) - 0.811));
+        return pow
+        (
+            (KsPlus - 2.25)/87.75 + Cs*KsPlus,
+            sin(0.4258*(log(KsPlus) - 0.811))
+        );
     }
     else
     {
-        deltaB = 1.0/kappa*log(1.0 + Cs*KsPlus);
+        return (1.0 + Cs*KsPlus);
     }
-
-    return exp(min(deltaB*kappa, 50.0));
 }
 
 
@@ -183,7 +182,7 @@ void nutRoughWallFunctionFvPatchScalarField::updateCoeffs()
     const scalar Cmu25 = pow(Cmu, 0.25);
     const scalar kappa = ras.kappa().value();
     const scalar E = ras.E().value();
-    scalar yPlusLam = ras.yPlusLam();
+    const scalar yPlusLam = ras.yPlusLam();
 
     const scalarField& y = ras.y()[patch().index()];
 
@@ -199,17 +198,35 @@ void nutRoughWallFunctionFvPatchScalarField::updateCoeffs()
         label faceCellI = patch().faceCells()[faceI];
 
         scalar uStar = Cmu25*sqrt(k[faceCellI]);
-
         scalar yPlus = uStar*y[faceI]/nuw[faceI];
-
         scalar KsPlus = uStar*Ks_[faceI]/nuw[faceI];
 
         scalar Edash = E;
-        scalar yPlusLamNew = yPlusLam;
+
         if (KsPlus > 2.25)
         {
             Edash = E/fnRough(KsPlus, Cs_[faceI], kappa);
-            yPlusLam = ras.yPlusLam(kappa, Edash);
+        }
+
+        if (yPlus > yPlusLam)
+        {
+            scalar limitingNutw = max(nutw[faceI], nuw[faceI]);
+
+            // To avoid oscillations limit the change in the wall viscosity
+            // which is particularly important if it temporarily becomes zero
+            nutw[faceI] =
+                max
+                (
+                    min
+                    (
+                        nuw[faceI]*(yPlus*kappa/log(Edash*yPlus) - 1),
+                        2*limitingNutw
+                    ), 0.5*limitingNutw
+                );
+        }
+        else
+        {
+            nutw[faceI] = 0.0;
         }
 
         if (debug)
@@ -217,18 +234,9 @@ void nutRoughWallFunctionFvPatchScalarField::updateCoeffs()
             Info<< "yPlus = " << yPlus
                 << ", KsPlus = " << KsPlus
                 << ", Edash = " << Edash
-                << ", yPlusLam = " << yPlusLam
+                << ", nutw = " << nutw[faceI]
                 << endl;
         }
-
-        if (yPlus > yPlusLamNew)
-        {
-            nutw[faceI] = nuw[faceI]*(yPlus*kappa/log(Edash*yPlus) - 1);
-        }
-        else
-        {
-            nutw[faceI] = 0.0;
-        }
     }
 }
 

From a624d515616603f0345d3e9de4f33b2f87416fa1 Mon Sep 17 00:00:00 2001
From: henry <Henry Weller h.weller@opencfd.co.uk>
Date: Fri, 15 May 2009 14:24:42 +0100
Subject: [PATCH 6/7] Added further stabilisation for the case where ksStar is
 very large.

---
 .../nutRoughWallFunctionFvPatchScalarField.C                   | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutRoughWallFunction/nutRoughWallFunctionFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutRoughWallFunction/nutRoughWallFunctionFvPatchScalarField.C
index 42e537295e..7b6e14493e 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutRoughWallFunction/nutRoughWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutRoughWallFunction/nutRoughWallFunctionFvPatchScalarField.C
@@ -219,7 +219,8 @@ void nutRoughWallFunctionFvPatchScalarField::updateCoeffs()
                 (
                     min
                     (
-                        nuw[faceI]*(yPlus*kappa/log(Edash*yPlus) - 1),
+                        nuw[faceI]
+                       *(yPlus*kappa/log(max(Edash*yPlus, 1+1e-4)) - 1),
                         2*limitingNutw
                     ), 0.5*limitingNutw
                 );

From 55e766bd6760d1b7aa08364f2f9f4b392718eace Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Tue, 19 May 2009 15:48:54 +0100
Subject: [PATCH 7/7] Fix: edge lookup from edgeLabels_ by index rather than
 just from edges

---
 src/meshTools/indexedOctree/treeDataEdge.C | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/meshTools/indexedOctree/treeDataEdge.C b/src/meshTools/indexedOctree/treeDataEdge.C
index 12f347c164..3f7ae86ce8 100644
--- a/src/meshTools/indexedOctree/treeDataEdge.C
+++ b/src/meshTools/indexedOctree/treeDataEdge.C
@@ -134,7 +134,7 @@ void Foam::treeDataEdge::findNearest
     {
         label index = indices[i];
 
-        const edge& e = edges_[index];
+        const edge& e = edges_[edgeLabels_[index]];
 
         pointHit nearHit = e.line(points_).nearestDist(sample);
 
@@ -170,7 +170,7 @@ void Foam::treeDataEdge::findNearest
     {
         label index = indices[i];
 
-        const edge& e = edges_[index];
+        const edge& e = edges_[edgeLabels_[index]];
 
         // Note: could do bb test ? Worthwhile?