Merge commit 'OpenCFD/master' into olesenm

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

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)
         );
     }

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)
             );
         }
 

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);

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();

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

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();

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

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();

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

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
 

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

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>());
+}
+
+
 // ************************************************************************* //

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
 );
 
 }

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];
     }
 }
 

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
             );

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?
 

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;

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);

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));
 }
 
 
@@ -177,15 +176,15 @@ void nutRoughWallFunctionFvPatchScalarField::rmap
 
 void nutRoughWallFunctionFvPatchScalarField::updateCoeffs()
 {
-    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
+    const RASModel& ras = db().lookupObject<RASModel>("RASProperties");
 
-    const scalar Cmu = rasModel.Cmu().value();
+    const scalar Cmu = ras.Cmu().value();
     const scalar Cmu25 = pow(Cmu, 0.25);
-    const scalar kappa = rasModel.kappa().value();
-    const scalar E = rasModel.E().value();
-    scalar yPlusLam = rasModel.yPlusLam();
+    const scalar kappa = ras.kappa().value();
+    const scalar E = ras.E().value();
+    const scalar yPlusLam = ras.yPlusLam();
 
-    const scalarField& y = rasModel.y()[patch().index()];
+    const scalarField& y = ras.y()[patch().index()];
 
     const scalarField& k = db().lookupObject<volScalarField>(kName_);
 
@@ -199,17 +198,36 @@ 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 = rasModel.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(max(Edash*yPlus, 1+1e-4)) - 1),
+                        2*limitingNutw
+                    ), 0.5*limitingNutw
+                );
+        }
+        else
+        {
+            nutw[faceI] = 0.0;
         }
 
         if (debug)
@@ -217,18 +235,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;
-        }
     }
 }