Merge branch 'master' of ssh://noisy/home/noisy3/OpenFOAM/OpenFOAM-dev

applications/utilities/parallelProcessing/decomposePar/decomposePar.C

@@ -566,7 +566,7 @@ int main(int argc, char *argv[])
             IOobject
             (
                 "cellProcAddressing",
-                "constant",
+                procMesh.facesInstance(),
                 procMesh.meshSubDir,
                 procMesh,
                 IOobject::MUST_READ,
@@ -579,7 +579,7 @@ int main(int argc, char *argv[])
             IOobject
             (
                 "boundaryProcAddressing",
-                "constant",
+                procMesh.facesInstance(),
                 procMesh.meshSubDir,
                 procMesh,
                 IOobject::MUST_READ,
@@ -603,7 +603,7 @@ int main(int argc, char *argv[])
                 IOobject
                 (
                     "faceProcAddressing",
-                    "constant",
+                    procMesh.facesInstance(),
                     procMesh.meshSubDir,
                     procMesh,
                     IOobject::MUST_READ,
@@ -645,7 +645,7 @@ int main(int argc, char *argv[])
                 IOobject
                 (
                     "pointProcAddressing",
-                    "constant",
+                    procMesh.facesInstance(),
                     procMesh.meshSubDir,
                     procMesh,
                     IOobject::MUST_READ,

applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C

@@ -269,7 +269,7 @@ bool domainDecomposition::writeDecomposition()
             IOobject
             (
                 this->polyMesh::name(),  // region name of undecomposed mesh
-                "constant",
+                pointsInstance(),
                 processorDb
             ),
             xferMove(procPoints),
@@ -620,7 +620,7 @@ bool domainDecomposition::writeDecomposition()
             IOobject
             (
                 "pointProcAddressing",
-                "constant",
+                procMesh.facesInstance(),
                 procMesh.meshSubDir,
                 procMesh,
                 IOobject::NO_READ,
@@ -635,7 +635,7 @@ bool domainDecomposition::writeDecomposition()
             IOobject
             (
                 "faceProcAddressing",
-                "constant",
+                procMesh.facesInstance(),
                 procMesh.meshSubDir,
                 procMesh,
                 IOobject::NO_READ,
@@ -650,7 +650,7 @@ bool domainDecomposition::writeDecomposition()
             IOobject
             (
                 "cellProcAddressing",
-                "constant",
+                procMesh.facesInstance(),
                 procMesh.meshSubDir,
                 procMesh,
                 IOobject::NO_READ,
@@ -665,7 +665,7 @@ bool domainDecomposition::writeDecomposition()
             IOobject
             (
                 "boundaryProcAddressing",
-                "constant",
+                procMesh.facesInstance(),
                 procMesh.meshSubDir,
                 procMesh,
                 IOobject::NO_READ,

applications/utilities/postProcessing/wall/yPlusLES/Make/options

@@ -1,4 +1,5 @@
 EXE_INC = \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/turbulenceModels \
     -I$(LIB_SRC)/turbulenceModels/incompressible/LES/LESModel \
     -I$(LIB_SRC)/turbulenceModels/LES/LESdeltas/lnInclude \

applications/utilities/postProcessing/wall/yPlusLES/yPlusLES.C

@@ -34,6 +34,7 @@ Description
 #include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H"
 #include "LESModel.H"
 #include "nearWallDist.H"
+#include "wallDist.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -49,7 +50,18 @@ int main(int argc, char *argv[])
     {
         runTime.setTime(timeDirs[timeI], timeI);
         Info<< "Time = " << runTime.timeName() << endl;
-        mesh.readUpdate();
+        fvMesh::readUpdateState state = mesh.readUpdate();
+
+        // Wall distance
+        if (timeI == 0 || state != fvMesh::UNCHANGED)
+        {
+            Info<< "Calculating wall distance\n" << endl;
+            wallDist y(mesh, true);
+            Info<< "Writing wall distance to field "
+                << y.name() << nl << endl;
+            y.write();
+        }
+
 
         volScalarField yPlus
         (
@@ -116,6 +128,9 @@ int main(int argc, char *argv[])
             }
         }
 
+        Info<< "Writing yPlus to field "
+            << yPlus.name() << nl << endl;
+
         yPlus.write();
     }
 

applications/utilities/postProcessing/wall/yPlusRAS/Make/options

@@ -1,4 +1,5 @@
 EXE_INC = \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/turbulenceModels \
     -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \

applications/utilities/postProcessing/wall/yPlusRAS/yPlusRAS.C

@@ -34,6 +34,7 @@ Description
 #include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H"
 #include "RASModel.H"
 #include "wallFvPatch.H"
+#include "wallDist.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -49,7 +50,17 @@ int main(int argc, char *argv[])
     {
         runTime.setTime(timeDirs[timeI], timeI);
         Info<< "Time = " << runTime.timeName() << endl;
-        mesh.readUpdate();
+        fvMesh::readUpdateState state = mesh.readUpdate();
+
+        // Wall distance
+        if (timeI == 0 || state != fvMesh::UNCHANGED)
+        {
+            Info<< "Calculating wall distance\n" << endl;
+            wallDist y(mesh, true);
+            Info<< "Writing wall distance to field "
+                << y.name() << nl << endl;
+            y.write();
+        }
 
         volScalarField yPlus
         (
@@ -106,6 +117,9 @@ int main(int argc, char *argv[])
             }
         }
 
+        Info<< "Writing yPlus to field "
+            << yPlus.name() << nl << endl;
+
         yPlus.write();
     }
 

src/decompositionMethods/decompositionMethods/scotchDecomp/scotchDecomp.C

@@ -153,6 +153,8 @@ Foam::label Foam::scotchDecomp::decompose
         const dictionary& scotchCoeffs =
             decompositionDict_.subDict("scotchCoeffs");
 
+        if (scotchCoeffs.found("writeGraph"))
+        {
             Switch writeGraph(scotchCoeffs.lookup("writeGraph"));
 
             if (writeGraph)
@@ -187,6 +189,7 @@ Foam::label Foam::scotchDecomp::decompose
                 }
             }
         }
+    }
 
 
     // Architecture

src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C

@@ -480,7 +480,7 @@ void Foam::fvMatrix<Type>::setReference
     const bool forceReference
 )
 {
-    if (celli >= 0 && (psi_.needReference() || forceReference))
+    if ((forceReference || psi_.needReference()) && celli >= 0)
     {
         source()[celli] += diag()[celli]*value;
         diag()[celli] += diag()[celli];

src/finiteVolume/fvMesh/extendedStencil/extendedStencil.C

@@ -31,168 +31,6 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-// Calculates per face a list of global cell/face indices.
-void Foam::extendedStencil::calcFaceStencil
-(
-    const labelListList& globalCellCells,
-    labelListList& faceStencil
-)
-{
-    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-    const label nBnd = mesh_.nFaces()-mesh_.nInternalFaces();
-    const labelList& own = mesh_.faceOwner();
-    const labelList& nei = mesh_.faceNeighbour();
-
-
-    // Determine neighbouring global cell Cells
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    labelListList neiGlobalCellCells(nBnd);
-    forAll(patches, patchI)
-    {
-        const polyPatch& pp = patches[patchI];
-
-        if (pp.coupled())
-        {
-            label faceI = pp.start();
-
-            forAll(pp, i)
-            {
-                neiGlobalCellCells[faceI-mesh_.nInternalFaces()] =
-                    globalCellCells[own[faceI]];
-                faceI++;
-            }
-        }
-    }
-    syncTools::swapBoundaryFaceList(mesh_, neiGlobalCellCells, false);
-
-
-
-    // Construct stencil in global numbering
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    faceStencil.setSize(mesh_.nFaces());
-
-    labelHashSet faceStencilSet;
-
-    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
-    {
-        faceStencilSet.clear();
-
-        const labelList& ownCCells = globalCellCells[own[faceI]];
-        label globalOwn = ownCCells[0];
-        // Insert cellCells
-        forAll(ownCCells, i)
-        {
-            faceStencilSet.insert(ownCCells[i]);
-        }
-
-        const labelList& neiCCells = globalCellCells[nei[faceI]];
-        label globalNei = neiCCells[0];
-        // Insert cellCells
-        forAll(neiCCells, i)
-        {
-            faceStencilSet.insert(neiCCells[i]);
-        }
-
-        // Guarantee owner first, neighbour second.
-        faceStencil[faceI].setSize(faceStencilSet.size());
-        label n = 0;
-        faceStencil[faceI][n++] = globalOwn;
-        faceStencil[faceI][n++] = globalNei;
-        forAllConstIter(labelHashSet, faceStencilSet, iter)
-        {
-            if (iter.key() != globalOwn && iter.key() != globalNei)
-            {
-                faceStencil[faceI][n++] = iter.key();
-            }
-        }
-        //Pout<< "internalface:" << faceI << " toc:" << faceStencilSet.toc()
-        //    << " faceStencil:" << faceStencil[faceI] << endl;
-    }
-    forAll(patches, patchI)
-    {
-        const polyPatch& pp = patches[patchI];
-        label faceI = pp.start();
-
-        if (pp.coupled())
-        {
-            forAll(pp, i)
-            {
-                faceStencilSet.clear();
-
-                const labelList& ownCCells = globalCellCells[own[faceI]];
-                label globalOwn = ownCCells[0];
-                forAll(ownCCells, i)
-                {
-                    faceStencilSet.insert(ownCCells[i]);
-                }
-
-                // And the neighbours of the coupled cell
-                const labelList& neiCCells =
-                    neiGlobalCellCells[faceI-mesh_.nInternalFaces()];
-                label globalNei = neiCCells[0];
-                forAll(neiCCells, i)
-                {
-                    faceStencilSet.insert(neiCCells[i]);
-                }
-
-                // Guarantee owner first, neighbour second.
-                faceStencil[faceI].setSize(faceStencilSet.size());
-                label n = 0;
-                faceStencil[faceI][n++] = globalOwn;
-                faceStencil[faceI][n++] = globalNei;
-                forAllConstIter(labelHashSet, faceStencilSet, iter)
-                {
-                    if (iter.key() != globalOwn && iter.key() != globalNei)
-                    {
-                        faceStencil[faceI][n++] = iter.key();
-                    }
-                }
-
-                //Pout<< "coupledface:" << faceI
-                //    << " toc:" << faceStencilSet.toc()
-                //    << " faceStencil:" << faceStencil[faceI] << endl;
-
-                faceI++;
-            }
-        }
-        else if (!isA<emptyPolyPatch>(pp))
-        {
-            forAll(pp, i)
-            {
-                faceStencilSet.clear();
-
-                const labelList& ownCCells = globalCellCells[own[faceI]];
-                label globalOwn = ownCCells[0];
-                forAll(ownCCells, i)
-                {
-                    faceStencilSet.insert(ownCCells[i]);
-                }
-
-                // Guarantee owner first, neighbour second.
-                faceStencil[faceI].setSize(faceStencilSet.size());
-                label n = 0;
-                faceStencil[faceI][n++] = globalOwn;
-                forAllConstIter(labelHashSet, faceStencilSet, iter)
-                {
-                    if (iter.key() != globalOwn)
-                    {
-                        faceStencil[faceI][n++] = iter.key();
-                    }
-                }
-
-                //Pout<< "boundaryface:" << faceI
-                //    << " toc:" << faceStencilSet.toc()
-                //    << " faceStencil:" << faceStencil[faceI] << endl;
-
-                faceI++;
-            }
-        }
-    }
-}
-
-
 Foam::autoPtr<Foam::mapDistribute> Foam::extendedStencil::calcDistributeMap
 (
     const globalIndex& globalNumbering,

src/finiteVolume/fvMesh/extendedStencil/extendedStencil.H

@@ -75,13 +75,6 @@ protected:
 
     // Protected Member Functions
 
-        //- Collect cell neighbours into extended stencil
-        void calcFaceStencil
-        (
-            const labelListList& globalCellCells,
-            labelListList& faceStencil
-        );
-
         //- Calculate distribute map
         autoPtr<mapDistribute> calcDistributeMap
         (

src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C

@@ -191,6 +191,13 @@ void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::calcFit
         );
     }
 
+    // Additional weighting for constant and linear terms
+    for(label i = 0; i < B.n(); i++)
+    {
+        B[i][0] *= wts[0];
+        B[i][1] *= wts[0];
+    }
+
     // Set the fit
     label stencilSize = C.size();
     coeffsi.setSize(stencilSize);
@@ -205,7 +212,7 @@ void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::calcFit
 
         for(label i=0; i<stencilSize; i++)
         {
-            coeffsi[i] = wts[i]*svd.VSinvUt()[0][i];
+            coeffsi[i] = wts[0]*wts[i]*svd.VSinvUt()[0][i];
             if (mag(coeffsi[i]) > maxCoeff)
             {
                 maxCoeff = mag(coeffsi[i]);
@@ -229,7 +236,8 @@ void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::calcFit
             //     << "    sing vals " << svd.S() << endl;
         // }
 
-        if (!goodFit) // (not good fit so increase weight in the centre)
+        if (!goodFit) // (not good fit so increase weight in the centre and weight
+                      //  for constant and linear terms)
         {
             // if (iIt == 7)
             // {
@@ -237,7 +245,10 @@ void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::calcFit
             //     (
             //         "FitData<Polynomial>::calcFit"
             //         "(const List<point>& C, const label facei"
-            //     )   << "Cannot fit face " << facei
+            //     )   << "Cannot fit face " << facei << " iteration " << iIt
+            //         << " with sum of weights " << sum(coeffsi) << nl
+            //         << "    Weights " << coeffsi << nl
+            //         << "    Linear weights " << wLin << " " << 1 - wLin << nl
             //         << "    sing vals " << svd.S() << endl;
             // }
 
@@ -249,6 +260,12 @@ void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::calcFit
                 B[0][j] *= 10;
                 B[1][j] *= 10;
             }
+
+            for(label i = 0; i < B.n(); i++)
+            {
+                B[i][0] *= 10;
+                B[i][1] *= 10;
+            }
         }
     }
 

src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutRoughWallFunction/mutRoughWallFunctionFvPatchScalarField.C

@@ -45,25 +45,23 @@ namespace RASModels
 scalar mutRoughWallFunctionFvPatchScalarField::fnRough
 (
     const scalar KsPlus,
-    const scalar Cs,
+    const scalar Cs
-    const scalar kappa
 ) const
 {
-    // Set deltaB based on non-dimensional roughness height
+    // Return fn based on non-dimensional roughness height
-    scalar deltaB = 0.0;
+
     if (KsPlus < 90.0)
     {
-        deltaB =
+        return pow
-            1.0/kappa
+        (
-            *log((KsPlus - 2.25)/87.75 + Cs*KsPlus)
+            (KsPlus - 2.25)/87.75 + Cs*KsPlus,
-            *sin(0.4258*(log(KsPlus) - 0.811));
+            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));
 }
 
 
@@ -216,8 +214,8 @@ void mutRoughWallFunctionFvPatchScalarField::updateCoeffs()
         scalar yPlusLamNew = yPlusLam;
         if (KsPlus > 2.25)
         {
-            Edash = E/fnRough(KsPlus, Cs_[faceI], kappa);
+            Edash = E/fnRough(KsPlus, Cs_[faceI]);
-            yPlusLam = rasModel.yPlusLam(kappa, Edash);
+            yPlusLamNew = rasModel.yPlusLam(kappa, Edash);
         }
 
         if (debug)
@@ -231,7 +229,9 @@ void mutRoughWallFunctionFvPatchScalarField::updateCoeffs()
 
         if (yPlus > yPlusLamNew)
         {
-            mutw[faceI] = muw[faceI]*(yPlus*kappa/log(Edash*yPlus) - 1);
+            mutw[faceI] =
+                muw[faceI]
+               *(yPlus*kappa/log(max(Edash*yPlus, 1+1e-4)) - 1);
         }
         else
         {

src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutRoughWallFunction/mutRoughWallFunctionFvPatchScalarField.H

@@ -83,12 +83,7 @@ class mutRoughWallFunctionFvPatchScalarField
     // Private member functions
 
         //- Compute the roughness function
-        scalar fnRough
+        scalar fnRough(const scalar KsPlus, const scalar Cs) const;
-        (
-            const scalar KsPlus,
-            const scalar Cs,
-            const scalar kappa
-        ) const;
 
 
 public:

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutRoughWallFunction/nutRoughWallFunctionFvPatchScalarField.C

@@ -45,8 +45,7 @@ namespace RASModels
 scalar nutRoughWallFunctionFvPatchScalarField::fnRough
 (
     const scalar KsPlus,
-    const scalar Cs,
+    const scalar Cs
-    const scalar kappa
 ) const
 {
     // Return fn based on non-dimensional roughness height
@@ -205,7 +204,7 @@ void nutRoughWallFunctionFvPatchScalarField::updateCoeffs()
 
         if (KsPlus > 2.25)
         {
-            Edash = E/fnRough(KsPlus, Cs_[faceI], kappa);
+            Edash = E/fnRough(KsPlus, Cs_[faceI]);
         }
 
         if (yPlus > yPlusLam)

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutRoughWallFunction/nutRoughWallFunctionFvPatchScalarField.H

@@ -80,12 +80,7 @@ class nutRoughWallFunctionFvPatchScalarField
     // Private member functions
 
         //- Compute the roughness function
-        scalar fnRough
+        scalar fnRough(const scalar KsPlus, const scalar Cs) const;
-        (
-            const scalar KsPlus,
-            const scalar Cs,
-            const scalar kappa
-        ) const;
 
 
 public: