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

applications/solvers/Lagrangian/reactingParcelFoam/Make/files

@@ -1,3 +1,3 @@
 reactingParcelFoam.C
 
-EXE = $(FOAM_USER_APPBIN)/reactingParcelFoam
+EXE = $(FOAM_APPBIN)/reactingParcelFoam

applications/solvers/combustion/coalChemistryFoam/Make/files

@@ -1,3 +1,3 @@
 coalChemistryFoam.C
 
-EXE = $(FOAM_USER_APPBIN)/coalChemistryFoam
+EXE = $(FOAM_APPBIN)/coalChemistryFoam

applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H

@@ -3,7 +3,11 @@
     (
         fvm::ddt(alpha1)
       + fvm::div(phi, alpha1)
-      - fvm::laplacian(Dab, alpha1)
+      - fvm::laplacian
+        (
+            Dab + alphatab*turbulence->nut(), alpha1,
+            "laplacian(Dab,alpha1)"
+        )
     );
 
     alpha1Eqn.solve();

applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H

@@ -50,6 +50,9 @@
 
     dimensionedScalar Dab(twoPhaseProperties.lookup("Dab"));
 
+    // Read the reciprocal of the turbulent Schmidt number
+    dimensionedScalar alphatab(twoPhaseProperties.lookup("alphatab"));
+
     // Need to store rho for ddt(rho, U)
     volScalarField rho("rho", alpha1*rho1 + (scalar(1) - alpha1)*rho2);
     rho.oldTime();

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/decomposeParDict

@@ -55,7 +55,17 @@ metisCoeffs
 }
 
 scotchCoeffs
-{}
+{
+    //processorWeights
+    //(
+    //    1
+    //    1
+    //    1
+    //    1
+    //);
+    //writeGraph  true;
+    //strategy "b";
+}
 
 manualCoeffs
 {

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/patch/patchIntegrate/patchIntegrate.C

@@ -78,7 +78,7 @@ int main(int argc, char *argv[])
             // Give patch area
             if (isType<cyclicPolyPatch>(mesh.boundaryMesh()[patchi]))
             {
-                Info<< "    Cyclic patch area: " << nl;
+                Info<< "    Cyclic patch vector area: " << nl;
                 label nFaces = mesh.boundaryMesh()[patchi].size();
                 vector sum1 = vector::zero;
                 vector sum2 = vector::zero;
@@ -92,12 +92,18 @@ int main(int argc, char *argv[])
                 Info<< "    - half 1 = " << sum1 << ", " << mag(sum1) << nl
                     << "    - half 2 = " << sum2 << ", " << mag(sum2) << nl
                     << "    - total  = " << (sum1 + sum2) << ", "
-                    << mag(sum1 + sum2) << endl;;
+                    << mag(sum1 + sum2) << endl;
+                Info<< "    Cyclic patch area magnitude = "
+                    << gSum(mesh.magSf().boundaryField()[patchi])/2.0 << endl;
             }
             else
             {
-                Info<< "    Patch area = "
+                Info<< "    Area vector of patch "
+                    << patchName << '[' << patchi << ']' << " = "
                     << gSum(mesh.Sf().boundaryField()[patchi]) << endl;
+                Info<< "    Area magnitude of patch "
+                    << patchName << '[' << patchi << ']' << " = "
+                    << gSum(mesh.magSf().boundaryField()[patchi]) << endl;
             }
 
             // Read field and calc integral
@@ -107,15 +113,26 @@ int main(int argc, char *argv[])
                     << fieldName << endl;
 
                 volScalarField field(fieldHeader, mesh);
-                vector sumField = gSum
-                (
-                    mesh.Sf().boundaryField()[patchi]
-                   *field.boundaryField()[patchi]
-                );
 
-                Info<< "    Integral of " << fieldName << " over patch "
+                Info<< "    Integral of " << fieldName
+                    << " over vector area of patch "
                     << patchName << '[' << patchi << ']' << " = "
-                    << sumField << nl;
+                    << gSum
+                       (
+                           mesh.Sf().boundaryField()[patchi]
+                          *field.boundaryField()[patchi]
+                       )
+                    << nl;
+
+                Info<< "    Integral of " << fieldName
+                    << " over area magnitude of patch "
+                    << patchName << '[' << patchi << ']' << " = "
+                    << gSum
+                       (
+                           mesh.magSf().boundaryField()[patchi]
+                          *field.boundaryField()[patchi]
+                       )
+                    << nl;
             }
             else if
             (

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/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndex.H

@@ -84,6 +84,9 @@ public:
         //- Start of procI+1 data
         inline const labelList& offsets() const;
 
+        //- my local size
+        inline label localSize() const;
+
         //- Global sum of localSizes
         inline label size() const;
 

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndexI.H

@@ -34,6 +34,17 @@ inline const Foam::labelList& Foam::globalIndex::offsets() const
 }
 
 
+inline Foam::label Foam::globalIndex::localSize() const
+{
+    return
+    (
+        Pstream::myProcNo() == 0
+      ? offsets_[Pstream::myProcNo()]
+      : offsets_[Pstream::myProcNo()] - offsets_[Pstream::myProcNo()-1]
+    );
+}
+
+
 inline Foam::label Foam::globalIndex::size() const
 {
     return offsets_[Pstream::nProcs()-1];

src/decompositionMethods/decompositionMethods/scotchDecomp/scotchDecomp.C

@@ -75,6 +75,24 @@ extern "C"
 }
 
 
+// Hack: scotch generates floating point errors so need to switch of error
+//       trapping!
+#if defined(linux) || defined(linuxAMD64) || defined(linuxIA64)
+#    define LINUX
+#endif
+
+#if defined(LINUX) && defined(__GNUC__)
+#    define LINUX_GNUC
+#endif
+
+#ifdef LINUX_GNUC
+#   ifndef __USE_GNU
+#       define __USE_GNU
+#   endif
+#   include <fenv.h>
+#endif
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -113,13 +131,30 @@ Foam::label Foam::scotchDecomp::decompose
 {
     // Strategy
     // ~~~~~~~~
+
     // Default.
     SCOTCH_Strat stradat;
     check(SCOTCH_stratInit(&stradat), "SCOTCH_stratInit");
-    //SCOTCH_stratGraphMap(&stradat, &argv[i][2]);
+
-    //fprintf(stdout, "S\tStrat=");
+    if (decompositionDict_.found("scotchCoeffs"))
-    //SCOTCH_stratSave(&stradat, stdout);
+    {
-    //fprintf(stdout, "\n");
+        const dictionary& scotchCoeffs =
+            decompositionDict_.subDict("scotchCoeffs");
+
+
+        string strategy;
+        if (scotchCoeffs.readIfPresent("strategy", strategy))
+        {
+            if (debug)
+            {
+                Info<< "scotchDecomp : Using strategy " << strategy << endl;
+            }
+            SCOTCH_stratGraphMap(&stradat, strategy.c_str());
+            //fprintf(stdout, "S\tStrat=");
+            //SCOTCH_stratSave(&stradat, stdout);
+            //fprintf(stdout, "\n");
+        }
+    }
 
 
     // Graph
@@ -153,37 +188,40 @@ Foam::label Foam::scotchDecomp::decompose
         const dictionary& scotchCoeffs =
             decompositionDict_.subDict("scotchCoeffs");
 
-        Switch writeGraph(scotchCoeffs.lookup("writeGraph"));
+        if (scotchCoeffs.found("writeGraph"))
-
-        if (writeGraph)
         {
-            OFstream str(mesh_.time().path() / mesh_.name() + ".grf");
+            Switch writeGraph(scotchCoeffs.lookup("writeGraph"));
 
-            Info<< "Dumping Scotch graph file to " << str.name() << endl
+            if (writeGraph)
-                << "Use this in combination with gpart." << endl;
-
-            label version = 0;
-            str << version << nl;
-            // Numer of vertices
-            str << xadj.size()-1 << ' ' << adjncy.size() << nl;
-            // Numbering starts from 0
-            label baseval = 0;
-            // Has weights?
-            label hasEdgeWeights = 0;
-            label hasVertexWeights = 0;
-            label numericflag = 10*hasEdgeWeights+hasVertexWeights;
-            str << baseval << ' ' << numericflag << nl;
-            for (label cellI = 0; cellI < xadj.size()-1; cellI++)
             {
-                label start = xadj[cellI];
+                OFstream str(mesh_.time().path() / mesh_.name() + ".grf");
-                label end = xadj[cellI+1];
-                str << end-start;
 
-                for (label i = start; i < end; i++)
+                Info<< "Dumping Scotch graph file to " << str.name() << endl
+                    << "Use this in combination with gpart." << endl;
+
+                label version = 0;
+                str << version << nl;
+                // Numer of vertices
+                str << xadj.size()-1 << ' ' << adjncy.size() << nl;
+                // Numbering starts from 0
+                label baseval = 0;
+                // Has weights?
+                label hasEdgeWeights = 0;
+                label hasVertexWeights = 0;
+                label numericflag = 10*hasEdgeWeights+hasVertexWeights;
+                str << baseval << ' ' << numericflag << nl;
+                for (label cellI = 0; cellI < xadj.size()-1; cellI++)
                 {
-                    str << ' ' << adjncy[i];
+                    label start = xadj[cellI];
+                    label end = xadj[cellI+1];
+                    str << end-start;
+
+                    for (label i = start; i < end; i++)
+                    {
+                        str << ' ' << adjncy[i];
+                    }
+                    str << nl;
                 }
-                str << nl;
             }
         }
     }
@@ -195,12 +233,36 @@ Foam::label Foam::scotchDecomp::decompose
 
     SCOTCH_Arch archdat;
     check(SCOTCH_archInit(&archdat), "SCOTCH_archInit");
-    check
+
-    (
+    List<label> processorWeights;
-        // SCOTCH_archCmpltw for weighted.
+    if (decompositionDict_.found("scotchCoeffs"))
-        SCOTCH_archCmplt(&archdat, nProcessors_),
+    {
-        "SCOTCH_archCmplt"
+        const dictionary& scotchCoeffs =
-    );
+            decompositionDict_.subDict("scotchCoeffs");
+
+        scotchCoeffs.readIfPresent("processorWeights", processorWeights);
+    }
+    if (processorWeights.size())
+    {
+        if (debug)
+        {
+            Info<< "scotchDecomp : Using procesor weights " << processorWeights
+                << endl;
+        }
+        check
+        (
+            SCOTCH_archCmpltw(&archdat, nProcessors_, processorWeights.begin()),
+            "SCOTCH_archCmpltw"
+        );
+    }
+    else
+    {
+        check
+        (
+            SCOTCH_archCmplt(&archdat, nProcessors_),
+            "SCOTCH_archCmplt"
+        );
+    }
 
 
     //SCOTCH_Mapping mapdat;
@@ -209,6 +271,16 @@ Foam::label Foam::scotchDecomp::decompose
     //SCOTCH_graphMapExit(&grafdat, &mapdat);
 
 
+    // Hack:switch off fpu error trapping
+#   ifdef LINUX_GNUC
+    int oldExcepts = fedisableexcept
+    (
+        FE_DIVBYZERO
+      | FE_INVALID
+      | FE_OVERFLOW
+    );
+#   endif
+
     finalDecomp.setSize(xadj.size()-1);
     finalDecomp = 0;
     check
@@ -223,6 +295,11 @@ Foam::label Foam::scotchDecomp::decompose
         "SCOTCH_graphMap"
     );
 
+#   ifdef LINUX_GNUC
+    feenableexcept(oldExcepts);
+#   endif
+
+
 
     //finalDecomp.setSize(xadj.size()-1);
     //check

src/finiteVolume/fields/fvPatchFields/derived/fixedFluxBuoyantPressure/fixedFluxBuoyantPressureFvPatchScalarField.C

@@ -113,7 +113,17 @@ void fixedFluxBuoyantPressureFvPatchScalarField::updateCoeffs()
     const fvPatchField<scalar>& rho =
         patch().lookupPatchField<volScalarField, scalar>("rho");
 
-    gradient() = -rho.snGrad()*(g.value() & patch().Cf());
+    // If the variable name is "pd" assume it is p - rho*g.h
+    // and set the gradient appropriately.
+    // Otherwise assume the variable is the static pressure.
+    if (dimensionedInternalField().name() == "pd")
+    {
+        gradient() = -rho.snGrad()*(g.value() & patch().Cf());
+    }
+    else
+    {
+        gradient() = rho*(g.value() & patch().nf());
+    }
 
     fixedGradientFvPatchScalarField::updateCoeffs();
 }

src/finiteVolume/fields/fvPatchFields/derived/fixedFluxBuoyantPressure/fixedFluxBuoyantPressureFvPatchScalarField.H

@@ -26,7 +26,10 @@ Class
     Foam::fixedFluxBuoyantPressureFvPatchScalarField
 
 Description
-    Foam::fixedFluxBuoyantPressureFvPatchScalarField
+    Set the pressure gradient boundary condition appropriately for buoyant flow.
+
+    If the variable name is "pd" assume it is p - rho*g.h and set the gradient
+    appropriately.  Otherwise assume the variable is the static pressure.
 
 SourceFiles
     fixedFluxBuoyantPressureFvPatchScalarField.C

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/lagrangian/coalCombustion/Make/files

@@ -19,4 +19,4 @@ submodels/surfaceReactionModel/COxidationKineticDiffusionLimitedRate/COxidationK
 submodels/surfaceReactionModel/COxidationMurphyShaddix/COxidationMurphyShaddix.C
 
 
-LIB = $(FOAM_USER_LIBBIN)/libcoalCombustion
+LIB = $(FOAM_LIBBIN)/libcoalCombustion

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: