Merge branch 'issues-2021-1' into 'develop'

BUG: 2021-1: Various bug fixes and developments

See merge request Development/openfoam!489

applications/utilities/preProcessing/applyBoundaryLayer/applyBoundaryLayer.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2019 OpenCFD Ltd.
+    Copyright (C) 2015-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -333,10 +333,10 @@ int main(int argc, char *argv[])
     // u/U0 = (y/ybl)^(1/7)
     // assumes U0 is the same as the current cell velocity
     Info<< "Setting boundary layer velocity" << nl << endl;
-    scalar yblv = ybl.value();
+    const scalar yblv = ybl.value();
     forAll(U, celli)
     {
-        if (y[celli] <= yblv)
+        if ((y[celli] > 0) && (y[celli] <= yblv))
         {
             mask[celli] = 1;
             U[celli] *= ::pow(y[celli]/yblv, (1.0/7.0));

bin/foamMonitor

@@ -7,7 +7,7 @@
 #    \\/     M anipulation  |
 #-------------------------------------------------------------------------------
 #    Copyright (C) 2015 OpenFOAM Foundation
-#    Copyright (C) 2019 OpenCFD Ltd.
+#    Copyright (C) 2019-2021 OpenCFD Ltd.
 #------------------------------------------------------------------------------
 # License
 #     This file is part of OpenFOAM.
@@ -34,35 +34,49 @@
 #     - requires gnuplot, gnuplot_x11
 #
 #------------------------------------------------------------------------------
-usage() {
+printHelp() {
-    exec 1>&2
-    while [ "$#" -ge 1 ]; do echo "$1"; shift; done
     cat<<USAGE
 
-Usage: ${0##*/} [OPTION] <file>
+Usage: ${0##*/} [OPTIONS] <file>
-options:
+Options:
-  -h | -help            prints the usage
+  -g | -grid            Draw grid lines
-  -i | -idle <time>     stops if <file> unchanging for <time> sec (default = 60)
+  -i | -idle <time>     Stop if <file> unchanging for <time> sec (default = 60)
-  -l | -logscale        plots data (y-axis) on log scale, e.g. for residuals
+  -l | -logscale        Plot y-axis data on log scale
-  -r | -refresh <time>  refreshes display every <time> sec (default = 10)
+  -r | -refresh <time>  Refresh display every <time> sec (default = 10)
-  -y | -yrange <range>  sets data (y-axis) <range>, format "[0:1]"
+  -x | -xrange <range>  Set <range> of x-axis data, format "[0:1]"
-  -g | -grid            draws grid lines on the plot
+  -y | -yrange <range>  Set <range> of y-axis data, format "[0:1]"
+  -h | -help            Display short help and exit
 
 Monitor data with Gnuplot from time-value(s) graphs written by OpenFOAM
-e.g. by functionObjects
+e.g. by functionObjects. For example,
-- requires gnuplot, gnuplot_x11
 
-Example:
     foamMonitor -l postProcessing/residuals/0/residuals.dat
 
 USAGE
+    exit 0  # A clean exit
+}
+
+
+# Report error and exit
+die()
+{
+    exec 1>&2
+    echo
+    echo "Error encountered:"
+    while [ "$#" -ge 1 ]; do echo "    $1"; shift; done
+    echo
+    echo "See '${0##*/} -help' for usage"
+    echo
     exit 1
 }
 
+
+# Set Gnuplot header
 plotFileHeader() {
     cat<<EOF
 set term x11 1 font "helvetica,17" linewidth 1.5 persist noraise
 $LOGSCALE
+$XRANGE
 $YRANGE
 $GRID
 set title "Data Monitoring"
@@ -71,34 +85,49 @@ plot \\
 EOF
 }
 
+
+# Set Gnuplot footer
 plotFileFooter() {
     cat<<EOF
-
 pause $REFRESH
 reread
 EOF
 }
 
-howMany() ( set -f; set -- $1; echo $# )
 
+# Count number of tokens in a variable
+howMany() {
+    ( set -f; set -- $1; echo $# )
+}
+
+
+#-------------------------------------------------------------------------------
 IDLE=60
 REFRESH=10
 LOGSCALE=""
+XRANGE=""
 YRANGE=""
 GRID=""
 GNUPLOT=$(which gnuplot)
-! [ "x$GNUPLOT" = "x" ] || usage "Gnuplot not installed"
+[ ! "$GNUPLOT" = "" ] || die "foamMonitor requires Gnuplot installed"
 
-# parse options
+#-------------------------------------------------------------------------------
+
+# Parse options
 while [ "$#" -gt 0 ]
 do
     case "$1" in
     -h | -help*)
-        usage
+        printHelp
         ;;
     -i | -idle)
-        [ "$#" -ge 2 ] || usage "'$1' option requires an argument"
+        [ "$#" -ge 2 ] || die "'$1' option requires an argument"
-        [ ! -z "${2##*[!0-9]*}" ] && IDLE=$2 || usage "Argument of '$1' is not an integer: '$2'"
+        if [ -n "${2##*[!0-9]*}" ]
+        then
+            IDLE=$2
+        else
+            die "Argument of '$1' is not an integer: '$2'"
+        fi
         shift 2
         ;;
     -l | -logscale)
@@ -106,12 +135,22 @@ do
         shift 1
         ;;
     -r | -refresh)
-        [ "$#" -ge 2 ] || usage "'$1' option requires an argument"
+        [ "$#" -ge 2 ] || die "'$1' option requires an argument"
-        [ ! -z "${2##*[!0-9]*}" ] && REFRESH=$2 || usage "Argument of '$1' is not an integer: '$2'"
+        if [ -n "${2##*[!0-9]*}" ]
+        then
+            REFRESH=$2
+        else
+            die "Argument of '$1' is not an integer: '$2'"
+        fi
+        shift 2
+        ;;
+    -x | -xrange)
+        [ "$#" -ge 2 ] || die "'$1' option requires an argument"
+        XRANGE="set xrange $2"
         shift 2
         ;;
     -y | -yrange)
-        [ "$#" -ge 2 ] || usage "'$1' option requires an argument"
+        [ "$#" -ge 2 ] || die "'$1' option requires an argument"
         YRANGE="set yrange $2"
         shift 2
         ;;
@@ -120,7 +159,7 @@ do
         shift 1
         ;;
     -*)
-        usage "unknown option: '$*'"
+        die "unknown option: '$*'"
         ;;
     *)
         break
@@ -128,28 +167,28 @@ do
     esac
 done
 
-[ $# -eq 1 ] || usage "Incorrect arguments specified"
+[ "$#" -eq 1 ] || die "Incorrect arguments specified"
-[ -f $1 ]    || usage "File $1 does not exit"
+[ -f "$1" ]    || die "File $1 does not exit"
-FILE=$1
+FILE="$1"
 
 # Get KEYS from header
-KEYS=$(grep -E '^#' $FILE | tail -1)
+KEYS=$(grep -E '^#' "$FILE" | tail -1)
 
-[ "x$KEYS" = "x" ] && KEYS="# Step"
+[ "$KEYS" = "" ] && KEYS="# Step"
 NKEYS=$(howMany "$KEYS")
-NCOLS=$(tail -1 $FILE | awk '{ print NF}')
+NCOLS=$(grep -m 1 '^[^#]' "$FILE" | awk '{ print NF }')
 
 # With full column labels, NKEYS = NCOLS + 1, since it includes "#"
 
 # If NKEYS > NCOLS + 1, REMOVE EXCESS KEYS
-NCOLSPONE=$(expr $NCOLS + 1)
+NCOLSPONE=$((NCOLS+1))
 [ "$NKEYS" -gt "$NCOLSPONE" ] && KEYS=$(echo $KEYS | cut -d" " -f1-$NCOLSPONE)
 NKEYS=$(howMany "$KEYS")
 
 i=0
 while [ "$NKEYS" -le "$NCOLS" ]
 do
-    i=$(expr $i + 1)
+    i=$((i+1))
     KEYS="$KEYS data$i"
     NKEYS=$(howMany "$KEYS")
 done
@@ -159,33 +198,33 @@ XLABEL=$(echo $KEYS | cut -d " " -f2)
 KEYS=$(echo $KEYS | cut -d " " -f3-)
 
 GPFILE=$(mktemp)
-plotFileHeader > $GPFILE
+plotFileHeader > "$GPFILE"
 i=1
 for field in $KEYS
 do
-    i=$(expr $i + 1)
+    i=$((i+1))
     PLOTLINE="\"$FILE\" using 1:${i} with lines title \"$field\""
     if [ $i -lt $NCOLS ]
     then
        PLOTLINE="$PLOTLINE, \\"
     fi
-    echo $PLOTLINE >> $GPFILE
+    echo $PLOTLINE >> "$GPFILE"
 done
-plotFileFooter >> $GPFILE
+plotFileFooter >> "$GPFILE"
 
-touch $FILE
+touch "$FILE"
-$GNUPLOT $GPFILE &
+$GNUPLOT "$GPFILE" &
 PID=$!
 
 while true
 do
     MODTIME=$(stat --format=%Y $FILE)
-    IDLEAGO=$(expr $(date +%s) - $IDLE)
+    IDLEAGO=$(($(date +%s)-IDLE))
     test "$MODTIME" -gt "$IDLEAGO" || break
     sleep $REFRESH
 done
 
 kill -9 $PID
-rm $GPFILE
+rm -f "$GPFILE"
 
 #------------------------------------------------------------------------------

src/functionObjects/field/pressure/pressure.H

@@ -96,7 +96,7 @@ Usage
         mode            \<option\>;
 
         // Optional entries (runtime modifiable)
-        p               \<pName\>;
+        field           \<pName\>;
         U               \<UName\>;
         rho             <rhoName>;
         rhoInf          1.0; // enabled if rho=rhoInf
@@ -118,7 +118,7 @@ Usage
       type         | Type name: pressure                | word |  yes  | -
       libs         | Library name: fieldFunctionObjects | word |  yes  | -
       mode         | Calculation mode (see below)       | word |  yes  | -
-      p            | Name of the pressure field         | word |  no   | p
+      field        | Name of the pressure field         | word |  no   | p
       U            | Name of the velocity field         | word |  no   | U
       rho          | Name of the density field          | word |  no   | rho
       rhoInf | Freestream density for coefficient calculation | scalar <!--

src/lagrangian/spray/submodels/BreakupModel/PilchErdman/PilchErdman.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2013 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -57,13 +58,6 @@ Foam::PilchErdman<CloudType>::PilchErdman(const PilchErdman<CloudType>& bum)
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::PilchErdman<CloudType>::~PilchErdman()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -93,13 +87,13 @@ bool Foam::PilchErdman<CloudType>::update
 )
 {
     // Weber number - eq (1)
-    scalar We = rhoc*sqr(Urmag)*d/sigma;
+    const scalar We = rhoc*sqr(Urmag)*d/sigma;
 
     // Ohnesorge number - eq (2)
-    scalar Oh = mu/sqrt(rho*d*sigma);
+    const scalar Oh = mu/sqrt(rho*d*sigma);
 
     // Critical Weber number - eq (5)
-    scalar Wec = 12.0*(1.0 + 1.077*pow(Oh, 1.6));
+    const scalar Wec = 12.0*(1.0 + 1.077*pow(Oh, 1.6));
 
     if (We > Wec)
     {
@@ -116,7 +110,7 @@ bool Foam::PilchErdman<CloudType>::update
             else if (We > 45)
             {
                 // bag-and-stamen breakup  - eq (10)
-                taubBar = 14.1*pow(We - 12.0, 0.25);
+                taubBar = 14.1*pow(We - 12.0, -0.25);
             }
             else if (We > 18)
             {
@@ -135,15 +129,14 @@ bool Foam::PilchErdman<CloudType>::update
             }
         }
 
-        scalar rho12 = sqrt(rhoc/rho);
+        const scalar rho12 = sqrt(rhoc/rho);
 
         // velocity of fragmenting drop - eq (20)
-        scalar Vd = Urmag*rho12*(B1_*taubBar + B2_*sqr(taubBar));
+        const scalar Vd = Urmag*rho12*(B1_*taubBar + B2_*sqr(taubBar));
 
         // maximum stable diameter - eq (33)
-        scalar Vd1 = sqr(1.0 - Vd/Urmag);
+        const scalar Vd1 = max(sqr(1.0 - Vd/Urmag), SMALL);
-        Vd1 = max(Vd1, SMALL);
+        const scalar dStable = Wec*sigma/(Vd1*rhoc*sqr(Urmag));
-        scalar dStable = Wec*sigma/(Vd1*rhoc*sqr(Urmag));
 
         if (d < dStable)
         {
@@ -153,13 +146,13 @@ bool Foam::PilchErdman<CloudType>::update
         }
         else
         {
-            scalar semiMass = nParticle*pow3(d);
+            const scalar semiMass = nParticle*pow3(d);
 
             // invert eq (3) to create a dimensional break-up time
-            scalar taub = taubBar*d/(Urmag*rho12);
+            const scalar taub = taubBar*d/(Urmag*rho12);
 
             // update droplet diameter according to the rate eq (implicitly)
-            scalar frac = dt/taub;
+            const scalar frac = dt/taub;
             d = (d + frac*dStable)/(1.0 + frac);
 
             // correct the number of particles to conserve mass

src/lagrangian/spray/submodels/BreakupModel/PilchErdman/PilchErdman.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -30,33 +31,42 @@ Group
     grpLagrangianIntermediateBreakupSubModels
 
 Description
-    Particle secondary breakup model, based on the reference:
+    Particle secondary breakup model based on
+    Pilch-Erdman total droplet breakup model.
 
+    Reference:
     \verbatim
-    Pilch, M. and Erdman, C.A.
+        Pilch, M., & Erdman, C. A. (1987).
-    "Use of breakup time data and velocity history data
+        Use of breakup time data and velocity history data
-     to predict the maximum size of stable fragments for acceleration
+        to predict the maximum size of stable fragments for
-     induced breakup of a liquid drop."
+        acceleration-induced breakup of a liquid drop.
-    Int. J. Multiphase Flows 13 (1987), 741-757
+        International journal of multiphase flow, 13(6), 741-757.
+        DOI:10.1016/0301-9322(87)90063-2
     \endverbatim
 
     The droplet fragment velocity is described by the equation:
 
     \f[
-        V_d = V sqrt(epsilon)(B1 T + B2 T^2)
+        V_d = V \sqrt(\epsilon)(B_1 T + B_2 T^2)
     \f]
 
+    where
+    \vartable
+        V_d      | Fragment velocity
+        V        | Magnitude of the relative velocity
+        \epsilon | Density ratio (rho_carrier/rho_droplet)
+        T        | characteristic break-up time
+        B_1      | Model coefficient
+        B_2      | Model coefficient
+    \endvartable
     Where:
-        V_d     : fragment velocity
-        V       : magnitude of the relative velocity
-        epsilon : density ratio (rho_carrier/rho_droplet)
-        T       : characteristic break-up time
-        B1, B2  : model input coefficients
 
     The authors suggest that:
         compressible flow   : B1 = 0.75*1.0; B2 = 3*0.116
         incompressible flow : B1 = 0.75*0.5; B2 = 3*0.0758
 
+SourceFiles
+    PilchErdman.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -79,11 +89,12 @@ class PilchErdman
 :
     public BreakupModel<CloudType>
 {
-private:
+    // Private Data
-
-    // Private data
 
+        //- Model coefficient
         scalar B1_;
+
+        //- Model coefficient
         scalar B2_;
 
 
@@ -110,9 +121,12 @@ public:
             );
         }
 
+        //- No copy assignment
+        void operator=(const PilchErdman<CloudType>&) = delete;
+
 
     //- Destructor
-    virtual ~PilchErdman();
+    virtual ~PilchErdman() = default;
 
 
     // Member Functions

src/sampling/sampledSet/uniform/uniformSet.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -40,8 +41,6 @@ namespace Foam
     addToRunTimeSelectionTable(sampledSet, uniformSet, word);
 }
 
-const Foam::scalar Foam::uniformSet::tol = 1e-3;
-
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -56,14 +55,14 @@ bool Foam::uniformSet::nextSample
 {
     bool pointFound = false;
 
-    const vector normOffset = offset/mag(offset);
+    const vector normOffset(offset/mag(offset));
 
     samplePt += offset;
     ++sampleI;
 
     for (; sampleI < nPoints_; ++sampleI)
     {
-        scalar s = (samplePt - currentPt) & normOffset;
+        const scalar s = (samplePt - currentPt) & normOffset;
 
         if (s > -smallDist)
         {
@@ -92,9 +91,8 @@ bool Foam::uniformSet::trackToBoundary
 ) const
 {
     // distance vector between sampling points
-    const vector offset = (end_ - start_)/(nPoints_ - 1);
+    const vector offset((end_ - start_)/(nPoints_ - 1));
-    const vector smallVec = tol*offset;
+    const scalar smallDist = mag(tol_*offset);
-    const scalar smallDist = mag(smallVec);
 
     point trackPt = singleParticle.position();
 
@@ -205,9 +203,9 @@ void Foam::uniformSet::calcSamples
             << exit(FatalError);
     }
 
-    const vector offset = (end_ - start_)/(nPoints_ - 1);
+    const vector offset((end_ - start_)/(nPoints_ - 1));
-    const vector normOffset = offset/mag(offset);
+    const vector normOffset(offset/mag(offset));
-    const vector smallVec = tol*offset;
+    const vector smallVec(tol_*offset);
     const scalar smallDist = mag(smallVec);
 
     // Force calculation of cloud addressing on all processors
@@ -236,7 +234,7 @@ void Foam::uniformSet::calcSamples
     label trackCelli = -1;
     label trackFacei = -1;
 
-    bool isSample =
+    const bool isSample =
         getTrackingPoint
         (
             start_,
@@ -325,7 +323,7 @@ void Foam::uniformSet::calcSamples
 
         while (bHitI < bHits.size())
         {
-            scalar dist =
+            const scalar dist =
                 (bHits[bHitI].hitPoint() - singleParticle.position())
               & normOffset;
 
@@ -341,6 +339,8 @@ void Foam::uniformSet::calcSamples
             if (dist > smallDist)
             {
                 // hitpoint is past tracking position
+                bPoint = bHits[bHitI].hitPoint();
+                bFacei = bHits[bHitI].index();
                 foundValidB = true;
                 break;
             }
@@ -427,6 +427,7 @@ Foam::uniformSet::uniformSet
     sampledSet(name, mesh, searchEngine, axis),
     start_(start),
     end_(end),
+    tol_(1e-3),
     nPoints_(nPoints)
 {
     genSamples();
@@ -444,6 +445,7 @@ Foam::uniformSet::uniformSet
     sampledSet(name, mesh, searchEngine, dict),
     start_(dict.get<point>("start")),
     end_(dict.get<point>("end")),
+    tol_(dict.getCheckOrDefault<scalar>("tol", 1e-3, scalarMinMax::ge(0))),
     nPoints_(dict.get<label>("nPoints"))
 {
     genSamples();

src/sampling/sampledSet/uniform/uniformSet.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,17 +28,51 @@ Class
     Foam::uniformSet
 
 Description
+    A sampler type which provides a uniform distribution of \c nPoints
+    sample locations along a straight line specified between a given
+    \c start and an \c end points.
 
-    For a dictionary specification:
+Usage
+    Example specification:
+    \verbatim
+    sets
+    (
+        <set>
+        {
+            // Mandatory entries
+            type        uniform;
+            axis        <options>;
+            start       <vector>;
+            end         <vector>;
+            nPoints     <label>;
+
+            // Optional entries
+            tol         <scalar>;
+        }
+
+    );
+    \endverbatim
+
+    where the entries mean:
     \table
-        Property | Description                             | Required | Default
+      Property  | Description                        | Type | Reqd  | Deflt
-        type     | uniform                                 | yes      |
+      type      | Type name: uniform                 | word | yes   | -
-        axis     | x, y, z, xyz, distance                  | yes      |
+      axis      | Output type of sample locations    | word | yes   | -
-        start    | The start point                         | yes      |
+      start     | Start point of sample line         | vector | yes | -
-        end      | The end point                           | yes      |
+      end       | End point of sample line           | vector | yes | -
-        nPoints  | The number of points between start/end  | yes
+      nPoints   | Number of points between start/end | label  | yes | -
+      tol       | Relative tolerance                 | scalar | no  | 1e-3
     \endtable
 
+    Options for the \c axis entry:
+    \verbatim
+      x        | x-ordinate of a sample
+      y        | y-ordinate of a sample
+      z        | z-ordinate of a sample
+      xyz      | x-y-z coordinates of a sample
+      distance | Normal distance to the first point of a sample
+    \endverbatim
+
 SourceFiles
     uniformSet.C
 
@@ -63,17 +98,22 @@ class uniformSet
 :
     public sampledSet
 {
-    // Private data
+    // Private Data
 
-        //- Starting point
+        //- Start point of sample line
         point start_;
 
-        //- End point
+        //- End point of sample line
         point end_;
 
-        //- Number of points
+        //- Relative tolerance when comparing points
+        //- relative to difference between start_ and end_
+        scalar tol_;
+
+        //- Number of sampling points
         label nPoints_;
 
+
     // Private Member Functions
 
         //- Calculates - starting at samplePt - the first sampling point
@@ -125,13 +165,6 @@ public:
     TypeName("uniform");
 
 
-    // Static data
-
-        //- Tolerance when comparing points relative to difference between
-        //  start_ and end_
-        static const scalar tol;
-
-
     // Constructors
 
         //- Construct from components