diff --git a/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C b/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C
index 21ed696555fa6f3d560f53b9dc020d4f4c68cc68..971c546c29f3d60c51b1b157de11a2d78fa6decc 100644
--- a/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C
+++ b/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C
@@ -881,6 +881,7 @@ int main(int argc, char *argv[])
 
         // Put all modifications into meshMod
         bool anyChange = collapser.setRefinement(allPointInfo, meshMod);
+        reduce(anyChange, orOp<bool>());
 
         if (anyChange)
         {
diff --git a/applications/utilities/parallelProcessing/decomposePar/decomposePar.C b/applications/utilities/parallelProcessing/decomposePar/decomposePar.C
index 4c94937c20df47f0beb2ad3a19e2390ce2e59c5e..e6ec359d3c199c68f6e594f06f2a2a16540f1a8e 100644
--- a/applications/utilities/parallelProcessing/decomposePar/decomposePar.C
+++ b/applications/utilities/parallelProcessing/decomposePar/decomposePar.C
@@ -35,41 +35,45 @@ Description
     execution of OpenFOAM.
 
 Usage
-    \b decomposePar [OPTION]
+    \b decomposePar [OPTIONS]
 
     Options:
-      - \par -dry-run
-        Test without actually decomposing
+      - \par -allRegions
+        Decompose all regions in regionProperties. Does not check for
+        existence of processor*.
+
+      - \par - case <dir>
+        Specify case directory to use (instead of the cwd).
 
       - \par -cellDist
         Write the cell distribution as a labelList, for use with 'manual'
         decomposition method and as a volScalarField for visualization.
 
-      - \par -region \<regionName\>
-        Decompose named region. Does not check for existence of processor*.
+      - \par -constant
+        Include the 'constant/' dir in the times list.
 
-      - \par -allRegions
-        Decompose all regions in regionProperties. Does not check for
-        existence of processor*.
+      - \par -copyUniform
+        Copy any \a uniform directories too.
 
       - \par -copyZero
         Copy \a 0 directory to processor* rather than decompose the fields.
 
-      - \par -copyUniform
-        Copy any \a uniform directories too.
+      - \par -debug-switch <name=val>
+        Specify the value of a registered debug switch. Default is 1
+        if the value is omitted. (Can be used multiple times)
 
-      - \par -constant
+      - \par -decomposeParDict <file>
+        Use specified file for decomposePar dictionary.
 
-      - \par -time xxx:yyy
-        Override controlDict settings and decompose selected times. Does not
-        re-decompose the mesh i.e. does not handle moving mesh or changing
-        mesh cases.
+      - \par -dry-run
+        Test without writing the decomposition. Changes -cellDist to
+        only write volScalarField.
 
       - \par -fields
         Use existing geometry decomposition and convert fields only.
 
-      - \par -noSets
-        Skip decomposing cellSets, faceSets, pointSets.
+      - \par fileHandler <handler>
+        Override the file handler type.
 
       - \par -force
         Remove any existing \a processor subdirectories before decomposing the
@@ -83,6 +87,58 @@ Usage
         be used with caution when the underlying (serial) geometry or the
         decomposition method etc. have been changed between decompositions.
 
+      - \par -info-switch <name=val>
+        Specify the value of a registered info switch. Default is 1
+        if the value is omitted. (Can be used multiple times)
+
+      - \par -latestTime
+        Select the latest time.
+
+      - \par -lib <name>
+        Additional library or library list to load (can be used multiple times).
+
+      - \par -noFunctionObjects
+        Do not execute function objects.
+
+      - \par -noSets
+        Skip decomposing cellSets, faceSets, pointSets.
+
+      - \par -noZero
+        Exclude the \a 0 dir from the times list.
+
+      - \par -opt-switch <name=val>
+        Specify the value of a registered optimisation switch (int/bool).
+        Default is 1 if the value is omitted. (Can be used multiple times)
+
+      - \par -region \<regionName\>
+        Decompose named region. Does not check for existence of processor*.
+
+      - \par -time <ranges>
+        Override controlDict settings and decompose selected times. Does not
+        re-decompose the mesh i.e. does not handle moving mesh or changing
+        mesh cases. Eg, ':10,20 40:70 1000:', 'none', etc.
+
+      - \par -verbose
+        Additional verbosity.
+
+      - \par -doc
+        Display documentation in browser.
+
+      - \par -doc-source
+        Display source code in browser.
+
+      - \par -help
+        Display short help and exit.
+
+      - \par -help-man
+        Display full help (manpage format) and exit.
+
+      - \par -help-notes
+        Display help notes (description) and exit.
+
+      - \par -help-full
+        Display full help and exit.
+
 \*---------------------------------------------------------------------------*/
 
 #include "OSspecific.H"
@@ -328,10 +384,12 @@ int main(int argc, char *argv[])
         times = timeSelector::selectIfPresent(runTime, args);
     }
 
-
     // Allow override of decomposeParDict location
-    const fileName decompDictFile =
-        args.getOrDefault<fileName>("decomposeParDict", "");
+    fileName decompDictFile(args.get<fileName>("decomposeParDict", ""));
+    if (!decompDictFile.empty() && !decompDictFile.isAbsolute())
+    {
+        decompDictFile = runTime.globalPath()/decompDictFile;
+    }
 
     // Get all region names
     wordList regionNames;
diff --git a/etc/caseDicts/annotated/blockMeshDict b/etc/caseDicts/annotated/blockMeshDict
new file mode 100644
index 0000000000000000000000000000000000000000..c088a1fd168a98711120c5ce20d1173f80bcc460
--- /dev/null
+++ b/etc/caseDicts/annotated/blockMeshDict
@@ -0,0 +1,96 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+scale    1;
+
+
+x0       0;
+x1       #eval "$x0 + 1.0";
+y0       0;
+y1       1;
+z0       0;
+z1       1;
+
+nx       2;
+ny       2;
+nz       1;
+
+gr       40.0;
+grInv    #eval "1/$gr";
+
+vertices
+(
+    ($x0 $y0 $z0)
+    ($x1 $y0 $z0)
+    ($x1 $y1 $z0)
+    ($x0 $y1 $z0)
+    ($x0 $y0 $z1)
+    ($x1 $y0 $z1)
+    ($x1 $y1 $z1)
+    ($x0 $y1 $z1)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) ($nx $ny $nz) simpleGrading (1 $grInv 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    upperWall
+    {
+        type symmetry;
+        faces
+        (
+            (3 7 6 2)
+        );
+    }
+    lowerWall
+    {
+        type wall;
+        faces
+        (
+            (1 5 4 0)
+        );
+    }
+    inlet
+    {
+        type patch;
+        faces
+        (
+            (0 4 7 3)
+        );
+    }
+    outlet
+    {
+        type patch;
+        faces
+        (
+            (2 6 5 1)
+        );
+    }
+);
+
+mergePatchPairs
+(
+);
+
+
+// ************************************************************************* //
diff --git a/etc/caseDicts/annotated/boxTurbDict b/etc/caseDicts/annotated/boxTurbDict
new file mode 100644
index 0000000000000000000000000000000000000000..b42bfc0388b9021f23da2dae93ce39bc2fcee086
--- /dev/null
+++ b/etc/caseDicts/annotated/boxTurbDict
@@ -0,0 +1,23 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      boxTurbDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Ea              10;
+
+k0              5;
+
+
+// ************************************************************************* //
diff --git a/tutorials/verificationAndValidation/turbulentInflow/system/sampling b/etc/caseDicts/annotated/createBoxTurbDict
similarity index 53%
rename from tutorials/verificationAndValidation/turbulentInflow/system/sampling
rename to etc/caseDicts/annotated/createBoxTurbDict
index ef34447e5a8c7c44d3b0836dbec5ae450b96e662..35b426f065698815c96e300b12be29a98ea4413d 100644
--- a/tutorials/verificationAndValidation/turbulentInflow/system/sampling
+++ b/etc/caseDicts/annotated/createBoxTurbDict
@@ -10,39 +10,43 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "system";
-    object      sampling;
+    location    "constant";
+    object      createBoxTurbDict;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-inletSampling
-{
-    type                sets;
-    libs                (sampling);
-    writeControl        writeTime;
-    timeStart           $/timeStart;
+N           (64 64 64);
 
-    interpolationScheme cellPoint;
-    setFormat           raw;
-    fields              (UPrime2Mean);
+// Suggested box size of 9*2*pi [cm]
+L           (0.56548667765 0.56548667765 0.56548667765);
+
+nModes      5000;
+
+// Energy as a function of wave number
+// Here using Comte-Bellot and Corrsin data at t.U_0/M = 42 (see Ref. table 3)
+Ek          table
+(
+    (15 0)
+    (20 0.000129)
+    (25 0.00023)
+    (30 0.000322)
+    (40 0.000435)
+    (50 0.000457)
+    (70 0.00038)
+    (100 0.00027)
+    (150 0.000168)
+    (200 0.00012)
+    (250 8.9e-05)
+    (300 7.03e-05)
+    (400 4.7e-05)
+    (600 2.47e-05)
+    (800 1.26e-05)
+    (1000 7.42e-06)
+    (1250 3.96e-06)
+    (1500 2.33e-06)
+    (1750 1.34e-06)
+    (2000 8e-07)
+);
 
-    sets
-    (
-        inletPatch
-        {
-            type        face;
-            axis        y;
-            start       (0.0 0 1.57);
-            end         (0.0 2 1.57);
-        }
-        inletCell
-        {
-            type        midPoint;
-            axis        y;
-            start       (0.062832 0 1.57);
-            end         (0.062832 2 1.57);
-        }
-    );
-}
 
 // ************************************************************************* //
diff --git a/etc/caseDicts/annotated/dsmcInitialiseDict b/etc/caseDicts/annotated/dsmcInitialiseDict
new file mode 100644
index 0000000000000000000000000000000000000000..d37723266528be6c840d59a8cfa94e47f4a9f1a3
--- /dev/null
+++ b/etc/caseDicts/annotated/dsmcInitialiseDict
@@ -0,0 +1,29 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      dsmcInitialiseDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberDensities
+{
+    N2          0.777e20;
+    O2          0.223e20;
+};
+
+temperature     300;
+
+velocity        (1325 -352 823);
+
+
+// ************************************************************************* //
diff --git a/etc/caseDicts/annotated/extrudeMeshDict b/etc/caseDicts/annotated/extrudeMeshDict
index a5d944d885ec6fb7a9f993df4f6aec702cac668f..81742b86ab4a946eb0145a41d45799d71d144e6e 100644
--- a/etc/caseDicts/annotated/extrudeMeshDict
+++ b/etc/caseDicts/annotated/extrudeMeshDict
@@ -24,7 +24,7 @@ constructFrom patch;
 //constructFrom surface;
 
 // If construct from patch/mesh:
-sourceCase "../cavity";
+sourceCase    "$FOAM_CASE";
 sourcePatches (movingWall);
 
 // If construct from patch: patch to use for back (can be same as sourcePatch)
diff --git a/etc/caseDicts/annotated/mdEquilibrationDict b/etc/caseDicts/annotated/mdEquilibrationDict
new file mode 100644
index 0000000000000000000000000000000000000000..6c2a3a4fc8c92958390ad3e81414d694fabc26c8
--- /dev/null
+++ b/etc/caseDicts/annotated/mdEquilibrationDict
@@ -0,0 +1,20 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      mdEquilibrationDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+targetTemperature  298;
+
+
+// ************************************************************************* //
diff --git a/etc/caseDicts/annotated/mdInitialiseDict b/etc/caseDicts/annotated/mdInitialiseDict
new file mode 100644
index 0000000000000000000000000000000000000000..f56f236ecfe1cebdf1ef3ab15d14941daba7ba50
--- /dev/null
+++ b/etc/caseDicts/annotated/mdInitialiseDict
@@ -0,0 +1,82 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      mdInitialiseDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Euler angles, expressed in degrees as phi, theta, psi, see
+// http://mathworld.wolfram.com/EulerAngles.html
+
+sectionA
+{
+    massDensity             1004;
+    temperature             298;
+    bulkVelocity            (0.0 0.0 0.0);
+    latticeIds
+    (
+        water
+    );
+    tetherSiteIds           ();
+    latticePositions
+    (
+        (0 0 0)
+    );
+    anchor                  (0 0 0);
+    orientationAngles       (0 0 0);
+    latticeCellShape        (1 1 1);
+}
+
+sectionB
+{
+    massDensity             1004;
+    temperature             298;
+    bulkVelocity            (0.0 0.0 0.0);
+    latticeIds
+    (
+        Ar
+    );
+    tetherSiteIds           ();
+    latticePositions
+    (
+        (0 0 0)
+    );
+    anchor                  (0 0 0);
+    orientationAngles       (0 0 0);
+    latticeCellShape        (1 1 1);
+}
+
+sectionC
+{
+    massDensity             1004;
+    temperature             298;
+    bulkVelocity            (0.0 0.0 0.0);
+    latticeIds
+    (
+        water1
+        water2
+    );
+    tetherSiteIds           ();
+    latticePositions
+    (
+        (0 0 0)
+        (0 0.5 0.5)
+        (0.5 0 0.5)
+        (0.5 0.5 0)
+    );
+    anchor                  (0 0 0);
+    orientationAngles       (0 0 0);
+    latticeCellShape        (1 1 1);
+}
+
+
+// ************************************************************************* //
diff --git a/etc/caseDicts/annotated/obstaclesDict b/etc/caseDicts/annotated/obstaclesDict
new file mode 100644
index 0000000000000000000000000000000000000000..87c2ee10d9f3f3875b33cbb74e8b886d67751a09
--- /dev/null
+++ b/etc/caseDicts/annotated/obstaclesDict
@@ -0,0 +1,39 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      obstaclesDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+scale   1.0;
+
+verbose 0;
+
+#include "groups/group4"
+
+_01
+{
+    groupId 1;
+
+    zpipe { direction z; length 0.947; diameter 0.026; }
+
+    obstacles
+    (
+        box { point (0 0 0); size (0.05 0.05 2); }
+        box { point (1 0 0); size (0.05 0.05 2); }
+        box { point (2 0 0); size (0.05 0.05 2); }
+        cyl { point (1.031   0.975   0.056); $zpipe; }
+    );
+}
+
+
+// ************************************************************************* //
diff --git a/etc/caseDicts/annotated/optimisationDict b/etc/caseDicts/annotated/optimisationDict
new file mode 100644
index 0000000000000000000000000000000000000000..9887a5921bee194084b6a6a836b02d7444b8274c
--- /dev/null
+++ b/etc/caseDicts/annotated/optimisationDict
@@ -0,0 +1,393 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      optimisationDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+optimisationManager     steadyOptimisation;  // singleRun;
+
+
+primalSolvers
+{
+    p1
+    {
+        active                 true;
+        type                   incompressible;
+        solver                 simple;
+
+        solutionControls
+        {
+            nIters 3000;
+            residualControl
+            {
+                "p.*"       1.e-7;
+                "U.*"       1.e-7;
+            }
+        }
+    }
+}
+
+
+adjointManagers
+{
+    am1
+    {
+        primalSolver             p1;
+        adjointSolvers
+        {
+            as1
+            {
+                // choose adjoint solver
+                //----------------------
+                active                 true;
+                type                   incompressible;
+                solver                 adjointSimple;
+                useSolverNameForFields true;
+
+                // manage objectives
+                //------------------
+                objectives
+                {
+                    type                incompressible;
+                    objectiveNames
+                    {
+                        losses
+                        {
+                            weight          1;
+                            type            PtLosses;
+                            patches         (Inlet Outlet);
+                        }
+                    }
+                }
+                // ATC treatment
+                //--------------
+                ATCModel
+                {
+                    ATCModel        standard;
+                }
+                // solution control
+                //------------------
+                solutionControls
+                {
+                    nIters 3000;
+                    residualControl
+                    {
+                        "pa.*"       1.e-7;
+                        "Ua.*"       1.e-7;
+                    }
+                }
+            }
+            vol
+            {
+                // choose adjoint solver
+                //----------------------
+                active                 true;
+                type                   incompressible;
+                solver                 adjointSimple;
+                useSolverNameForFields true;
+                isConstraint           true;
+                // manage objectives
+                //------------------
+                objectives
+                {
+                    type                incompressible;
+
+                    objectiveNames
+                    {
+                        vol
+                        {
+                            weight              1;
+                            type                partialVolume;
+                            patches             (lower upper);
+                        }
+                    }
+                }
+
+                // ATC treatment
+                //--------------
+                ATCModel
+                {
+                    ATCModel        standard;
+                }
+
+                // solution control
+                //------------------
+                solutionControls
+                {
+                    nIters 3000;
+                    residualControl
+                    {
+                        "pa.*"       1.e-7;
+                        "Ua.*"       1.e-7;
+                    }
+                }
+            }
+        }
+    }
+}
+
+
+optimisation
+{
+    optimisationType
+    {
+        type             shapeOptimisation;
+        writeEachMesh    true;
+    }
+
+    sensitivities
+    {
+        type            volumetricBSplinesFI;
+        patches         (lower upper);
+    }
+
+    updateMethod
+    {
+        method  SQP;
+        SQP
+        {
+            etaHessian        0.8;
+            nSteepestDescent  1;
+            scaleFirstHessian true;
+        }
+    }
+
+    meshMovement
+    {
+        type                   volumetricBSplines;
+        maxAllowedDisplacement 2.e-3;
+    }
+
+/*
+    updateMethod
+    {
+        method    BFGS;
+        BFGS
+        {
+            etaHessian        0.8;
+            nSteepestDescent  1;
+            scaleFirstHessian true;
+            activeDesignVariables
+            (
+                141 142 144 145 147 148
+            );
+        }
+    }
+*/
+
+/*
+    sensitivities
+    {
+        type            multiple; // used to compute many kinds of sensitivities at the same time
+        patches         (lower upper);
+        sensTypes
+        {
+            FIVolSplines
+            {
+                type                volumetricBSplinesFI;
+                patches             (lower upper);
+                includeDistance     true;
+                adjointEikonalSolver
+                {
+                    iters  1000;
+                    tolerance 1.e-6;
+                }
+            }
+            ESIVolSplines
+            {
+                type                         volumetricBSplines;
+                patches                      (lower upper);
+                includeObjectiveContribution true; // one of this or the equivalent flag in
+                                                   // surfaceSensitivities has to be set to true
+                                                   // with this being the prefered one
+                surfaceSensitivities
+                {
+                    patches                      (lower upper);
+                    includeSurfaceArea           true;
+                    includeMeshMovement          true;
+                    includeDistance              true;
+                    includeObjectiveContribution false;
+
+                    // adjointEikonal and adjointMeshMovement solvers should be always nested
+                    // within the dictionary of the sensitivity type they correspond to.
+                    // For (E)SI based sensitivities, this means the surfaceSensitivities dict
+                    // Default values are provided, so the dictionaries can be skipped
+                    adjointEikonalSolver
+                    {
+                        iters  1000;
+                        tolerance 1.e-6;
+                    }
+                    adjointMeshMovementSolver
+                    {
+                        iters  10000;
+                        tolerance 1.e-6;
+                    }
+                }
+            }
+            SIVolSplines
+            {
+                type                         volumetricBSplines;
+                patches                      (lower upper);
+                includeObjectiveContribution true; // same comment as above
+                surfaceSensitivities
+                {
+                    patches                      (lower upper);
+                    includeSurfaceArea           true;
+                    includeMeshMovement          false;
+                    includeDistance              true;
+                    includeObjectiveContribution false;
+                    adjointEikonalSolver
+                    {
+                        iters  1000;
+                        tolerance 1.e-6;
+                    }
+                }
+            }
+            FIBezier
+            {
+                type            BezierFI;
+                includeDistance true;
+                patches         (lower upper);
+                dxdbSolver
+                {
+                    iters           1000;
+                    tolerance       1.e-6;
+                }
+                adjointEikonalSolver
+                {
+                    iters  1000;
+                    tolerance 1.e-6;
+                }
+            }
+            ESIBezier
+            {
+                type            Bezier;
+                includeObjectiveContribution true; // same comment as above
+                surfaceSensitivities
+                {
+                    patches             (lower upper);
+                    includeSurfaceArea  true;
+                    includeMeshMovement true;
+                    includeDistance     true;
+                    includeObjectiveContribution false;
+                    adjointEikonalSolver
+                    {
+                        iters  1000;
+                        tolerance 1.e-6;
+                    }
+                    adjointMeshMovementSolver
+                    {
+                        iters  10000;
+                        tolerance 1.e-6;
+                    }
+                }
+                patches         (lower upper);
+            }
+            SIBezier
+            {
+                type                Bezier;
+                includeObjectiveContribution true; // same comment as above
+                surfaceSensitivities
+                {
+                    patches             (lower upper);
+                    includeSurfaceArea  true;
+                    includeMeshMovement false;
+                    includeDistance     true;
+                    includeObjectiveContribution false;
+                    adjointEikonalSolver
+                    {
+                        iters  1000;
+                        tolerance 1.e-6;
+                    }
+                }
+                patches             (lower upper);
+            }
+        }
+    }
+*/
+
+/*
+    sensitivities
+    {
+        type            multiple; // used to compute many kinds of sensitivities at the same time
+        patches         (lower upper);
+        sensTypes
+        {
+            FIVolSplines
+            {
+                type                volumetricBSplinesFI;
+                patches             (lower upper);
+            }
+            ESIVolSplines
+            {
+                type                volumetricBSplines;
+                patches             (lower upper);
+            }
+            SIVolSplines
+            {
+                type                volumetricBSplines;
+                patches             (lower upper);
+                surfaceSensitivities
+                {
+                    includeMeshMovement false;
+                }
+            }
+            FIBezier
+            {
+                type                BezierFI;
+                patches             (lower upper);
+            }
+            ESIBezier
+            {
+                type                Bezier;
+                patches             (lower upper);
+            }
+            SIBezier
+            {
+                type                Bezier;
+                patches             (lower upper);
+                surfaceSensitivities
+                {
+                    includeMeshMovement false;
+                }
+            }
+        }
+    }
+*/
+}
+
+
+/*
+Bezier
+{
+    nBezier 24;
+    confineXmovement
+    (
+        true false false false false false false false false false false true
+        true false false false false false false false false false false true
+    );
+    confineYmovement
+    (
+        true false false false false false false false false false false true
+        true false false false false false false false false false false true
+    );
+    confineZmovement
+    (
+        true true true true true true true true true true true true
+        true true true true true true true true true true true true
+    );
+}
+*/
+
+
+// ************************************************************************* //
diff --git a/etc/caseDicts/annotated/potentialDict b/etc/caseDicts/annotated/potentialDict
new file mode 100644
index 0000000000000000000000000000000000000000..71552faf68331c3ca9b5fab5a9d218a0daa2fe4a
--- /dev/null
+++ b/etc/caseDicts/annotated/potentialDict
@@ -0,0 +1,118 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      potentialDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Subdictionaries specifying types of intermolecular potential.
+// Sub-sub dictionaries specify the potentials themselves.
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Removal order
+
+// This is the order in which to remove overlapping pairs if more than one
+// type of molecule is present.  The most valuable molecule type is at the
+// right hand end, the molecule that will be removed 1st is 1st on the list.
+// Not all types need to be present, a molecule that is not present is
+// automatically less valuable than any on the list.  For molecules of the
+// same type there is no control over which is removed.
+
+removalOrder ( water );
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Potential Energy Limit
+
+// Maximum permissible pair energy allowed at startup.  Used to remove
+// overlapping molecules created during preprocessing.
+
+potentialEnergyLimit 1e-18;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Pair potentials
+
+// If a pair are not present here it is assumed that they do not interact.
+
+// Electrostatic pair interactions are not listed here - they are handled
+// separately.
+
+// If there are r different type of molecules, and a pair force is required
+// between all combinations, then there are C = r(r+1)/2 combinations,
+// i.e. for r = {1,2,3,4}, C = {1,3,6,10} (sum of triangular numbers).
+
+// Pair potentials are specified by the combination of their ids,
+// for MOLA and MOLB, "MOLA-MOLB" OR "MOLB-MOLA" is acceptable
+// (strictly OR, both or neither will throw an error)
+
+pair
+{
+    O-O
+    {
+        pairPotential   lennardJones;
+        rCut            1.0e-9;
+        rMin            0.1e-9;
+        dr              1e-13;
+        lennardJonesCoeffs
+        {
+            sigma       3.154e-10;
+            epsilon     1.07690722e-21;
+        }
+        energyScalingFunction   noScaling;
+        writeTables     yes;
+    }
+
+    electrostatic
+    {
+        pairPotential   dampedCoulomb;
+        rCut            1e-9;
+        rMin            2e-11;
+        dr              2e-12;
+        dampedCoulombCoeffs
+        {
+            alpha       2e9;
+        }
+        energyScalingFunction   shiftedForce;
+        writeTables     yes;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Tethering Potentials
+
+tether
+{
+    O
+    {
+        tetherPotential restrainedHarmonicSpring;
+        restrainedHarmonicSpringCoeffs
+        {
+            springConstant  0.277;
+            rR              1.2e-9;
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// External Forces
+
+// Bulk external forces (namely gravity) will be specified as forces rather
+// than potentials to allow their direction to be controlled.
+
+external
+{
+    gravity             (0 0 0);
+}
+
+
+// ************************************************************************* //
diff --git a/etc/caseDicts/annotated/probesDict b/etc/caseDicts/annotated/probesDict
new file mode 100644
index 0000000000000000000000000000000000000000..81f8194456f6293d2f6f304aab97e8c2a04a97c3
--- /dev/null
+++ b/etc/caseDicts/annotated/probesDict
@@ -0,0 +1,31 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version         2.0;
+    format          ascii;
+    class           dictionary;
+    location        system;
+    object          probesDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Fields to be probed (runTime modifiable)
+fields
+(
+    T H2O p kT
+);
+
+// Locations to be probed (runTime modifiable)
+probeLocations
+(
+    (0.005 0.0 0.0)
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/etc/caseDicts/annotated/sampleDict b/etc/caseDicts/annotated/sampleDict
new file mode 100644
index 0000000000000000000000000000000000000000..8695d9983a26e075b5caaa6010cd1fb4d1a0569e
--- /dev/null
+++ b/etc/caseDicts/annotated/sampleDict
@@ -0,0 +1,325 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      sampleDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Set output format : choice of
+//      xmgr
+//      jplot
+//      gnuplot
+//      raw
+//      vtk
+//      ensight
+//      csv
+setFormat raw;
+
+// Surface output format. Choice of
+//      null        : suppress output
+//      ensight     : Ensight Gold format, one field per case file
+//      foamFile    : separate points, faces and values file
+//      dx          : DX scalar or vector format
+//      vtk         : VTK ascii format
+//      raw         : x y z value format for use with e.g. gnuplot 'splot'.
+//
+// Note:
+// other formats such as obj, stl, etc can also be written (by proxy)
+// but without any values!
+surfaceFormat vtk;
+
+// optionally define extra controls for the output formats
+formatOptions
+{
+    ensight
+    {
+        format  ascii;
+    }
+}
+
+// interpolationScheme. choice of
+//      cell          : use cell-centre value only; constant over cells
+//                      (default)
+//      cellPoint     : use cell-centre and vertex values
+//      cellPointFace : use cell-centre, vertex and face values.
+//      pointMVC      : use point values only (Mean Value Coordinates)
+//      cellPatchConstrained : like 'cell' but uses cell-centre except on
+//                             boundary faces where it uses the boundary value.
+//                             For use with e.g. patchCloudSet.
+// 1] vertex values determined from neighbouring cell-centre values
+// 2] face values determined using the current face interpolation scheme
+//    for the field (linear, gamma, etc.)
+interpolationScheme cellPoint;
+
+// Fields to sample.
+fields
+(
+    p
+    U
+);
+
+// Set sampling definition: choice of
+//      uniform             evenly distributed points on line
+//      face                one point per face intersection
+//      midPoint            one point per cell, inbetween two face intersections
+//      midPointAndFace     combination of face and midPoint
+//
+//      polyLine            specified points, not nessecary on line, uses
+//                          tracking
+//      cloud               specified points, uses findCell
+//      triSurfaceMeshPointSet  points of triSurface
+//
+// axis: how to write point coordinate. Choice of
+// - x/y/z: x/y/z coordinate only
+// - xyz: three columns
+//  (probably does not make sense for anything but raw)
+// - distance: distance from start of sampling line (if uses line) or
+//             distance from first specified sampling point
+//
+// type specific:
+//      uniform, face, midPoint, midPointAndFace : start and end coordinate
+//      uniform: extra number of sampling points
+//      polyLine, cloud: list of coordinates
+//      patchCloud: list of coordinates and set of patches to look for nearest
+//      patchSeed: random sampling on set of patches. Points slightly off
+//                 face centre.
+sets
+(
+    lineX1
+    {
+        type        uniform;
+        axis        distance;
+
+        //- cavity. Slightly perturbed so not to align with face or edge.
+        start       (0.0201 0.05101 0.00501);
+        end         (0.0601 0.05101 0.00501);
+        nPoints     10;
+    }
+
+    lineX2
+    {
+        type        face;
+        axis        x;
+
+        //- cavity
+        start       (0.0001 0.0525 0.00501);
+        end         (0.0999 0.0525 0.00501);
+    }
+
+    somePoints
+    {
+        type    cloud;
+        axis    xyz;
+        points  ((0.049 0.049 0.00501)(0.051 0.049 0.00501));
+    }
+
+    somePatchPoints
+    {
+        // Sample nearest points on selected patches. Looks only up to
+        // maxDistance away. Any sampling point not found will get value
+        // pTraits<Type>::max (usually VGREAT)
+        // Use with interpolations:
+        // - cell (cell value)
+        // - cellPatchConstrained (boundary value)
+        // - cellPoint (interpolated boundary value)
+        type        patchCloud;
+        axis        xyz;
+        points      ((0.049 0.099 0.005)(0.051 0.054 0.005));
+        maxDistance 0.1;    // maximum distance to search
+        patches     (".*Wall.*");
+    }
+
+    patchSeed
+    {
+        type        patchSeed;
+        axis        xyz;
+        patches     (".*Wall.*");
+        // Number of points to seed. Divided amongst all processors according
+        // to fraction of patches they hold.
+        maxPoints   100;
+    }
+
+);
+
+
+// Surface sampling definition
+//
+// 1] patches are not triangulated by default
+// 2] planes are always triangulated
+// 3] iso-surfaces are always triangulated
+surfaces
+(
+    constantPlane
+    {
+        type            plane;    // always triangulated
+        basePoint       (0.0501 0.0501 0.005);
+        normalVector    (0.1 0.1 1);
+
+        //- Optional: restrict to a particular zone
+        // zone        zone1;
+
+        //- Optional: do not triangulate (only for surfaceFormats that support
+        //            polygons)
+        //triangulate     false;
+    }
+
+    interpolatedPlane
+    {
+        type            plane;    // always triangulated
+
+        // Make plane relative to the coordinateSystem (Cartesian)
+        coordinateSystem
+        {
+            origin      (0.0501 0.0501 0.005);
+
+            // Add a coordinate rotation
+            // (required, so here one that doesn't change anything)
+            coordinateRotation
+            {
+                type    axesRotation;
+                e1      (1 0 0);
+                e2      (0 1 0);
+            }
+        }
+        basePoint       (0 0 0);
+        normalVector    (0.1 0.1 1);
+        interpolate     true;
+    }
+
+    walls_constant
+    {
+        type            patch;
+        patches         ( ".*Wall.*" );
+        // Optional: whether to leave as faces (=default) or triangulate
+        // triangulate     false;
+    }
+
+    walls_interpolated
+    {
+        type            patch;
+        patches         ( ".*Wall.*" );
+        interpolate     true;
+        // Optional: whether to leave as faces (=default) or triangulate
+        // triangulate     false;
+    }
+
+    nearWalls_interpolated
+    {
+        // Sample cell values off patch. Does not need to be the near-wall
+        // cell, can be arbitrarily far away.
+        type            patchInternalField;
+        patches         ( ".*Wall.*" );
+        interpolate     true;
+
+
+        // Optional: specify how to obtain sampling points from the patch
+        //           face centres (default is 'normal')
+        //
+        //  //- Specify distance to offset in normal direction
+        offsetMode  normal;
+        distance    0.1;
+        //
+        //  //- Specify single uniform offset
+        //  offsetMode  uniform;
+        //  offset      (0 0 0.0001);
+        //
+        //  //- Specify offset per patch face
+        //  offsetMode  nonuniform;
+        //  offsets     ((0 0 0.0001) (0 0 0.0002));
+
+
+        // Optional: whether to leave as faces (=default) or triangulate
+        // triangulate     false;
+    }
+
+    interpolatedIso
+    {
+        // Iso surface for interpolated values only
+        type            isoSurface;    // always triangulated
+        isoField        rho;
+        isoValue        0.5;
+        interpolate     true;
+
+        //zone            ABC;          // Optional: zone only
+        //exposedPatchName fixedWalls;  // Optional: zone only
+
+        // regularise      false;    // Optional: do not simplify
+        // mergeTol        1e-10;    // Optional: fraction of mesh bounding box
+                                     // to merge points (default=1e-6)
+    }
+    constantIso
+    {
+        // Iso surface for constant values.
+        // Triangles guaranteed not to cross cells.
+        type            isoSurfaceCell;    // always triangulated
+        isoField        rho;
+        isoValue        0.5;
+        interpolate     false;
+        regularise      false;              // do not simplify
+        // mergeTol        1e-10;    // Optional: fraction of mesh bounding box
+                                     // to merge points (default=1e-6)
+    }
+
+    triangleCut
+    {
+        // Cutingplane using iso surface
+        type            cuttingPlane;
+        planeType       pointAndNormal;
+        pointAndNormalDict
+        {
+            basePoint       (0.4 0 0.4);
+            normalVector    (1 0.2 0.2);
+        }
+        interpolate     true;
+
+        //zone            ABC;          // Optional: zone only
+        //exposedPatchName fixedWalls;  // Optional: zone only
+
+        // regularise      false;    // Optional: do not simplify
+        // mergeTol        1e-10;    // Optional: fraction of mesh bounding box
+                                     // to merge points (default=1e-6)
+    }
+
+    distance
+    {
+        // Isosurface from signed/unsigned distance to surface
+        type            distanceSurface;
+        signed          true;
+
+        // Definition of surface
+        surfaceType     triSurfaceMesh;
+        surfaceName     integrationPlane.stl;
+        // Distance to surface
+        distance        0.0;
+
+        //cell            false;// optional: use isoSurface instead
+                                // of isoSurfaceCell
+        interpolate     false;
+        regularise      false;  // Optional: do not simplify
+        // mergeTol 1e-10;      // Optional: fraction of mesh bounding box
+                                // to merge points (default=1e-6)
+    }
+
+    triSurfaceSampling
+    {
+        // Sampling on triSurface
+        type        sampledTriSurfaceMesh;
+        surface     integrationPlane.stl;
+        source      boundaryFaces;  // What to sample: cells (nearest cell)
+                                    // insideCells (only triangles inside cell)
+                                    // boundaryFaces (nearest boundary face)
+        interpolate true;
+    }
+);
+
+
+// *********************************************************************** //
diff --git a/tutorials/verificationAndValidation/turbulentInflow/system/fieldAverage b/etc/caseDicts/annotated/setExprBoundaryFieldsDict
similarity index 71%
rename from tutorials/verificationAndValidation/turbulentInflow/system/fieldAverage
rename to etc/caseDicts/annotated/setExprBoundaryFieldsDict
index 2af894e3917d6ff65a1b6abb7b3cd7d150dd6fab..6412e4f64b29a30fd4887b49dc0b4e68119781e2 100644
--- a/tutorials/verificationAndValidation/turbulentInflow/system/fieldAverage
+++ b/etc/caseDicts/annotated/setExprBoundaryFieldsDict
@@ -10,27 +10,23 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "system";
-    object      fieldAverage;
+    object      setExprBoundaryFieldsDict;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-fieldAverage1
+pattern
 {
-    type                fieldAverage;
-    libs                (fieldFunctionObjects);
-    writeControl        writeTime;
-    timeStart           $/timeStart;
+    field   T;
 
-    fields
+    expressions
     (
-        U
         {
-            mean        on;
-            prime2Mean  on;
-            base        time;
+            patch   outlet2;
+            target  something;
+            expression #{ (pos().x() < 1e-4 ? 60 : 120) #};
         }
     );
 }
 
+
 // ************************************************************************* //
diff --git a/tutorials/verificationAndValidation/turbulentInflow/system/controlDict.template b/etc/caseDicts/annotated/setExprFieldsDict
similarity index 56%
rename from tutorials/verificationAndValidation/turbulentInflow/system/controlDict.template
rename to etc/caseDicts/annotated/setExprFieldsDict
index 91209ef2c3316b26e7c1fd7b82dfbd84d943cb55..2e06d9a26e68afb3aa9611eb1ae267a712b10dd1 100644
--- a/tutorials/verificationAndValidation/turbulentInflow/system/controlDict.template
+++ b/etc/caseDicts/annotated/setExprFieldsDict
@@ -10,50 +10,43 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "system";
-    object      controlDict;
+    object      setExprFieldsDict;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-application     pimpleFoam;
+expressions
+(
+    T
+    {
+        field       T;
+        dimensions  [0 0 0 1 0 0 0];
+
+        constants
+        {
+            centre (0.21 0 0.01);
+        }
+
+        variables
+        (
+            "radius = 0.1"
+        );
+
+        condition
+        #{
+            // Within the radius
+            (mag(pos() - $[(vector)constants.centre]) < radius)
+
+            // but only +ve y!
+          && pos((pos() - $[(vector)constants.centre]).y()) > 0
+        #};
+
+        expression
+        #{
+            300
+          + 200 * (1 - mag(pos() - $[(vector)constants.centre]) / radius)
+        #};
+    }
+);
 
-startFrom       latestTime;
-
-startTime       0;
-
-stopAt          endTime;
-
-endTime         END_TIME;
-
-deltaT          4e-3;
-
-writeControl    timeStep;
-
-writeInterval   1250;
-
-purgeWrite      3;
-
-writeFormat     ascii;
-
-writePrecision  8;
-
-writeCompression off;
-
-timeFormat      general;
-
-timePrecision   8;
-
-runTimeModifiable false;
-
-adjustTimeStep  false;
-
-// Allow 10% of time for initialisation before sampling
-timeStart       #eval #{ 0.1 * ${/endTime} #};
-
-functions
-{
-    #include "fieldAverage"
-    #include "sampling"
-}
 
 // ************************************************************************* //
diff --git a/etc/caseDicts/annotated/topoSetSourcesDict b/etc/caseDicts/annotated/topoSetSourcesDict
index 8ba14c4bc9a5bb925ad7e8db64d1f30580b067d7..58bb81929da906a744762a0645df4af03a490f51 100644
--- a/etc/caseDicts/annotated/topoSetSourcesDict
+++ b/etc/caseDicts/annotated/topoSetSourcesDict
@@ -241,6 +241,12 @@ cellSet_doc
         //     name    ".*Zone";
         // }
     }
+
+    //- Cells attached to the outside of the input cellSet
+    {
+        source      haloToCell;
+        steps       3;          // Number of grow/shrink steps to use
+    }
 }
 
 
diff --git a/etc/caseDicts/postProcessing/fields/log b/etc/caseDicts/postProcessing/fields/log
new file mode 100644
index 0000000000000000000000000000000000000000..ceb852f5473368c7206ef513edc25cd8367f9f86
--- /dev/null
+++ b/etc/caseDicts/postProcessing/fields/log
@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Version:  v1912
+    \\  /    A nd           | Website:  www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+Description
+    Calculates the natural logarithm of an input volScalarField
+\*---------------------------------------------------------------------------*/
+
+type            log;
+libs            (fieldFunctionObjects);
+
+field           <fieldName>;
+
+executeControl  writeTime;
+writeControl    writeTime;
+scale           <scalar>;
+translate       <scalar>;
+
+// ************************************************************************* //
diff --git a/src/TurbulenceModels/turbulenceModels/RAS/derivedFvPatchFields/turbulentMixingLengthDissipationRateInlet/turbulentMixingLengthDissipationRateInletFvPatchScalarField.C b/src/TurbulenceModels/turbulenceModels/RAS/derivedFvPatchFields/turbulentMixingLengthDissipationRateInlet/turbulentMixingLengthDissipationRateInletFvPatchScalarField.C
index 01d8846f1e7c3ea1fedd54834899439c55d694c2..9f4f0a2a16b92cb56a65d991c21dcc086812a7b0 100644
--- a/src/TurbulenceModels/turbulenceModels/RAS/derivedFvPatchFields/turbulentMixingLengthDissipationRateInlet/turbulentMixingLengthDissipationRateInletFvPatchScalarField.C
+++ b/src/TurbulenceModels/turbulenceModels/RAS/derivedFvPatchFields/turbulentMixingLengthDissipationRateInlet/turbulentMixingLengthDissipationRateInletFvPatchScalarField.C
@@ -48,8 +48,9 @@ turbulentMixingLengthDissipationRateInletFvPatchScalarField
 )
 :
     inletOutletFvPatchScalarField(p, iF),
+    kName_("k"),
     mixingLength_(0.0),
-    kName_("k")
+    Cmu_(0.0)
 {
     this->refValue() = 0.0;
     this->refGrad() = 0.0;
@@ -67,8 +68,9 @@ turbulentMixingLengthDissipationRateInletFvPatchScalarField
 )
 :
     inletOutletFvPatchScalarField(ptf, p, iF, mapper),
+    kName_(ptf.kName_),
     mixingLength_(ptf.mixingLength_),
-    kName_(ptf.kName_)
+    Cmu_(ptf.Cmu_)
 {}
 
 
@@ -81,8 +83,12 @@ turbulentMixingLengthDissipationRateInletFvPatchScalarField
 )
 :
     inletOutletFvPatchScalarField(p, iF),
-    mixingLength_(dict.get<scalar>("mixingLength")),
-    kName_(dict.getOrDefault<word>("k", "k"))
+    kName_(dict.getOrDefault<word>("k", "k")),
+    mixingLength_
+    (
+        dict.getCheck<scalar>("mixingLength", scalarMinMax::ge(SMALL))
+    ),
+    Cmu_(dict.getCheckOrDefault<scalar>("Cmu", 0.09, scalarMinMax::ge(SMALL)))
 {
     this->phiName_ = dict.getOrDefault<word>("phi", "phi");
 
@@ -101,8 +107,9 @@ turbulentMixingLengthDissipationRateInletFvPatchScalarField
 )
 :
     inletOutletFvPatchScalarField(ptf),
+    kName_(ptf.kName_),
     mixingLength_(ptf.mixingLength_),
-    kName_(ptf.kName_)
+    Cmu_(ptf.Cmu_)
 {}
 
 
@@ -114,8 +121,9 @@ turbulentMixingLengthDissipationRateInletFvPatchScalarField
 )
 :
     inletOutletFvPatchScalarField(ptf, iF),
+    kName_(ptf.kName_),
     mixingLength_(ptf.mixingLength_),
-    kName_(ptf.kName_)
+    Cmu_(ptf.Cmu_)
 {}
 
 
@@ -138,10 +146,9 @@ void turbulentMixingLengthDissipationRateInletFvPatchScalarField::updateCoeffs()
         )
     );
 
-    const scalar Cmu =
-        turbModel.coeffDict().getOrDefault<scalar>("Cmu", 0.09);
+    Cmu_ = turbModel.coeffDict().getOrDefault<scalar>("Cmu", Cmu_);
 
-    const scalar Cmu75 = pow(Cmu, 0.75);
+    const scalar Cmu75 = pow(Cmu_, 0.75);
 
     const fvPatchScalarField& kp =
         patch().lookupPatchField<volScalarField, scalar>(kName_);
@@ -149,7 +156,7 @@ void turbulentMixingLengthDissipationRateInletFvPatchScalarField::updateCoeffs()
     const fvsPatchScalarField& phip =
         patch().lookupPatchField<surfaceScalarField, scalar>(this->phiName_);
 
-    this->refValue() = Cmu75*kp*sqrt(kp)/mixingLength_;
+    this->refValue() = (Cmu75/mixingLength_)*pow(kp, 1.5);
     this->valueFraction() = 1.0 - pos0(phip);
 
     inletOutletFvPatchScalarField::updateCoeffs();
diff --git a/src/TurbulenceModels/turbulenceModels/RAS/derivedFvPatchFields/turbulentMixingLengthDissipationRateInlet/turbulentMixingLengthDissipationRateInletFvPatchScalarField.H b/src/TurbulenceModels/turbulenceModels/RAS/derivedFvPatchFields/turbulentMixingLengthDissipationRateInlet/turbulentMixingLengthDissipationRateInletFvPatchScalarField.H
index bab77523126cee5af6636136c24d8417e1118a38..9246a2853632b85660aa888ec03d1fc64333dc3b 100644
--- a/src/TurbulenceModels/turbulenceModels/RAS/derivedFvPatchFields/turbulentMixingLengthDissipationRateInlet/turbulentMixingLengthDissipationRateInletFvPatchScalarField.H
+++ b/src/TurbulenceModels/turbulenceModels/RAS/derivedFvPatchFields/turbulentMixingLengthDissipationRateInlet/turbulentMixingLengthDissipationRateInletFvPatchScalarField.H
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -30,43 +31,63 @@ Group
     grpRASBoundaryConditions grpInletBoundaryConditions
 
 Description
-    This boundary condition provides a turbulence dissipation, \f$\epsilon\f$
-    (epsilon) inlet condition based on a specified mixing length.  The patch
-    values are calculated using:
+    This boundary condition provides an inlet condition for turbulent kinetic
+    energy dissipation rate, i.e. \c epsilon, based on a specified mixing
+    length.  The patch values are calculated using:
 
         \f[
             \epsilon_p = \frac{C_{\mu}^{0.75} k^{1.5}}{L}
         \f]
 
     where
-
     \vartable
-        \epsilon_p | patch epsilon values
-        C_{\mu} | Model coefficient, set to 0.09
-        k       | turbulence kinetic energy
-        L       | length scale
+      \epsilon_p | Patch epsilon values     [m2/s3]
+      C_\mu      | Empirical model constant retrived from turbulence model
+      k          | Turbulent kinetic energy [m2/s2]
+      L          | Mixing length scale      [m]
     \endvartable
 
 Usage
-    \table
-        Property     | Description             | Required    | Default value
-        mixingLength | Length scale [m]        | yes         |
-        phi          | flux field name         | no          | phi
-        k            | turbulence kinetic energy field name | no | k
-    \endtable
-
     Example of the boundary condition specification:
     \verbatim
     <patchName>
     {
+        // Mandatory entries (unmodifiable)
         type            turbulentMixingLengthDissipationRateInlet;
+
+        // Mandatory entries (runtime modifiable)
         mixingLength    0.005;
-        value           uniform 200;   // placeholder
+
+        // Optional entries (runtime modifiable)
+        Cmu             0.09;
+        k               k;
+        phi             phi;
+
+        // Placeholder
+        value           uniform 200;
     }
     \endverbatim
 
+    where the entries mean:
+    \table
+      Property     | Description                  | Type   | Req'd | Dflt
+      mixingLength | Mixing length scale [m]      | scalar |  yes  | -
+      Cmu          | Empirical model constant     | scalar |  no   | 0.09
+      phi          | Name of flux field           | word   |  no   | phi
+      k            | Name of turbulent kinetic energy field | word | no | k
+    \endtable
+
 Note
-    In the event of reverse flow, a zero-gradient condition is applied
+    - The boundary condition is derived from \c inletOutlet condition.
+      Therefore, in the event of reverse flow, a zero-gradient condition
+      is applied.
+    - The order of precedence to input the empirical model constant \c Cmu is:
+      turbulence model, boundary condition dictionary, and default value=0.09.
+    - The empirical model constant \c Cmu is not a spatiotemporal variant field.
+      Therefore, the use of the boundary condition may not be fully consistent
+      with the turbulence models where \c Cmu is a variant field, such as
+      \c realizableKE closure model in this respect. Nevertheless, workflow
+      observations suggest that the matter poses no importance.
 
 See also
     Foam::inletOutletFvPatchField
@@ -94,14 +115,17 @@ class turbulentMixingLengthDissipationRateInletFvPatchScalarField
 :
     public inletOutletFvPatchScalarField
 {
-    // Private data
-
-        //- Turbulent length scale
-        scalar mixingLength_;
+    // Private Data
 
         //- Name of the turbulent kinetic energy field
         word kName_;
 
+        //- Mixing length scale
+        scalar mixingLength_;
+
+        //- Empirical model constant
+        scalar Cmu_;
+
 
 public:
 
@@ -127,8 +151,8 @@ public:
         );
 
         //- Construct by mapping given
-        //  turbulentMixingLengthDissipationRateInletFvPatchScalarField
-        //  onto a new patch
+        //- turbulentMixingLengthDissipationRateInletFvPatchScalarField
+        //- onto a new patch
         turbulentMixingLengthDissipationRateInletFvPatchScalarField
         (
             const turbulentMixingLengthDissipationRateInletFvPatchScalarField&,
@@ -179,7 +203,7 @@ public:
         }
 
 
-    // Member functions
+    // Member Functions
 
         //- Update the coefficients associated with the patch field
         virtual void updateCoeffs();
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/turbulentDigitalFilterInletFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/turbulentDigitalFilterInletFvPatchVectorField.C
index bcbf9caa899e35a8e53c6fcc57f134d51a9be298..2b95c58c82e3a29ff234dc2c61dcde36cc72293f 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/turbulentDigitalFilterInletFvPatchVectorField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/turbulentDigitalFilterInletFvPatchVectorField.C
@@ -30,21 +30,6 @@ License
 #include "mathematicalConstants.H"
 #include "addToRunTimeSelectionTable.H"
 
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-const Foam::Enum
-<
-    Foam::turbulentDigitalFilterInletFvPatchVectorField::variantType
->
-Foam::turbulentDigitalFilterInletFvPatchVectorField::variantNames
-({
-    { variantType::DIGITAL_FILTER, "digitalFilter" },
-    { variantType::DIGITAL_FILTER, "DFM" },
-    { variantType::FORWARD_STEPWISE, "forwardStepwise" },
-    { variantType::FORWARD_STEPWISE, "reducedDigitalFilter" },
-    { variantType::FORWARD_STEPWISE, "FSM" }
-});
-
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 const Foam::pointToPointPlanarInterpolation&
@@ -114,10 +99,10 @@ Foam::turbulentDigitalFilterInletFvPatchVectorField::patchMapper() const
 
 
 Foam::List<Foam::Pair<Foam::label>>
-Foam::turbulentDigitalFilterInletFvPatchVectorField::patchIndexPairs()
+Foam::turbulentDigitalFilterInletFvPatchVectorField::indexPairs()
 {
     // Get patch normal direction into the domain
-    const vector nf(computePatchNormal());
+    const vector nf((-gAverage(patch().nf())).normalise());
 
     // Find the second local coordinate direction
     direction minCmpt = 0;
@@ -163,7 +148,7 @@ Foam::turbulentDigitalFilterInletFvPatchVectorField::patchIndexPairs()
     const boundBox bb(localPos);
 
     // Normalise to (i, j) coordinates
-    const Vector2D<label> n(planeDivisions_.first(), planeDivisions_.second());
+    const Vector2D<label> n(n_.first(), n_.second());
     invDelta_[0] = n[0]/bb.span()[0];
     invDelta_[1] = n[1]/bb.span()[1];
     const point localMinPt(bb.min());
@@ -187,8 +172,7 @@ Foam::turbulentDigitalFilterInletFvPatchVectorField::patchIndexPairs()
 }
 
 
-void Foam::turbulentDigitalFilterInletFvPatchVectorField::
-checkRTensorRealisable() const
+void Foam::turbulentDigitalFilterInletFvPatchVectorField::checkR() const
 {
     const vectorField& faceCentres = this->patch().patch().faceCentres();
 
@@ -198,65 +182,65 @@ checkRTensorRealisable() const
         {
             FatalErrorInFunction
                 << "Reynolds stress tensor component Rxx cannot be negative"
-                   "or zero, where Rxx = " << R_[facei].xx() << " at the face "
-                   "centre = " << faceCentres[facei] << exit(FatalError);
+                << " or zero, where Rxx = " << R_[facei].xx()
+                << " at the face centre = " << faceCentres[facei]
+                << exit(FatalError);
         }
 
         if (R_[facei].yy() < 0 || R_[facei].zz() < 0)
         {
             FatalErrorInFunction
                 << "Reynolds stress tensor components Ryy or Rzz cannot be"
-                << "negative where Ryy = " << R_[facei].yy() << ", and Rzz = "
-                << R_[facei].zz() << " at the face centre = "
-                << faceCentres[facei] << exit(FatalError);
+                << " negative where Ryy = " << R_[facei].yy()
+                << ", and Rzz = " << R_[facei].zz()
+                << " at the face centre = " << faceCentres[facei]
+                << exit(FatalError);
         }
 
-        scalar term0 = R_[facei].xx()*R_[facei].yy() - sqr(R_[facei].xy());
+        const scalar x0 = R_[facei].xx()*R_[facei].yy() - sqr(R_[facei].xy());
 
-        if (term0 <= 0)
+        if (x0 <= 0)
         {
             FatalErrorInFunction
                 << "Reynolds stress tensor component group, Rxx*Ryy - Rxy^2"
-                << "cannot be negative or zero at the face centre = "
-                << faceCentres[facei] << exit(FatalError);
+                << " cannot be negative or zero"
+                << " at the face centre = " << faceCentres[facei]
+                << exit(FatalError);
         }
 
-        scalar term1 = R_[facei].zz() - sqr(R_[facei].xz())/R_[facei].xx();
-        scalar term2 =
-            sqr(R_[facei].yz() - R_[facei].xy()*R_[facei].xz()
-           /(R_[facei].xx()*term0));
-        scalar term3 = term1 - term2;
+        const scalar x1 = R_[facei].zz() - sqr(R_[facei].xz())/R_[facei].xx();
+        const scalar x2 = sqr(R_[facei].yz() - R_[facei].xy()*R_[facei].xz()
+                         /(R_[facei].xx()*x0));
+        const scalar x3 = x1 - x2;
 
-        if (term3 < 0)
+        if (x3 < 0)
         {
             FatalErrorInFunction
-                << "Reynolds stress tensor component group,"
+                << "Reynolds stress tensor component group, "
                 << "Rzz - Rxz^2/Rxx - (Ryz - Rxy*Rxz/(Rxx*(Rxx*Ryy - Rxy^2)))^2"
-                << "cannot be negative at the face centre = "
-                << faceCentres[facei] << exit(FatalError);
+                << " cannot be negative at the face centre = "
+                << faceCentres[facei]
+                << exit(FatalError);
         }
     }
 
-    #ifdef FULLDEBUG
-    Info<< "Ends: checkRTensorRealisable()."
-        << " Reynolds tensor (on patch) is consistent." << nl;
-    #endif
+    Info<< "  # Reynolds stress tensor on patch is consistent #" << endl;
 }
 
 
 Foam::symmTensorField Foam::turbulentDigitalFilterInletFvPatchVectorField::
-computeLundWuSquires() const
+calcLund() const
 {
-    checkRTensorRealisable();
+    checkR();
 
-    symmTensorField LundWuSquires(symmTensorField(R_.size()));
+    symmTensorField Lund(symmTensorField(R_.size()));
 
     forAll(R_, facei)
     {
         const symmTensor& R = R_[facei];
-        symmTensor& lws = LundWuSquires[facei];
+        symmTensor& lws = Lund[facei];
 
-        // (Klein et al., 2003, Eq. 5)
+        // (KSJ:Eq. 5)
         lws.xx() = Foam::sqrt(R.xx());
         lws.xy() = R.xy()/lws.xx();
         lws.xz() = R.xz()/lws.xx();
@@ -265,128 +249,96 @@ computeLundWuSquires() const
         lws.zz() = Foam::sqrt(R.zz() - sqr(lws.xz()) - sqr(lws.yz()));
     }
 
-    #ifdef FULLDEBUG
-    Info<< "Ends: computeLundWuSquires()." << nl;
-    #endif
-
-    return LundWuSquires;
-}
-
-
-Foam::vector Foam::turbulentDigitalFilterInletFvPatchVectorField::
-computePatchNormal() const
-{
-    vector patchNormal(-gAverage(patch().nf()));
-    return patchNormal.normalise();
+    return Lund;
 }
 
 
 Foam::scalar Foam::turbulentDigitalFilterInletFvPatchVectorField::
-computeInitialFlowRate() const
+calcFlowRate() const
 {
-    const vector patchNormal(computePatchNormal());
+    const vector patchNormal((-gAverage(patch().nf())).normalise());
     return gSum((UMean_ & patchNormal)*patch().magSf());
 }
 
 
-void Foam::turbulentDigitalFilterInletFvPatchVectorField::convertToTimeScale
-(
-    tensor& L
-) const
-{
-    if (isTaylorHypot_)
-    {
-        forAll(L.x(), i)
-        {
-            L[i] /= patchNormalSpeed_;
-        }
-
-    #ifdef FULLDEBUG
-    Info<< "Ends: convertToTimeScale()."
-        << "Streamwise integral length scales are converted to time scales via"
-        << "Taylor's frozen turbulence hypothesis" << nl;
-    #endif
-    }
-}
-
-
 Foam::tensor Foam::turbulentDigitalFilterInletFvPatchVectorField::
-convertScalesToGridUnits
+meterToCell
 (
     const tensor& L
 ) const
 {
     tensor Ls(L);
-    convertToTimeScale(Ls);
+
+    // Convert streamwise integral length scales to integral
+    // time scales by using Taylor's frozen turbulence hypothesis
+    forAll(Ls.x(), i)
+    {
+        Ls[i] /= Ubulk_;
+    }
+
     const scalar invDeltaT = 1.0/db().time().deltaTValue();
 
+    //  (KSJ:Eq. 13)
     Ls.row(0, Ls.x()*invDeltaT);
     Ls.row(1, Ls.y()*invDelta_[0]);
     Ls.row(2, Ls.z()*invDelta_[1]);
 
-    #ifdef FULLDEBUG
-    Info<< "Ends: convertScalesToGridUnits()."
-        << "Units of input length scales are converted from metres to"
-        << "virtual-patch cell size/time-step" << nl;
-    #endif
-
     return Ls;
 }
 
 
 Foam::List<Foam::label> Foam::turbulentDigitalFilterInletFvPatchVectorField::
-initLenRandomBox() const
+initBox() const
 {
     label initValue = 0;
     label rangeModifier = 0;
 
-    if (variant_ == variantType::FORWARD_STEPWISE)
+    if (fsm_)
     {
-        // Initialise with 1 since x-dir possess 1 node with this variant
-        initValue = pTraits<label>::nComponents;
+        // Initialise with 1 since streamwise-dir possesses 1 cell in FSM
+        initValue = 1;
         rangeModifier = pTraits<vector>::nComponents;
     }
 
-    List<label> lenRandomBox(pTraits<tensor>::nComponents, initValue);
-    Vector<label> lenGrid
+    List<label> szBox(pTraits<tensor>::nComponents, initValue);
+    Vector<label> szPlane
     (
-        pTraits<label>::nComponents,
-        planeDivisions_.first(),
-        planeDivisions_.second()
+        1,
+        n_.first(),
+        n_.second()
     );
 
-    // Initialise: Main convenience factor, lenRandomBox_
     for
     (
-        const label& i
+        const label i
       : labelRange(rangeModifier, pTraits<tensor>::nComponents - rangeModifier)
     )
     {
         // Slicing index
-        const label sliceI = label(i/pTraits<vector>::nComponents);
+        const label slicei = label(i/pTraits<vector>::nComponents);
 
-        // Refer to 'computeFilterCoeffs()'
+        // Refer to "calcFilterCoeffs()"
         const label n = ceil(L_[i]);
         const label twiceN = 4*n;
 
-        // Initialise: Random-number set sizes
-        lenRandomBox[i] = lenGrid[sliceI] + twiceN;
+        // Size of random-number sets
+        szBox[i] = szPlane[slicei] + twiceN;
     }
 
-    return lenRandomBox;
+    return szBox;
 }
 
 
 Foam::List<Foam::label> Foam::turbulentDigitalFilterInletFvPatchVectorField::
-initBoxFactors2D() const
+initFactors2D() const
 {
     List<label> boxFactors2D(pTraits<vector>::nComponents);
 
     for (direction dir = 0; dir < pTraits<vector>::nComponents; ++dir)
     {
         boxFactors2D[dir] =
-            lenRandomBox_[pTraits<vector>::nComponents + dir]
-           *lenRandomBox_[pTraits<symmTensor>::nComponents + dir];
+            szBox_[pTraits<vector>::nComponents + dir]
+           *szBox_[pTraits<symmTensor>::nComponents + dir];
     }
 
     return boxFactors2D;
@@ -394,13 +346,13 @@ initBoxFactors2D() const
 
 
 Foam::List<Foam::label> Foam::turbulentDigitalFilterInletFvPatchVectorField::
-initBoxFactors3D() const
+initFactors3D() const
 {
     List<label> boxFactors3D(pTraits<vector>::nComponents);
 
     for (direction dir = 0; dir < pTraits<vector>::nComponents; ++dir)
     {
-        boxFactors3D[dir] = randomBoxFactors2D_[dir]*lenRandomBox_[dir];
+        boxFactors3D[dir] = boxFactors2D_[dir]*szBox_[dir];
     }
 
     return boxFactors3D;
@@ -408,14 +360,13 @@ initBoxFactors3D() const
 
 
 Foam::List<Foam::label> Foam::turbulentDigitalFilterInletFvPatchVectorField::
-initBoxPlaneFactors() const
+initPlaneFactors() const
 {
     List<label> planeFactors(pTraits<vector>::nComponents);
 
     for (direction dir = 0; dir < pTraits<vector>::nComponents; ++dir)
     {
-        planeFactors[dir] =
-            randomBoxFactors2D_[dir]*(lenRandomBox_[dir] - pTraits<label>::one);
+        planeFactors[dir] = boxFactors2D_[dir]*(szBox_[dir] - 1);
     }
 
     return planeFactors;
@@ -423,29 +374,23 @@ initBoxPlaneFactors() const
 
 
 Foam::List<Foam::List<Foam::scalar>>
-Foam::turbulentDigitalFilterInletFvPatchVectorField::fillRandomBox()
+Foam::turbulentDigitalFilterInletFvPatchVectorField::fillBox()
 {
-    List<List<scalar>> randomBox(pTraits<vector>::nComponents, List<scalar>());
+    List<List<scalar>> box(pTraits<vector>::nComponents, List<scalar>());
 
     // Initialise: Remaining convenience factors for (e1 e2 e3)
     for (direction dir = 0; dir < pTraits<vector>::nComponents; ++dir)
     {
         // Initialise: Random-box content with random-number sets
-        randomBox[dir] =
-            generateRandomSet<List<scalar>, scalar>(randomBoxFactors3D_[dir]);
+        box[dir] = randomSet<List<scalar>, scalar>(boxFactors3D_[dir]);
     }
 
-    #ifdef FULLDEBUG
-    Info<< "Ends: fillRandomBox()."
-        << "Random-number box is created." << nl;
-    #endif
-
-    return randomBox;
+    return box;
 }
 
 
 Foam::List<Foam::List<Foam::scalar>>
-Foam::turbulentDigitalFilterInletFvPatchVectorField::computeFilterCoeffs() const
+Foam::turbulentDigitalFilterInletFvPatchVectorField::calcFilterCoeffs() const
 {
     List<List<scalar>> filterCoeffs
     (
@@ -455,7 +400,7 @@ Foam::turbulentDigitalFilterInletFvPatchVectorField::computeFilterCoeffs() const
 
     label initValue = 0;
 
-    if(variant_ == variantType::FORWARD_STEPWISE)
+    if (fsm_)
     {
         initValue = pTraits<vector>::nComponents;
     }
@@ -463,11 +408,11 @@ Foam::turbulentDigitalFilterInletFvPatchVectorField::computeFilterCoeffs() const
     for (direction dir = initValue; dir < pTraits<tensor>::nComponents; ++dir)
     {
         // The smallest filter width larger than length scale
-        // (Klein et al., 2003, 'n' in Eq. 15)
+        // (KSJ:'n' in Eq. 15)
         const label n  = ceil(L_[dir]);
 
         // Full filter-width (except mid-zero) according to the condition
-        // (Klein et al., 2003, Eq. 15 whereat N is minimum =2n)
+        // (KSJ:Eq. 15 whereat N is minimum =2n)
         const label twiceN = 4*n;
 
         // Initialise filter-coeffs containers with full filter-width size
@@ -477,8 +422,7 @@ Foam::turbulentDigitalFilterInletFvPatchVectorField::computeFilterCoeffs() const
         // First element: -N within [-N, N]
         const scalar initElem = -2*scalar(n);
 
-        // Container initialised with [-N, N]
-        // (Klein et al., 2003, p. 658, Item-b)
+        // Container initialised with [-N, N] (KSJ:p. 658, item-b)
         std::iota
         (
             filterCoeffs[dir].begin(),
@@ -486,50 +430,43 @@ Foam::turbulentDigitalFilterInletFvPatchVectorField::computeFilterCoeffs() const
             initElem
         );
 
-        // Compute filter-coeffs
-        // (Klein et al., 2003, Eq. 14 (Gaussian))
-        // (Bercin et al., 2018, Fig. 9 (Exponential))
+        // Compute filter-coeffs (KSJ:Eq. 14 (Gaussian))
         List<scalar> fTemp(filterCoeffs[dir]);
         scalar fSum = 0;
         const scalar nSqr = n*n;
 
-        if (isGaussian_)
+        if (Gaussian_)
         {
-            fTemp = sqr(exp(modelConst_*sqr(fTemp)/nSqr));
+            fTemp = sqr(exp(C1_*sqr(fTemp)/nSqr));
             fSum = Foam::sqrt(sum(fTemp));
 
             filterCoeffs[dir] =
-                exp(modelConst_*sqr(filterCoeffs[dir])/nSqr)/fSum;
+                exp(C1_*sqr(filterCoeffs[dir])/nSqr)/fSum;
         }
         else
         {
-            fTemp = sqr(exp(modelConst_*mag(fTemp)/n));
+            fTemp = sqr(exp(C1_*mag(fTemp)/n));
             fSum = Foam::sqrt(sum(fTemp));
 
-            filterCoeffs[dir] = exp(modelConst_*mag(filterCoeffs[dir])/n)/fSum;
+            filterCoeffs[dir] = exp(C1_*mag(filterCoeffs[dir])/n)/fSum;
         }
     }
 
-    #ifdef FULLDEBUG
-    Info<< "Ends: computeFilterCoeffs()."
-        << " Filter coefficients are computed." << nl;
-    #endif
-
     return filterCoeffs;
 }
 
 
-void Foam::turbulentDigitalFilterInletFvPatchVectorField::rndShiftRefill()
+void Foam::turbulentDigitalFilterInletFvPatchVectorField::shiftRefill()
 {
-    forAll(randomBox_, dir)
+    forAll(box_, dir)
     {
-        List<scalar>& r = randomBox_[dir];
+        List<scalar>& r = box_[dir];
 
         // Shift forward from the back to the front / First Out
-        inplaceRotateList(r, randomBoxFactors2D_[dir]);
+        inplaceRotateList(r, boxFactors2D_[dir]);
 
         // Refill the back with a new random-number set / First In
-        for (label i = 0; i < randomBoxFactors2D_[dir]; ++i)
+        for (label i = 0; i < boxFactors2D_[dir]; ++i)
         {
             r[i] = rndGen_.GaussNormal<scalar>();
         }
@@ -537,7 +474,7 @@ void Foam::turbulentDigitalFilterInletFvPatchVectorField::rndShiftRefill()
 }
 
 
-void Foam::turbulentDigitalFilterInletFvPatchVectorField::mapFilteredRandomBox
+void Foam::turbulentDigitalFilterInletFvPatchVectorField::mapFilteredBox
 (
     vectorField& U
 )
@@ -546,24 +483,24 @@ void Foam::turbulentDigitalFilterInletFvPatchVectorField::mapFilteredRandomBox
     {
         const label xf = x.first();
         const label xs = x.second();
-        U[xf][0] = filteredRandomBox_[0][xs];
-        U[xf][1] = filteredRandomBox_[1][xs];
-        U[xf][2] = filteredRandomBox_[2][xs];
+        U[xf][0] = filteredBox_[0][xs];
+        U[xf][1] = filteredBox_[1][xs];
+        U[xf][2] = filteredBox_[2][xs];
     }
 }
 
 
-void Foam::turbulentDigitalFilterInletFvPatchVectorField::embedOnePointCorrs
+void Foam::turbulentDigitalFilterInletFvPatchVectorField::onePointCorrs
 (
     vectorField& U
 ) const
 {
-    forAll(LundWuSquires_, facei)
+    forAll(Lund_, facei)
     {
         vector& Us = U[facei];
-        const symmTensor& lws = LundWuSquires_[facei];
+        const symmTensor& lws = Lund_[facei];
 
-        // (Klein et al. p. 658, Item-e)
+        // (KSJ:p. 658, item-e)
         Us.z() = Us.x()*lws.xz() + Us.y()*lws.yz() + Us.z()*lws.zz();
         Us.y() = Us.x()*lws.xy() + Us.y()*lws.yy();
         Us.x() = Us.x()*lws.xx();
@@ -571,42 +508,24 @@ void Foam::turbulentDigitalFilterInletFvPatchVectorField::embedOnePointCorrs
 }
 
 
-void Foam::turbulentDigitalFilterInletFvPatchVectorField::embedMeanVelocity
-(
-    vectorField& U
-) const
-{
-    U += UMean_;
-}
-
-
-void Foam::turbulentDigitalFilterInletFvPatchVectorField::correctFlowRate
-(
-    vectorField& U
-) const
-{
-    U *= (initialFlowRate_/gSum(U & -patch().Sf()));
-}
-
-
-void Foam::turbulentDigitalFilterInletFvPatchVectorField::embedTwoPointCorrs()
+void Foam::turbulentDigitalFilterInletFvPatchVectorField::twoPointCorrs()
 {
     for (direction dir = 0; dir < pTraits<vector>::nComponents; ++dir)
     {
-        List<scalar>& in = randomBox_[dir];
-        List<scalar>& out = filteredRandomBox_[dir];
+        List<scalar>& in = box_[dir];
+        List<scalar>& out = filteredBox_[dir];
         const List<scalar>& filter1 = filterCoeffs_[dir];
         const List<scalar>& filter2 = filterCoeffs_[3 + dir];
         const List<scalar>& filter3 = filterCoeffs_[6 + dir];
 
-        const label sz1 = lenRandomBox_[dir];
-        const label sz2 = lenRandomBox_[3 + dir];
-        const label sz3 = lenRandomBox_[6 + dir];
+        const label sz1 = szBox_[dir];
+        const label sz2 = szBox_[3 + dir];
+        const label sz3 = szBox_[6 + dir];
         const label szfilter1 = filterCoeffs_[dir].size();
         const label szfilter2 = filterCoeffs_[3 + dir].size();
         const label szfilter3 = filterCoeffs_[6 + dir].size();
-        const label sz23 = randomBoxFactors2D_[dir];
-        const label sz123 = randomBoxFactors3D_[dir];
+        const label sz23 = boxFactors2D_[dir];
+        const label sz123 = boxFactors3D_[dir];
         const label validSlice2 = sz2 - (szfilter2 - label(1));
         const label validSlice3 = sz3 - (szfilter3 - label(1));
 
@@ -705,116 +624,31 @@ void Foam::turbulentDigitalFilterInletFvPatchVectorField::embedTwoPointCorrs()
 }
 
 
-void Foam::turbulentDigitalFilterInletFvPatchVectorField::computeDFM
-(
-    vectorField& U
-)
-{
-    #ifdef FULLDEBUG
-    Info<< "Starts: computeDFM()" << nl;
-    #endif
-
-    if (Pstream::master())
-    {
-        embedTwoPointCorrs();
-        rndShiftRefill();
-    }
-
-    Pstream::scatter(filteredRandomBox_);
-
-    mapFilteredRandomBox(U);
-
-    embedOnePointCorrs(U);
-
-    embedMeanVelocity(U);
-
-    if (isCorrectedFlowRate_)
-    {
-        correctFlowRate(U);
-    }
-
-    #ifdef FULLDEBUG
-    Info<< "Ends: computeDFM()" << nl;
-    #endif
-}
-
-
-void Foam::turbulentDigitalFilterInletFvPatchVectorField::computeReducedDFM
-(
-    vectorField& U
-)
-{
-    #ifdef FULLDEBUG
-    Info<< "Starts: computeReducedDFM()" << nl;
-    #endif
-
-    if (Pstream::master())
-    {
-        embedTwoPointCorrs();
-        rndShiftRefill();
-    }
-
-    Pstream::scatter(filteredRandomBox_);
-
-    mapFilteredRandomBox(U);
-
-    computeFSM(U);
-
-    embedOnePointCorrs(U);
-
-    embedMeanVelocity(U);
-
-    if (isCorrectedFlowRate_)
-    {
-        correctFlowRate(U);
-    }
-
-    #ifdef FULLDEBUG
-    Info<< "Ends: computeReducedDFM()" << nl;
-    #endif
-}
-
-
 Foam::List<Foam::scalar> Foam::turbulentDigitalFilterInletFvPatchVectorField::
-computeConstList1FSM() const
+calcCoeffs1FSM() const
 {
-    List<scalar> constList1FSM(pTraits<vector>::nComponents);
+    List<scalar> coeffs1FSM(pTraits<vector>::nComponents);
 
     forAll(L_.x(), i)
     {
-        constList1FSM[i] = exp(const1FSM_/L_[i]);
+        coeffs1FSM[i] = exp(C1FSM_/L_[i]);
     }
 
-    return constList1FSM;
+    return coeffs1FSM;
 }
 
 
 Foam::List<Foam::scalar> Foam::turbulentDigitalFilterInletFvPatchVectorField::
-computeConstList2FSM() const
+calcCoeffs2FSM() const
 {
-    List<scalar> constList2FSM(pTraits<vector>::nComponents);
+    List<scalar> coeffs2FSM(pTraits<vector>::nComponents);
 
     forAll(L_.x(), i)
     {
-        constList2FSM[i] = Foam::sqrt(1.0 - exp(const2FSM_/L_[i]));
+        coeffs2FSM[i] = Foam::sqrt(1.0 - exp(C2FSM_/L_[i]));
     }
 
-    return constList2FSM;
-}
-
-
-void Foam::turbulentDigitalFilterInletFvPatchVectorField::computeFSM
-(
-    vectorField& U
-)
-{
-    for (direction dir = 0; dir < pTraits<vector>::nComponents; ++dir)
-    {
-        U.component(dir) = U0_.component(dir)*constList1FSM_[dir]
-                          +U.component(dir)*constList2FSM_[dir];
-    }
-
-    U0_ = U;
+    return coeffs2FSM;
 }
 
 
@@ -829,44 +663,40 @@ turbulentDigitalFilterInletFvPatchVectorField
 :
     fixedValueFvPatchField<vector>(p, iF),
     mapperPtr_(nullptr),
-    variant_(variantType::DIGITAL_FILTER),
-    isGaussian_(true),
-    isFixedSeed_(true),
-    isContinuous_(false),
-    isCorrectedFlowRate_(true),
-    interpolateR_(false),
-    interpolateUMean_(false),
-    isInsideMesh_(false),
-    isTaylorHypot_(true),
+    fsm_(false),
+    Gaussian_(true),
+    fixSeed_(true),
+    continuous_(false),
+    correctFlowRate_(true),
+    interpR_(false),
+    interpUMean_(false),
     mapMethod_("nearestCell"),
     curTimeIndex_(-1),
-    tiny_(1e-8),
-    patchNormalSpeed_(Zero),
-    modelConst_(-0.5*constant::mathematical::pi),
+    Ubulk_(0.0),
+    C1_(-0.5*constant::mathematical::pi),
     perturb_(1e-5),
-    initialFlowRate_(pTraits<scalar>::one),
+    flowRate_(1.0),
     rndGen_(1234),
-    planeDivisions_(Zero, Zero),
+    n_(0, 0),
     invDelta_(Zero),
     indexPairs_(Zero),
     R_(Zero),
-    LundWuSquires_(Zero),
+    Lund_(Zero),
     UMean_(Zero),
     Lbak_(Zero),
     L_(Zero),
-    const1FSM_(Zero),
-    const2FSM_(Zero),
-    constList1FSM_(Zero),
-    constList2FSM_(Zero),
-    lenRandomBox_(Zero),
-    randomBoxFactors2D_(Zero),
-    randomBoxFactors3D_(Zero),
+    C1FSM_(Zero),
+    C2FSM_(Zero),
+    coeffs1FSM_(Zero),
+    coeffs2FSM_(Zero),
+    szBox_(Zero),
+    boxFactors2D_(Zero),
+    boxFactors3D_(Zero),
     iNextToLastPlane_(Zero),
-    randomBox_(Zero),
+    box_(Zero),
     filterCoeffs_(Zero),
-    filteredRandomBox_(Zero),
-    U0_(Zero),
-    computeVariant(nullptr)
+    filteredBox_(Zero),
+    U0_(Zero)
 {}
 
 
@@ -881,44 +711,40 @@ turbulentDigitalFilterInletFvPatchVectorField
 :
     fixedValueFvPatchField<vector>(ptf, p, iF, mapper),
     mapperPtr_(nullptr),
-    variant_(ptf.variant_),
-    isGaussian_(ptf.isGaussian_),
-    isFixedSeed_(ptf.isFixedSeed_),
-    isContinuous_(ptf.isContinuous_),
-    isCorrectedFlowRate_(ptf.isCorrectedFlowRate_),
-    interpolateR_(ptf.interpolateR_),
-    interpolateUMean_(ptf.interpolateUMean_),
-    isInsideMesh_(ptf.isInsideMesh_),
-    isTaylorHypot_(ptf.isTaylorHypot_),
+    fsm_(ptf.fsm_),
+    Gaussian_(ptf.Gaussian_),
+    fixSeed_(ptf.fixSeed_),
+    continuous_(ptf.continuous_),
+    correctFlowRate_(ptf.correctFlowRate_),
+    interpR_(ptf.interpR_),
+    interpUMean_(ptf.interpUMean_),
     mapMethod_(ptf.mapMethod_),
     curTimeIndex_(-1),
-    tiny_(ptf.tiny_),
-    patchNormalSpeed_(ptf.patchNormalSpeed_),
-    modelConst_(ptf.modelConst_),
+    Ubulk_(ptf.Ubulk_),
+    C1_(ptf.C1_),
     perturb_(ptf.perturb_),
-    initialFlowRate_(ptf.initialFlowRate_),
+    flowRate_(ptf.flowRate_),
     rndGen_(ptf.rndGen_),
-    planeDivisions_(ptf.planeDivisions_),
+    n_(ptf.n_),
     invDelta_(ptf.invDelta_),
     indexPairs_(ptf.indexPairs_),
     R_(ptf.R_),
-    LundWuSquires_(ptf.LundWuSquires_),
+    Lund_(ptf.Lund_),
     UMean_(ptf.UMean_),
     Lbak_(ptf.Lbak_),
     L_(ptf.L_),
-    const1FSM_(ptf.const1FSM_),
-    const2FSM_(ptf.const2FSM_),
-    constList1FSM_(ptf.constList1FSM_),
-    constList2FSM_(ptf.constList2FSM_),
-    lenRandomBox_(ptf.lenRandomBox_),
-    randomBoxFactors2D_(ptf.randomBoxFactors2D_),
-    randomBoxFactors3D_(ptf.randomBoxFactors3D_),
+    C1FSM_(ptf.C1FSM_),
+    C2FSM_(ptf.C2FSM_),
+    coeffs1FSM_(ptf.coeffs1FSM_),
+    coeffs2FSM_(ptf.coeffs2FSM_),
+    szBox_(ptf.szBox_),
+    boxFactors2D_(ptf.boxFactors2D_),
+    boxFactors3D_(ptf.boxFactors3D_),
     iNextToLastPlane_(ptf.iNextToLastPlane_),
-    randomBox_(ptf.randomBox_),
+    box_(ptf.box_),
     filterCoeffs_(ptf.filterCoeffs_),
-    filteredRandomBox_(ptf.filteredRandomBox_),
-    U0_(ptf.U0_),
-    computeVariant(ptf.computeVariant)
+    filteredBox_(ptf.filteredBox_),
+    U0_(ptf.U0_)
 {}
 
 
@@ -932,163 +758,115 @@ turbulentDigitalFilterInletFvPatchVectorField
 :
     fixedValueFvPatchField<vector>(p, iF, dict),
     mapperPtr_(nullptr),
-    variant_
-    (
-        variantNames.getOrDefault
-        (
-            "variant",
-            dict,
-            variantType::DIGITAL_FILTER
-        )
-    ),
-    isGaussian_
-    (
-        (variant_ == variantType::FORWARD_STEPWISE)
-      ? false
-      : dict.getOrDefault("isGaussian", true)
-    ),
-    isFixedSeed_(dict.getOrDefault("isFixedSeed", true)),
-    isContinuous_(dict.getOrDefault("isContinuous", false)),
-    isCorrectedFlowRate_(dict.getOrDefault("isCorrectedFlowRate", true)),
-    interpolateR_(dict.getOrDefault("interpolateR", false)),
-    interpolateUMean_(dict.getOrDefault("interpolateUMean", false)),
-    isInsideMesh_(dict.getOrDefault("isInsideMesh", false)),
-    isTaylorHypot_(dict.getOrDefault("isTaylorHypot", true)),
+    fsm_(dict.getOrDefault("fsm", false)),
+    Gaussian_(fsm_ ? false : dict.getOrDefault("Gaussian", true)),
+    fixSeed_(dict.getOrDefault("fixSeed", true)),
+    continuous_(dict.getOrDefault("continuous", false)),
+    correctFlowRate_(dict.getOrDefault("correctFlowRate", true)),
+    interpR_(false),
+    interpUMean_(false),
     mapMethod_(dict.getOrDefault<word>("mapMethod", "nearestCell")),
     curTimeIndex_(-1),
-    tiny_(dict.getOrDefault<scalar>("threshold", 1e-8)),
-    patchNormalSpeed_
-    (
-        dict.getCheck<scalar>
-        (
-            "patchNormalSpeed",
-            scalarMinMax::ge(tiny_)
-        )
-    ),
-    modelConst_
+    Ubulk_(dict.getCheck<scalar>("Ubulk", scalarMinMax::ge(ROOTSMALL))),
+    C1_
     (
-        dict.getOrDefault<scalar>
-        (
-            "modelConst",
-           -0.5*constant::mathematical::pi
-        )
+        dict.getOrDefault<scalar>("C1", -0.5*constant::mathematical::pi)
     ),
-    perturb_(dict.getOrDefault<scalar>("perturb", 1e-5)),
-    initialFlowRate_(pTraits<scalar>::one),
-    rndGen_
+    perturb_
     (
-        isFixedSeed_
-      ? 1234
-      : time(0)
+        dict.getCheckOrDefault<scalar>("perturb", 1e-5, scalarMinMax::ge(SMALL))
     ),
-    planeDivisions_
+    flowRate_(1.0),
+    rndGen_(fixSeed_ ? 1234 : time(0)),
+    n_
     (
         dict.getCheck<Tuple2<label, label>>
         (
-            "planeDivisions",
-            [&](const Tuple2<label, label>& len)
+            "n",
+            [&](const Tuple2<label, label>& x)
             {
-                return tiny_ < min(len.first(), len.second()) ? true : false;
+                return min(x.first(), x.second()) > 0 ? true : false;
             }
         )
     ),
-    invDelta_(),
-    indexPairs_(patchIndexPairs()),
-    R_(interpolateOrRead<symmTensor>("R", dict, interpolateR_)),
-    LundWuSquires_(computeLundWuSquires()),
-    UMean_(interpolateOrRead<vector>("UMean", dict, interpolateUMean_)),
+    invDelta_(Zero),
+    indexPairs_(indexPairs()),
+    R_(interpolateOrRead<symmTensor>("R", dict, interpR_)),
+    Lund_(calcLund()),
+    UMean_(interpolateOrRead<vector>("UMean", dict, interpUMean_)),
     Lbak_
     (
         dict.getCheck<tensor>
         (
             "L",
-            [&](const tensor& l){return tiny_ < cmptMin(l) ? true : false;}
+            [&](const tensor& x){return cmptMin(x) > ROOTSMALL ? true : false;}
         )
     ),
-    L_(convertScalesToGridUnits(Lbak_)),
-    const1FSM_
+    L_(meterToCell(Lbak_)),
+    C1FSM_
     (
-        (variant_ == variantType::FORWARD_STEPWISE)
+        fsm_
       ? dict.getOrDefault<scalar>
         (
-            "const1FSM",
+            "C1FSM",
            -0.25*constant::mathematical::pi
         )
-      : Zero
+      : 0.0
     ),
-    const2FSM_
+    C2FSM_
     (
-        (variant_ == variantType::FORWARD_STEPWISE)
+        fsm_
       ? dict.getOrDefault<scalar>
         (
-            "const2FSM",
-           -0.5*constant::mathematical::pi
+            "C2FSM",
+            -0.5*constant::mathematical::pi
         )
-      : Zero
+      : 0.0
     ),
-    constList1FSM_
+    coeffs1FSM_(fsm_ ? calcCoeffs1FSM() : List<scalar>()),
+    coeffs2FSM_(fsm_ ? calcCoeffs2FSM() : List<scalar>()),
+    szBox_(initBox()),
+    boxFactors2D_(initFactors2D()),
+    boxFactors3D_(initFactors3D()),
+    iNextToLastPlane_(initPlaneFactors()),
+    box_
     (
-        (variant_ == variantType::FORWARD_STEPWISE)
-      ? computeConstList1FSM()
-      : List<scalar>()
-    ),
-    constList2FSM_
-    (
-        (variant_ == variantType::FORWARD_STEPWISE)
-      ? computeConstList2FSM()
-      : List<scalar>()
-    ),
-    lenRandomBox_(initLenRandomBox()),
-    randomBoxFactors2D_(initBoxFactors2D()),
-    randomBoxFactors3D_(initBoxFactors3D()),
-    iNextToLastPlane_(initBoxPlaneFactors()),
-    randomBox_
-    (
-        (isContinuous_ && Pstream::master())
+        (continuous_ && Pstream::master())
       ? dict.getOrDefault<List<List<scalar>>>
         (
-            "randomBox",
-            fillRandomBox() // First time-step
+            "box",
+            fillBox() // First time-step
         )
       :
-        fillRandomBox()
+        fillBox()
     ),
-    filterCoeffs_(computeFilterCoeffs()),
-    filteredRandomBox_
+    filterCoeffs_(calcFilterCoeffs()),
+    filteredBox_
     (
         pTraits<vector>::nComponents,
-        List<scalar>(planeDivisions_.first()*planeDivisions_.second(), Zero)
+        List<scalar>(n_.first()*n_.second(), Zero)
     ),
     U0_
     (
-        (variant_ == variantType::FORWARD_STEPWISE)
-      ? generateRandomSet<vectorField, vector>(patch().size())
+        fsm_
+      ? randomSet<vectorField, vector>(patch().size())
       : vectorField()
-    ),
-    computeVariant
-    (
-        (variant_ == variantType::FORWARD_STEPWISE)
-      ? &turbulentDigitalFilterInletFvPatchVectorField::computeReducedDFM
-      : &turbulentDigitalFilterInletFvPatchVectorField::computeDFM
     )
 {
-    if (isCorrectedFlowRate_)
+    if (correctFlowRate_)
     {
-        initialFlowRate_ = computeInitialFlowRate();
+        flowRate_ = calcFlowRate();
     }
 
     // Check if varying or fixed time-step computation
     if (db().time().isAdjustTimeStep())
     {
         WarningInFunction
-            << "Varying time-step computations are not fully supported"
-            << " for the moment."<< nl << nl;
+            << "  # Varying time-step computations are not fully supported #"
+            << endl;
     }
 
-    #ifdef FULLDEBUG
-    Info<< "Ends: Resource acquisition/initialisation for the"
-        << " synthetic turbulence boundary condition." << nl;
-    #endif
+    Info<< "  # turbulentDigitalFilterInlet initialisation completed #" << endl;
 }
 
 
@@ -1100,44 +878,40 @@ turbulentDigitalFilterInletFvPatchVectorField
 :
     fixedValueFvPatchField<vector>(ptf),
     mapperPtr_(nullptr),
-    variant_(ptf.variant_),
-    isGaussian_(ptf.isGaussian_),
-    isFixedSeed_(ptf.isFixedSeed_),
-    isContinuous_(ptf.isContinuous_),
-    isCorrectedFlowRate_(ptf.isCorrectedFlowRate_),
-    interpolateR_(ptf.interpolateR_),
-    interpolateUMean_(ptf.interpolateUMean_),
-    isInsideMesh_(ptf.isInsideMesh_),
-    isTaylorHypot_(ptf.isTaylorHypot_),
+    fsm_(ptf.fsm_),
+    Gaussian_(ptf.Gaussian_),
+    fixSeed_(ptf.fixSeed_),
+    continuous_(ptf.continuous_),
+    correctFlowRate_(ptf.correctFlowRate_),
+    interpR_(ptf.interpR_),
+    interpUMean_(ptf.interpUMean_),
     mapMethod_(ptf.mapMethod_),
     curTimeIndex_(ptf.curTimeIndex_),
-    tiny_(ptf.tiny_),
-    patchNormalSpeed_(ptf.patchNormalSpeed_),
-    modelConst_(ptf.modelConst_),
+    Ubulk_(ptf.Ubulk_),
+    C1_(ptf.C1_),
     perturb_(ptf.perturb_),
-    initialFlowRate_(ptf.initialFlowRate_),
+    flowRate_(ptf.flowRate_),
     rndGen_(ptf.rndGen_),
-    planeDivisions_(ptf.planeDivisions_),
+    n_(ptf.n_),
     invDelta_(ptf.invDelta_),
     indexPairs_(ptf.indexPairs_),
     R_(ptf.R_),
-    LundWuSquires_(ptf.LundWuSquires_),
+    Lund_(ptf.Lund_),
     UMean_(ptf.UMean_),
     Lbak_(ptf.Lbak_),
     L_(ptf.L_),
-    const1FSM_(ptf.const1FSM_),
-    const2FSM_(ptf.const2FSM_),
-    constList1FSM_(ptf.constList1FSM_),
-    constList2FSM_(ptf.constList2FSM_),
-    lenRandomBox_(ptf.lenRandomBox_),
-    randomBoxFactors2D_(ptf.randomBoxFactors2D_),
-    randomBoxFactors3D_(ptf.randomBoxFactors3D_),
+    C1FSM_(ptf.C1FSM_),
+    C2FSM_(ptf.C2FSM_),
+    coeffs1FSM_(ptf.coeffs1FSM_),
+    coeffs2FSM_(ptf.coeffs2FSM_),
+    szBox_(ptf.szBox_),
+    boxFactors2D_(ptf.boxFactors2D_),
+    boxFactors3D_(ptf.boxFactors3D_),
     iNextToLastPlane_(ptf.iNextToLastPlane_),
-    randomBox_(ptf.randomBox_),
+    box_(ptf.box_),
     filterCoeffs_(ptf.filterCoeffs_),
-    filteredRandomBox_(ptf.filteredRandomBox_),
-    U0_(ptf.U0_),
-    computeVariant(ptf.computeVariant)
+    filteredBox_(ptf.filteredBox_),
+    U0_(ptf.U0_)
 {}
 
 
@@ -1150,44 +924,40 @@ turbulentDigitalFilterInletFvPatchVectorField
 :
     fixedValueFvPatchField<vector>(ptf, iF),
     mapperPtr_(nullptr),
-    variant_(ptf.variant_),
-    isGaussian_(ptf.isGaussian_),
-    isFixedSeed_(ptf.isFixedSeed_),
-    isContinuous_(ptf.isContinuous_),
-    isCorrectedFlowRate_(ptf.isCorrectedFlowRate_),
-    interpolateR_(ptf.interpolateR_),
-    interpolateUMean_(ptf.interpolateUMean_),
-    isInsideMesh_(ptf.isInsideMesh_),
-    isTaylorHypot_(ptf.isTaylorHypot_),
+    fsm_(ptf.fsm_),
+    Gaussian_(ptf.Gaussian_),
+    fixSeed_(ptf.fixSeed_),
+    continuous_(ptf.continuous_),
+    correctFlowRate_(ptf.correctFlowRate_),
+    interpR_(ptf.interpR_),
+    interpUMean_(ptf.interpUMean_),
     mapMethod_(ptf.mapMethod_),
     curTimeIndex_(-1),
-    tiny_(ptf.tiny_),
-    patchNormalSpeed_(ptf.patchNormalSpeed_),
-    modelConst_(ptf.modelConst_),
+    Ubulk_(ptf.Ubulk_),
+    C1_(ptf.C1_),
     perturb_(ptf.perturb_),
-    initialFlowRate_(ptf.initialFlowRate_),
+    flowRate_(ptf.flowRate_),
     rndGen_(ptf.rndGen_),
-    planeDivisions_(ptf.planeDivisions_),
+    n_(ptf.n_),
     invDelta_(ptf.invDelta_),
     indexPairs_(ptf.indexPairs_),
     R_(ptf.R_),
-    LundWuSquires_(ptf.LundWuSquires_),
+    Lund_(ptf.Lund_),
     UMean_(ptf.UMean_),
     Lbak_(ptf.Lbak_),
     L_(ptf.L_),
-    const1FSM_(ptf.const1FSM_),
-    const2FSM_(ptf.const2FSM_),
-    constList1FSM_(ptf.constList1FSM_),
-    constList2FSM_(ptf.constList2FSM_),
-    lenRandomBox_(ptf.lenRandomBox_),
-    randomBoxFactors2D_(ptf.randomBoxFactors2D_),
-    randomBoxFactors3D_(ptf.randomBoxFactors3D_),
+    C1FSM_(ptf.C1FSM_),
+    C2FSM_(ptf.C2FSM_),
+    coeffs1FSM_(ptf.coeffs1FSM_),
+    coeffs2FSM_(ptf.coeffs2FSM_),
+    szBox_(ptf.szBox_),
+    boxFactors2D_(ptf.boxFactors2D_),
+    boxFactors3D_(ptf.boxFactors3D_),
     iNextToLastPlane_(ptf.iNextToLastPlane_),
-    randomBox_(ptf.randomBox_),
+    box_(ptf.box_),
     filterCoeffs_(ptf.filterCoeffs_),
-    filteredRandomBox_(ptf.filteredRandomBox_),
-    U0_(ptf.U0_),
-    computeVariant(ptf.computeVariant)
+    filteredBox_(ptf.filteredBox_),
+    U0_(ptf.U0_)
 {}
 
 
@@ -1204,7 +974,38 @@ void Foam::turbulentDigitalFilterInletFvPatchVectorField::updateCoeffs()
     {
         vectorField& U = *this;
 
-        computeVariant(this, U);
+        if (Pstream::master())
+        {
+            twoPointCorrs();
+            shiftRefill();
+        }
+
+        Pstream::scatter(filteredBox_);
+
+        mapFilteredBox(U);
+
+        if (fsm_)
+        {
+            //(XC:Eq. 14)
+            for (direction dir = 0; dir < pTraits<vector>::nComponents; ++dir)
+            {
+                U.component(dir) =
+                    U0_.component(dir)*coeffs1FSM_[dir]
+                  + U.component(dir)*coeffs2FSM_[dir];
+            }
+
+            U0_ = U;
+        }
+
+        onePointCorrs(U);
+
+        U += UMean_;
+
+        if (correctFlowRate_)
+        {
+            // (KCX:Eq. 8)
+            U *= (flowRate_/gSum(U & -patch().Sf()));
+        }
 
         curTimeIndex_ = db().time().timeIndex();
     }
@@ -1219,27 +1020,24 @@ void Foam::turbulentDigitalFilterInletFvPatchVectorField::write
 ) const
 {
     fixedValueFvPatchField<vector>::write(os),
-    os.writeEntry("variant", variantNames[variant_]);
-    os.writeEntry("planeDivisions", planeDivisions_);
+    os.writeEntry("n", n_);
     os.writeEntry("L", Lbak_);
 
-    if (!interpolateR_)
+    if (!interpR_)
     {
         R_.writeEntry("R", os);
     }
 
-    if (!interpolateUMean_)
+    if (!interpUMean_)
     {
         UMean_.writeEntry("UMean", os);
     }
 
-    os.writeEntryIfDifferent<bool>("isGaussian", true, isGaussian_);
-    os.writeEntryIfDifferent<bool>("isFixedSeed", true, isFixedSeed_);
-    os.writeEntryIfDifferent<bool>("isContinuous", false, isContinuous_);
-    os.writeEntryIfDifferent<bool>
-        ("isCorrectedFlowRate", true, isCorrectedFlowRate_);
-    os.writeEntryIfDifferent<bool>("isInsideMesh", false, isInsideMesh_);
-    os.writeEntryIfDifferent<bool>("isTaylorHypot", true, isTaylorHypot_);
+    os.writeEntryIfDifferent<bool>("fsm", false, fsm_);
+    os.writeEntryIfDifferent<bool>("Gaussian", true, Gaussian_);
+    os.writeEntryIfDifferent<bool>("fixSeed", true, fixSeed_);
+    os.writeEntryIfDifferent<bool>("continuous", false, continuous_);
+    os.writeEntryIfDifferent<bool>("correctFlowRate", true, correctFlowRate_);
 
     if (!mapMethod_.empty())
     {
@@ -1251,24 +1049,42 @@ void Foam::turbulentDigitalFilterInletFvPatchVectorField::write
         );
     }
 
-    os.writeEntry("threshold", tiny_);
-    os.writeEntry("patchNormalSpeed", patchNormalSpeed_);
-    os.writeEntry("modelConst", modelConst_);
+    os.writeEntry("Ubulk", Ubulk_);
+    os.writeEntry("C1", C1_);
     os.writeEntry("perturb", perturb_);
 
-    if (variant_ == variantType::FORWARD_STEPWISE)
+    if (fsm_)
     {
-        os.writeEntry("const1FSM", const1FSM_);
-        os.writeEntry("const2FSM", const2FSM_);
+        os.writeEntry("C1FSM", C1FSM_);
+        os.writeEntry("C2FSM", C2FSM_);
     }
 
-    if (isContinuous_ && Pstream::master())
+    if (continuous_ && Pstream::master())
     {
-        os.writeEntry("randomBox", randomBox_);
+        os.writeEntry("box", box_);
     }
 }
 
 
+void Foam::turbulentDigitalFilterInletFvPatchVectorField::autoMap
+(
+    const fvPatchFieldMapper& m
+)
+{
+    fixedValueFvPatchField<vector>::autoMap(m);
+}
+
+
+void Foam::turbulentDigitalFilterInletFvPatchVectorField::rmap
+(
+    const fvPatchVectorField& ptf,
+    const labelList& addr
+)
+{
+    fixedValueFvPatchField<vector>::rmap(ptf, addr);
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/turbulentDigitalFilterInletFvPatchVectorField.H b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/turbulentDigitalFilterInletFvPatchVectorField.H
index bbb0aa54f07b737d955557376b99ccd750fa45a5..3fe83e23713946afa653c37f020a58c54f5e8146 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/turbulentDigitalFilterInletFvPatchVectorField.H
+++ b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/turbulentDigitalFilterInletFvPatchVectorField.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -30,152 +30,179 @@ Group
     grpInletBoundaryConditions
 
 Description
-    Velocity boundary condition generating synthetic turbulence-alike
-    time-series for LES and DES turbulent flow computations.
-
-    To this end, two synthetic turbulence generators can be chosen:
-
-    - Digital-filter method-based generator (DFM)
-
-    \verbatim
-    Klein, M., Sadiki, A., and Janicka, J.
-        A digital filter based generation of inflow data for spatially
-        developing direct numerical or large eddy simulations,
-        Journal of Computational Physics (2003) 186(2):652-665.
-        doi:10.1016/S0021-9991(03)00090-1
-    \endverbatim
-
-    - Forward-stepwise method-based generator (FSM)
+    Digital-filter based boundary condition for velocity, i.e. \c U, to generate
+    synthetic turbulence-alike time-series for LES and DES turbulent flow
+    computations from input turbulence statistics.
 
+    References:
     \verbatim
-    Xie, Z.-T., and Castro, I.
-        Efficient generation of inflow conditions for large eddy simulation of
-        street-scale flows, Flow, Turbulence and Combustion (2008) 81(3):449-470
-        doi:10.1007/s10494-008-9151-5
+        Digital-filter method-based generator (DFM) (tag:KSJ):
+            Klein, M., Sadiki, A., & Janicka, J. (2003).
+            A digital filter based generation of inflow data for spatially
+            developing direct numerical or large eddy simulations.
+            Journal of computational Physics, 186(2), 652-665.
+            DOI:10.1016/S0021-9991(03)00090-1
+
+        Forward-stepwise method-based generator (FSM) (tag:XC)
+            Xie, Z. T., & Castro, I. P. (2008).
+            Efficient generation of inflow conditions for
+            large eddy simulation of street-scale flows.
+            Flow, turbulence and combustion, 81(3), 449-470.
+            DOI:10.1007/s10494-008-9151-5
+
+        Mass-inflow rate correction (tag:KCX):
+            Kim, Y., Castro, I. P., & Xie, Z. T. (2013).
+            Divergence-free turbulence inflow conditions for
+            large-eddy simulations with incompressible flow solvers.
+            Computers & Fluids, 84, 56-68.
+            DOI:10.1016/j.compfluid.2013.06.001
     \endverbatim
 
-    In DFM or FSM, a random number set (mostly white noise), and a group
+    In \c DFM or \c FSM, a random number set (mostly white noise), and a group
     of target statistics (mostly mean flow, Reynolds stress tensor profiles and
-    length-scale sets) are fused into a new number set (stochastic time-series,
+    length-scale sets) are merged into a new number set (stochastic time-series,
     yet consisting of the statistics) by a chain of mathematical operations
     whose characteristics are designated by the target statistics, so that the
     realised statistics of the new sets could match the target.
 
+    \verbatim
     Random number sets ---->-|
                              |
                          DFM or FSM ---> New stochastic time-series consisting
                              |           turbulence statistics
     Turbulence statistics ->-|
+    \endverbatim
 
-    The main difference between DFM and FSM is that the latter replaces the
-    streamwise convolution summation in DFM by a simpler and a quantitatively
-    justified equivalent procedure in order to reduce computational costs.
-    Accordingly, the latter potentially brings resource advantages for
-    computations involving relatively large length-scale sets and small
-    time-steps.
+    The main difference between \c DFM and \c FSM is that \c FSM replaces
+    the expensive-to-run streamwise convolution summation in \c DFM by a simpler
+    and an almost-equivalent-in-effect numerical procedure in order to reduce
+    computational costs. Accordingly, \c FSM potentially brings computational
+    resource advantages for computations involving relatively large streamwise
+    length-scale sets and small time-steps.
 
-    Synthetic turbulence is produced on a virtual rectangular structured-mesh
-    turbulence plane, which is parallel to the actual patch, and is mapped onto
-    the chosen patch by the selected mapping method.
+    Synthetic turbulence is generated on a virtual rectangular structured-mesh
+    plane, which is parallel to the chosen patch, and is mapped onto this patch
+    by the selected mapping method.
 
 Usage
-    \table
-     Property    | Description                     | Required    | Default value
-     planeDivisions   | Number of nodes on turbulence plane (e2, e3) [-] | yes |
-     L        | Integral length-scale set (9-comp):{e1,e2,e3}{u,v,w} [m] | yes |
-     R         | Reynolds stress tensor set (xx xy xz yy yz zz) [m2/s2]  | yes |
-     patchNormalSpeed    | Characteristic mean flow speed  [m/s] | yes |
-     isGaussian  | Autocorrelation function form                 | no | true
-     isFixedSeed | Flag to identify random-number seed is fixed  | no | true
-     isContinuous   | Flag for random-number restart behaviour   | no | false
-     isCorrectedFlowRate | Flag for mass flow rate correction    | no | true
-     interpolateR     | Placeholder flag: interpolate R field    | no | false
-     interpolateUMean | Placeholder flag: interpolate UMean field   | no | false
-     isInsideMesh | Placeholder flag: TP is inside mesh or on patch | no | false
-     isTaylorHypot | Placeholder flag: Taylor's hypothesis is on | no | true
-     mapMethod   | Method to map reference values             | no | nearestCell
-     threshold   | Threshold to avoid unintentional 'tiny' input | no | 1e-8
-     modelConst  | Model constant (Klein et al., Eq. 14)         | no | -0.5*PI
-     perturb     | Point perturbation for interpolation          | no | 1e-5
-     const1FSM | A model coefficient in FSM (Xie-Castro, Eq. 14) | no | -0.25*PI
-     const2FSM | A model coefficient in FSM (Xie-Castro, Eq. 14) | no | -0.5*PI
-    \endtable
-
-    Minimal example of the boundary condition specification with commented
-    options:
+    Example of the boundary condition specification:
     \verbatim
     <patchName>
     {
-        type                turbulentDigitalFilter;
-        variant             digitalFilter;          // reducedDigitalFilter;
-        planeDivisions      (<planeDivisionsHeight> <planeDivisionsWidth>);
-        L               (<Lxu> <Lxv> <Lxw> <Lyu> <Lyv> <Lyw> <Lzu> <Lzv> <Lzw>);
-        R                   (<Rxx> <Rxy> <Rxz> <Ryy> <Ryz> <Rzz>);
-        patchNormalSpeed    <characteristic flow speed>;
-        value               uniform (0 0 0);        // mandatory placeholder
-
-        // Optional entries with default input
-        isGaussian          true;           // false    // always false for FSM
-        isFixedSeed         true;           // false
-        isContinuous        false;          // true
-        isCorrectedFlowRate true;           // false
-        interpolateR        false;          // placeholder
-        interpolateUMean    false;          // placeholder
-        isInsideMesh        false;          // placeholder
-        isTaylorHypot       true;           // placeholder
-        mapMethod           nearestCell;    // planarInterpolation
-        threshold           1e-8;
-        modelConst          -1.5707;        //-0.5*PI;
+        // Mandatory entries (unmodifiable)
+        type                turbulentDigitalFilterInlet;
+        n                   (<nHeight> <nWidth>);
+        L                   (<L1> <L2> ... <L9>);
+        R                   uniform (<Rxx> <Rxy> <Rxz> <Ryy> <Ryz> <Rzz>);
+        UMean               uniform (1 0 0);
+        Ubulk               10.0;
+
+        // Optional entries (unmodifiable)
+        fsm                 false;
+        Gaussian            true;     // always false for FSM
+        fixSeed             true;
+        continuous          false;
+        correctFlowRate     true;
+        mapMethod           nearestCell;
         perturb             1e-5;
+        C1                  -1.5707;  //-0.5*PI;
+        C1FSM               -0.7854   //-0.25*PI;
+        C2FSM               -1.5707;  //-0.5*PI;
 
-        // Optional entries for only FSM with default input
-        const1FSM           -0.7854         //-0.25*PI;
-        const2FSM           -1.5707;        //-0.5*PI;
+        // Optional (inherited) entries
+        ...
     }
     \endverbatim
 
-    Among the dictionary entries, two entries can be input as patch profiles:
+    where the entries mean:
+    \table
+      Property | Description                            | Type  | Req'd | Dflt
+      type     | Type name: turbulentDigitalFilterInlet | word  | yes   | -
+      n        | Number of cells on turbulence generation plane <!--
+                                           --> | tuple of labels | yes   | -
+      L        | Integral length-scale set <!--
+              --> (Lxu Lxv Lxw Lyu Lyv Lyw Lzu Lzv Lzw) [m] | tensor | yes | -
+      R        | Reynolds stress tensor set <!--
+              --> (xx xy xz yy yz zz) [m2/s2]      | symmTensorField | yes | -
+      UMean    | Mean velocity profile [m/s]       | vectorField     | yes | -
+      Ubulk    | Characteristic patch-normal bulk flow speed  [m/s] <!--
+                                                        --> | scalar | yes | -
+      fsm      | Flag to turn on the forward-stepwise method | bool | no | false
+      Gaussian | Autocorrelation function form         | bool | no | true
+      fixSeed  | Flag to fix random-number generator seed to 1234 <!--
+             --> or generate a new seed based on clock-time per simulation <!--
+                                                     --> | bool | no | true
+      continuous | Flag to write random-number sets at output time, <!--
+             --> and to read them on restart. Otherwise, generate <!--
+             --> new random-number sets of restart       | bool | no  false
+      correctFlowRate | Flag to correct mass-inflow rate on turbulence <!--
+                    --> plane in (only) streamwise direction | bool | no | true
+      mapMethod  | Interpolation-to-patch method      | word | no | nearestCell
+      perturb    | Point perturbation for planarInterpolation mapMethod <!--
+                                                  --> | scalar | no | 1e-5
+      C1         | Model constant shaping autocorrelation function <!--
+                                     --> (KSJ:Eq. 14) | scalar | no | -0.5*PI
+      C1FSM | Model coefficient in FSM (XC:Eq. 14)    | scalar | no | -0.25*PI
+      C2FSM | Model coefficient in FSM (XC:Eq. 14)    | scalar | no | -0.5*PI
+    \endtable
+
+    The inherited entries are elaborated in:
+      - \link fixedValueFvPatchFields.H \endlink
 
+    Options for the \c fsm entry:
     \verbatim
-     - Reynolds stress tensor, R
-     - Mean velocity, UMean
+      false | Method due to (KSJ)
+      true  | Method due to (XC)
     \endverbatim
 
-    Profile data and corresponding coordinates are then input in the following
-    directories:
-
+    Options for the \c Gaussian entry:
     \verbatim
-     - $FOAM_CASE/constant/boundaryData/\<patchName\>/points
-     - $FOAM_CASE/constant/boundaryData/\<patchName\>/0/\{R|UMean\}
+      true  | Gaussian function
+      false | Exponential function (only option for FSM)
     \endverbatim
 
-    The profile data and corresponding coordinates take the same form used by
-    the \c timeVaryingMappedFixedValue and \c turbulentDFSEMInlet boundary
-    conditions, consisting of a \c points file containing a list of 3-D
-    coordinates, and profile data files providing a value per coordinate.
+    Options for the \c mapMethod entry:
+    \verbatim
+      nearestCell         | One-to-one direct map, no interpolation
+      planarInterpolation | Bilinear interpolation
+    \endverbatim
 
-    Reynolds stress tensor and mean velocity input are in the global coordinate
-    system whereas integral length scale set input is in the local patch
-    coordinate system.
+    Patch-profile input is available for two entries:
+    \verbatim
+      R     | Reynolds stress tensor
+      UMean | Mean velocity
+    \endverbatim
+    where the input profiles and profile coordinates are located in:
+    \verbatim
+      Coordinates | $FOAM_CASE/constant/boundaryData/\<patchName\>/points
+      R/UMean     | $FOAM_CASE/constant/boundaryData/\<patchName\>/0/\{R/UMean\}
+    \endverbatim
 
-    It is assumed that the patch normal direction is \c e1, and the remaining
-    patch plane directions are \c e2 and \c e3 following the right-handed
-    coordinate system, which should be taken into consideration when the
-    integral length scale set is input. The first three integral scale entries,
-    i.e. L_e1u, Le1v, Le2w, should always correspond to the length scales
-    that are in association with the convective mean flow direction.
+    \c points file contains a list of three-dimensional coordinates, and
+    profile data files provide a value corresponding to each coordinate.
 
 Note
-    - \c mapMethod \c planarInterpolation option requires point coordinates
-    which can form a plane, thus input point coordinates varying only in a
-    single direction will trigger error.
-    - \c adjustTimeStep = true option is not fully supported at the moment.
-
-SeeAlso
-    turbulentDFSEMInletFvPatchVectorField.C
+    - \c mapMethod=planarInterpolation option needs point coordinates that can
+    form a plane.
+    - \c adjustTimeStep=true option is currently not fully supported.
+    - In order to obtain Reynolds stress tensor information, experiments, RANS
+    simulations or engineering relations can be used.
+    - \c continuous=true means deterministic-statistical consistent restart
+    (relatively more expensive), and \c continuous=false means deterministic
+    discontinuity in synthetic turbulence time-series by keeping statistical
+    consistency (relatively cheaper).
+    - For \c L, the first three entries should always correspond to the
+    length scales in association with the convective (streamwise) mean flow
+    direction.
+    - Streamwise integral length scales are converted to integral time scales
+    by using Taylor's frozen turbulence hypothesis, and \c Ubulk.
+
+See also
+    - turbulentDFSEMInletFvPatchVectorField.C
 
 SourceFiles
     turbulentDigitalFilterInletFvPatchVectorField.C
+    turbulentDigitalFilterInletFvPatchVectorFieldTemplates.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -184,7 +211,6 @@ SourceFiles
 
 #include "fixedValueFvPatchFields.H"
 #include "Random.H"
-#include <functional>
 #include "fieldTypes.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -202,172 +228,125 @@ class turbulentDigitalFilterInletFvPatchVectorField
 :
     public fixedValueFvPatchVectorField
 {
-
-    // Private Enumerations
-
-        //- Options for the synthetic turbulence generator variant
-        enum variantType : uint8_t
-        {
-            DIGITAL_FILTER = 1,    //!< Digital-filter method (Klein et al.)
-            FORWARD_STEPWISE = 2,  //!< Forward-stepwise method (Xie-Castro)
-        };
-
-        //- Names for variant types
-        static const Enum<variantType> variantNames;
-
-
     // Private Data
 
-        //- 2D interpolation (for 'planarInterpolation' mapMethod)
+        //- Bilinear interpolation (for 'mapMethod=planarInterpolation')
         mutable autoPtr<pointToPointPlanarInterpolation> mapperPtr_;
 
-        //- Selected option for the synthetic turbulence generator variant
-        const enum variantType variant_;
-
-        //- Flag: correlation function form is Gaussian or Exponential
-        //  for variantType::DIGITAL_FILTER, default=Gaussian
-        //  for variantType::FORWARD_STEPWISE, default=Exponential (only option)
-        const bool isGaussian_;
+        //- Flag to enable the forward-stepwise method
+        const bool fsm_;
 
-        //- Flag: random-number generator seed is fixed (default=true) or
-        //- generated pseudo-randomly based on clock-time per simulation
-        const bool isFixedSeed_;
+        //- Flag to select correlation function form: Gaussian or exponential
+        const bool Gaussian_;
 
-        //- Flag: write non-manipulated random-number sets at output time, and
-        //- to read them on restart. Otherwise, generate new random-number sets
-        //- on restart. (default=false)
-        //  true: deterministic & statistically consistent, more expensive
-        //  false: deterministic discontinuity & statistically consistent, cheaper
-        const bool isContinuous_;
+        //- Flag to fix the random-number generator seed to 1234 or
+        //- generate a new seed based on clock-time per simulation
+        const bool fixSeed_;
 
-        //- Flag: mass flow rate is corrected on turbulence plane (default=true)
-        const bool isCorrectedFlowRate_;
+        //- Flag to write random-number sets at output time, and to read them
+        //- on restart. Otherwise, generate new random-number sets on restart
+        const bool continuous_;
 
-        //- Flag: interpolate R field (default=false)
-        bool interpolateR_;
+        //- Flag to correct mass flow rate on turbulence plane
+        const bool correctFlowRate_;
 
-        //- Flag: interpolate UMean field (default=false)
-        bool interpolateUMean_;
+        //- Internal flag to read R from data files
+        bool interpR_;
 
-        //- Flag: turbulence plane is inside mesh or on a patch (default=false)
-        //  Currently, true option is not available.
-        const bool isInsideMesh_;
-
-        //- Flag: convert streamwise (x) length scales to time scales by
-        //- Taylor's 'frozen turbulence' hypothesis (default=true)
-        //  Currently, false option is not available.
-        const bool isTaylorHypot_;
+        //- Internal flag to read UMean from data files
+        bool interpUMean_;
 
         //- Method for interpolation between a patch and turbulence plane
-        //- (default=nearestCell)
-        //  Options:
-        //   - nearestCell: one-to-one direct map, no interpolation
-        //   - planarInterpolation: bilinear interpolation
         const word mapMethod_;
 
         //- Current time index
         label curTimeIndex_;
 
-        //- Threshold to avoid unintentional 'tiny' input scalars (default=1e-8)
-        const scalar tiny_;
-
-        //- Characteristic (e.g. bulk) mean speed of flow in the patch normal
-        //- direction [m/s]
-        const scalar patchNormalSpeed_;
+        //- Characteristic patch-normal bulk flow speed [m/s]
+        const scalar Ubulk_;
 
-        //- Model constant shaping autocorr function [-] (default='-0.5*pi')
-        //  (Klein et al., 2003, Eq. 14)
-        const scalar modelConst_;
+        //- Model constant shaping autocorrelation function (KSJ:Eq. 14)
+        const scalar C1_;
 
-        //- Fraction of perturbation (fraction of bounding box) to add (for
-        //- 'planarInterpolation' mapMethod)
+        //- Fraction of perturbation (fraction of bounding box) to add
         const scalar perturb_;
 
-        //- Initial (first time-step) mass/vol flow rate [m^3/s]
-        scalar initialFlowRate_;
+        //- First time-step mass/volumetric flow rate
+        scalar flowRate_;
 
         //- Random number generator
         Random rndGen_;
 
-        //- Number of nodes on turbulence plane (e2 e3) [-]
-        const Tuple2<label, label> planeDivisions_;
+        //- Number of cells on turbulence plane (<nHeight> <nWidth>) [-]
+        const Tuple2<label, label> n_;
 
-        //- Turbulence plane mesh size (reversed) (e2 e3) [1/m]
+        //- Uniform mesh size on turbulence plane (reversed) [1/m]
         Vector2D<scalar> invDelta_;
 
-        //- Nearest cell mapping: Index pairs between patch and turbulence plane
+        //- Index pairs between patch and turbulence plane
+        //- for the nearest cell mapping
         const List<Pair<label>> indexPairs_;
 
-        //- Reynolds stress tensor profile (xx xy xz yy yz zz) in global
-        //- coordinates [m^2/s^2]
-        symmTensorField R_;
+        //- Reynolds stress tensor profile field in global coordinates [m2/s2]
+        const symmTensorField R_;
 
-        //- Lund-Wu-Squires transformation (Cholesky decomp.) [m/s]
+        //- Lund-Wu-Squires transformed R field (Cholesky decomp.) [m/s]
         //- Mapped onto actual mesh patch rather than turbulence plane
-        //  (Klein et al., 2003, Eq. 5)
-        symmTensorField LundWuSquires_;
+        //  (KSJ:Eq. 5)
+        const symmTensorField Lund_;
 
-        //- Mean inlet velocity profile in global coordinates [m/s]
+        //- Mean inlet velocity profile field in global coordinates [m/s]
         vectorField UMean_;
 
         //- Integral length-scale set per turbulence plane section in local
         //- coordinates (e1u, e1v, e1w, e2u, e2v, e2w, e3u, e3v, e3w) [m]
-        //  First three entries should always correspond to the length scales
-        //  in association with the convective mean flow direction
         //  Backup of L_ for restart purposes
         const tensor Lbak_;
 
-        //- Integral length-scale set in mesh units [node]
+        //- Integral length-scale set in mesh units [cell]
         const tensor L_;
 
         //- One of the two model coefficients in FSM
-        //  (Xie-Castro, 2008, the argument of the first exp func in Eq. 14)
-        const scalar const1FSM_;
+        //  (XC:The argument of the first exp func in Eq. 14)
+        const scalar C1FSM_;
 
         //- One of the two model coefficients in FSM
-        //  (Xie-Castro, 2008, the argument of the second exp func in Eq. 14)
-        const scalar const2FSM_;
+        //  (XC:The argument of the second exp func in Eq. 14)
+        const scalar C2FSM_;
 
         //- One of the two exponential functions in FSM
-        //  (Xie-Castro, 2008, the first exponential function in Eq. 14)
-        const List<scalar> constList1FSM_;
+        //  (XC:The first exponential function in Eq. 14)
+        const List<scalar> coeffs1FSM_;
 
         //- One of the two exponential functions in FSM
-        //  (Xie-Castro, 2008, the first exponential function in Eq. 14)
-        const List<scalar> constList2FSM_;
+        //  (XC:The first exponential function in Eq. 14)
+        const List<scalar> coeffs2FSM_;
 
-        //- Number of nodes in random-number box [node]
+        //- Number of cells in random-number box [cell]
         //- Random-number sets within box are filtered with filterCoeffs_
-        //- (e1u, e1v, e1w, e2u, e2v, e2w, e3u, e3v, e3w)
-        const List<label> lenRandomBox_;
+        const List<label> szBox_;
 
-        //- Convenience factors for 2-D random-number box [-]
-        const List<label> randomBoxFactors2D_;
+        //- Convenience factors for two-dimensional random-number box [-]
+        const List<label> boxFactors2D_;
 
-        //- Convenience factors for 3-D randomNum box [-]
-        const List<label> randomBoxFactors3D_;
+        //- Convenience factors for three-dimensional random-number box [-]
+        const List<label> boxFactors3D_;
 
         //- Index to the first elem of last plane of random-number box [-]
         const List<label> iNextToLastPlane_;
 
-        //- Random-number sets distributed over a 3-D box (u, v, w)
-        List<List<scalar>> randomBox_;
+        //- Random-number sets distributed over three-dimensional box (u, v, w)
+        List<List<scalar>> box_;
 
-        //- Filter coefficients corresponding to L_ [-]
+        //- Filter coefficients corresponding to L [-]
         const List<List<scalar>> filterCoeffs_;
 
-        //- Filter-applied random-number sets [m/s] (effectively turb plane)
-        List<List<scalar>> filteredRandomBox_;
+        //- Filter-applied random-number sets [m/s], i.e. turbulence plane
+        List<List<scalar>> filteredBox_;
 
         //- Filter-applied previous-time-step velocity field [m/s] used in FSM
         vectorField U0_;
 
-        //- Run-time function selector between the original and reduced methods
-        const std::function
-        <
-            void(turbulentDigitalFilterInletFvPatchVectorField*, vectorField&)
-        > computeVariant;
-
 
     // Private Member Functions
 
@@ -393,82 +372,59 @@ class turbulentDigitalFilterInletFvPatchVectorField
 
         //- Generate random-number sets obeying the standard normal distribution
         template<class Form, class Type>
-        Form generateRandomSet(const label len);
+        Form randomSet(const label len);
 
         //- Compute nearest cell index-pairs between turbulence plane and patch
-        List<Pair<label>> patchIndexPairs();
+        List<Pair<label>> indexPairs();
 
-        //- Check R_ (mapped on actual mesh) for mathematical domain errors
-        void checkRTensorRealisable() const;
+        //- Check R on patch for mathematical domain errors
+        void checkR() const;
 
         //- Compute Lund-Wu-Squires transformation
-        //  (Klein et al., 2003, Eq. 5)
-        symmTensorField computeLundWuSquires() const;
+        symmTensorField calcLund() const;
 
-        //- Compute patch-normal into the domain
-        vector computePatchNormal() const;
+        //- Compute the first time-step mass/
+        //- volumetric flow rate based on UMean
+        scalar calcFlowRate() const;
 
-        //- Compute initial (first time-step) mass/vol flow rate based on UMean_
-        scalar computeInitialFlowRate() const;
-
-        //- Convert streamwise integral length scales to integral time scales
-        //- via Taylor's frozen turbulence hypothesis
-        void convertToTimeScale(tensor& L) const;
-
-        //- Convert length scale phys. unit to turbulence plane mesh-size unit
-        //- (Klein et al., 2003, Eq. 13)
-        tensor convertScalesToGridUnits(const tensor& L) const;
+        //- Convert integral length scales in meters
+        //- to turbulence plane cell-size units
+        tensor meterToCell(const tensor& L) const;
 
         //- Resource allocation functions for the convenience factors
-        List<label> initLenRandomBox() const;
-        List<label> initBoxFactors2D() const;
-        List<label> initBoxFactors3D() const;
-        List<label> initBoxPlaneFactors() const;
+        List<label> initBox() const;
+        List<label> initFactors2D() const;
+        List<label> initFactors3D() const;
+        List<label> initPlaneFactors() const;
 
         //- Compute various convenience factors for random-number box
-        List<List<scalar>> fillRandomBox();
+        List<List<scalar>> fillBox();
 
         //- Compute filter coeffs once per simulation
-        //  (Klein et al., 2003, Eq. 14)
-        List<List<scalar>> computeFilterCoeffs() const;
+        List<List<scalar>> calcFilterCoeffs() const;
 
         //- Discard current time-step random-box plane (closest to patch) by
-        //- shifting from the back to the front, and dd new plane to the back
-        void rndShiftRefill();
-
-        //- Map two-point correlated random-number sets on patch based on chosen
-        //- mapping method
-        void mapFilteredRandomBox(vectorField& U);
+        //- shifting from the back to the front, and add new plane to the back
+        void shiftRefill();
 
-        //- Map R_ on patch
-        void embedOnePointCorrs(vectorField& U) const;
+        //- Map two-point correlated random-number sets
+        //- on patch based on chosen mapping method
+        void mapFilteredBox(vectorField& U);
 
-        //- Map UMean_ on patch
-        void embedMeanVelocity(vectorField& U) const;
+        //- Embed one-point correlations, i.e. R, on patch
+        void onePointCorrs(vectorField& U) const;
 
-        //- Correct mass/vol flow rate in (only) streamwise direction
-        void correctFlowRate(vectorField& U) const;
+        //- Embed two-point correlations, i.e. L, on box
+        //  Three-dimensional "valid"-type separable
+        //  convolution summation algorithm
+        //  (Based on Song Ho Ahn's two-dimensional "full"-type convolution)
+        void twoPointCorrs();
 
-        //- 3-D 'valid'-type 'separable' convolution summation algorithm
-        //  'Inspired' from Song Ho Ahn's 2-D 'full'-type convolution algorithm
-        //  with his permission
-        void embedTwoPointCorrs();
+        //- Compute coeffs1FSM_ once per simulation
+        List<scalar> calcCoeffs1FSM() const;
 
-        //- Compute the DFM (Klein et al., 2003)
-        void computeDFM(vectorField& U);
-
-        //- Compute the reduced DFM (i.e. hybrid FSM-DFM) (Xie-Castro, 2008)
-        void computeReducedDFM(vectorField& U);
-
-        //- Compute constList1FSM_ once per simulation
-        List<scalar> computeConstList1FSM() const;
-
-        //- Compute constList2FSM_ once per simulation
-        List<scalar> computeConstList2FSM() const;
-
-        //- Compute the forward-stepwise method
-        //  (Xie-Castro, 2008, Eq. 14)
-        void computeFSM(vectorField& U);
+        //- Compute coeffs2FSM_ once per simulation
+        List<scalar> calcCoeffs2FSM() const;
 
 
 public:
@@ -476,6 +432,7 @@ public:
    //- Runtime type information
    TypeName("turbulentDigitalFilterInlet");
 
+
     // Constructors
 
         //- Construct from patch and internal field
@@ -550,6 +507,19 @@ public:
             virtual void updateCoeffs();
 
 
+        // Mapping functions
+
+            //- Map (and resize as needed) from self given a mapping object
+            virtual void autoMap(const fvPatchFieldMapper& m);
+
+            //- Reverse map the given fvPatchField onto this fvPatchField
+            virtual void rmap
+            (
+                const fvPatchVectorField& ptf,
+                const labelList& addr
+            );
+
+
        //- Write
        virtual void write(Ostream&) const;
 };
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/turbulentDigitalFilterInletFvPatchVectorFieldTemplates.C b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/turbulentDigitalFilterInletFvPatchVectorFieldTemplates.C
index 03a17b787fe01c9c404b01ef9db629a9cf1238c5..706d7cd72040144b09a40ddd0f3fc54c11bc4955 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/turbulentDigitalFilterInletFvPatchVectorFieldTemplates.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/turbulentDigitalFilterInletFvPatchVectorFieldTemplates.C
@@ -5,8 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2015 OpenFOAM Foundation
-    Copyright (C) 2016 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -123,7 +122,7 @@ Foam::turbulentDigitalFilterInletFvPatchVectorField::interpolateBoundaryData
     const rawIOField<Type> vals(io, false);
 
 
-    Info<< "Turbulent DFM/FSM patch " << patchName
+    Info<< "turbulentDigitalFilterInlet patch " << patchName
         << ": Interpolating field " << fieldName
         << " from " << valsFile << endl;
 
@@ -132,7 +131,7 @@ Foam::turbulentDigitalFilterInletFvPatchVectorField::interpolateBoundaryData
 
 
 template<class Form, class Type>
-Form Foam::turbulentDigitalFilterInletFvPatchVectorField::generateRandomSet
+Form Foam::turbulentDigitalFilterInletFvPatchVectorField::randomSet
 (
     const label len
 )
@@ -149,4 +148,5 @@ Form Foam::turbulentDigitalFilterInletFvPatchVectorField::generateRandomSet
     return randomSet;
 }
 
+
 // ************************************************************************* //
diff --git a/src/functionObjects/field/Make/files b/src/functionObjects/field/Make/files
index acf12cb0a4b5c5d9719ad4a554ca7f302eb78884..fecf17a751472ba4b7c14247265d1f3d35e6ca23 100644
--- a/src/functionObjects/field/Make/files
+++ b/src/functionObjects/field/Make/files
@@ -79,6 +79,8 @@ pressure/pressure.C
 MachNo/MachNo.C
 Curle/Curle.C
 reference/reference.C
+log/log.C
+pow/pow.C
 
 fieldsExpression/fieldsExpression.C
 add/add.C
diff --git a/src/functionObjects/field/log/log.C b/src/functionObjects/field/log/log.C
new file mode 100644
index 0000000000000000000000000000000000000000..70ff35eace83dca0c69ff070a971f6c8e29735f1
--- /dev/null
+++ b/src/functionObjects/field/log/log.C
@@ -0,0 +1,124 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2018-2019 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
+------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "log.H"
+#include "volFields.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+    defineTypeNameAndDebug(log, 0);
+    addToRunTimeSelectionTable(functionObject, log, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+bool Foam::functionObjects::log::calc()
+{
+    if (foundObject<volScalarField>(fieldName_))
+    {
+        const volScalarField& x = lookupObject<volScalarField>(fieldName_);
+
+        // Cache the current debug setting for dimensionSet
+        const bool dimensionSetDebug = dimensionSet::debug;
+
+        // Switch-off dimension checking if requested
+        if (!checkDimensions_)
+        {
+            dimensionSet::debug = 0;
+        }
+
+        bool stored = store
+        (
+            resultName_,
+            scale_*Foam::log(max(x, clipValue_)) + offset_
+        );
+
+        // Reinstate dimension checking
+        if (!checkDimensions_)
+        {
+            dimensionSet::debug = dimensionSetDebug;
+        }
+
+        return stored;
+    }
+
+    return false;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::functionObjects::log::log
+(
+    const word& name,
+    const Time& runTime,
+    const dictionary& dict
+)
+:
+    fieldExpression(name, runTime, dict, typeName),
+    checkDimensions_(true),
+    clipValue_(SMALL),
+    scale_(1.0),
+    offset_(0.0)
+{
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::functionObjects::log::read(const dictionary& dict)
+{
+    if (fvMeshFunctionObject::read(dict) && fieldExpression::read(dict))
+    {
+        checkDimensions_ = dict.getOrDefault<Switch>("checkDimensions", true);
+        clipValue_ =
+            dict.getCheckOrDefault<scalar>
+            (
+                "clip",
+                SMALL,
+                scalarMinMax::ge(SMALL)
+            );
+        scale_ = dict.getOrDefault<scalar>("scale", 1.0);
+        offset_ = dict.getOrDefault<scalar>("offset", 0.0);
+
+        return true;
+    }
+
+    return false;
+}
+
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/log/log.H b/src/functionObjects/field/log/log.H
new file mode 100644
index 0000000000000000000000000000000000000000..c3eb23114d8a8fab801e30222eb973765e0400be
--- /dev/null
+++ b/src/functionObjects/field/log/log.H
@@ -0,0 +1,205 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2018-2019 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::functionObjects::log
+
+Group
+    grpFieldFunctionObjects
+
+Description
+    Computes the natural logarithm of an input \c volScalarField.
+
+    \f[
+        f = s \ln(max(f_0, a)) + t
+    \f]
+
+    where
+    \vartable
+      f     | Output volScalarField
+      f_0   | Input volScalarField
+      \ln   | Natural logarithm operator
+      a     | Clip scalar
+      s     | Scaling factor
+      t     | Offset factor
+    \endvartable
+
+    \table
+      Operand       | Type           | Location
+      input         | volScalarField | $FOAM_CASE/\<time\>/\<inpField\>
+      output file   | -              | -
+      output field  | volScalarField | $FOAM_CASE/\<time\>/\<outField\>
+    \endtable
+
+Usage
+    Minimal example by using \c system/controlDict.functions:
+    \verbatim
+    log1
+    {
+        // Mandatory entries (unmodifiable)
+        type            log;
+        libs            (fieldFunctionObjects);
+
+        // Mandatory (inherited) entry (runtime modifiable)
+        field           <inpField>;
+
+        // Optional entries (runtime modifiable)
+        clip            1e-3;
+        checkDimensions false;
+        scale           1.0;
+        offset          0.0;
+
+        // Optional (inherited) entries
+        ...
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property     | Description                        | Type | Req'd | Dflt
+      type         | Type name: log                     | word |  yes  | -
+      libs         | Library name: fieldFunctionObjects | word |  yes  | -
+      field        | Name of the operand field          | word |  yes  | -
+      clip         | Value to clip the operand field values <!--
+                 --> to prevent zero or negative input  | scalar | no  | SMALL
+      checkDimensions | Flag to check dimensions of the operand field <!--
+                                                    --> | bool   | no  | true
+      scale        | Scaling factor - \c s above        | scalar | no  | 1.0
+      offset       | Offset factor - \c t above         | scalar | no  | 0.0
+    \endtable
+
+    The inherited entries are elaborated in:
+     - \link functionObject.H \endlink
+     - \link fieldExpression.H \endlink
+
+    Minimal example by using the \c postProcess utility:
+    \verbatim
+        postProcess -func "log(<inpField>)" -scale 1.0 -offset 0.0
+    \endverbatim
+
+Note
+    - Performs \f$\ln(max(x, a))\f$ where \f$x\f$ is a \c volScalarField, and
+    \f$a\f$ a clip scalar, equals to \c SMALL by default. This prevents zero or
+    negative input \f$x\f$, hence the domain error in the natural logarithm.
+    - Dimension checking can optionally be suspended if \f$x\f$ is dimensioned.
+
+See also
+    - Foam::functionObject
+    - Foam::functionObjects::fieldExpression
+    - Foam::functionObjects::fvMeshFunctionObject
+    - ExtendedCodeGuide::functionObjects::field::log
+
+SourceFiles
+    log.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef functionObjects_log_H
+#define functionObjects_log_H
+
+#include "fieldExpression.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class log Declaration
+\*---------------------------------------------------------------------------*/
+
+class log
+:
+    public fieldExpression
+{
+    // Private Data
+
+        //- Flag to check dimensions of the operand
+        Switch checkDimensions_;
+
+        //- Value to clip the operand field values
+        //- to prevent zero or negative input
+        scalar clipValue_;
+
+        //- Scaling factor
+        scalar scale_;
+
+        //- Offset factor
+        scalar offset_;
+
+
+    // Private Member Functions
+
+        //- Calculate the log field and return true if successful
+        virtual bool calc();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("log");
+
+
+    // Constructors
+
+        //- Construct from Time and dictionary
+        log
+        (
+            const word& name,
+            const Time& runTime,
+            const dictionary& dict
+        );
+
+        //- No copy construct
+        log(const log&) = delete;
+
+        //- No copy assignment
+        void operator=(const log&) = delete;
+
+
+    //- Destructor
+    virtual ~log() = default;
+
+
+    // Member Functions
+
+        //- Read the randomise data
+        virtual bool read(const dictionary&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace functionObjects
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/pow/pow.C b/src/functionObjects/field/pow/pow.C
new file mode 100644
index 0000000000000000000000000000000000000000..6315653f93c28d20278ce540dc662efb4cb89ec4
--- /dev/null
+++ b/src/functionObjects/field/pow/pow.C
@@ -0,0 +1,117 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "pow.H"
+#include "volFields.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+    defineTypeNameAndDebug(pow, 0);
+    addToRunTimeSelectionTable(functionObject, pow, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+bool Foam::functionObjects::pow::calc()
+{
+    if (foundObject<volScalarField>(fieldName_))
+    {
+        const volScalarField& x = lookupObject<volScalarField>(fieldName_);
+
+        // Cache the current debug setting for dimensionSet
+        const bool dimensionSetDebug = dimensionSet::debug;
+
+        // Switch-off dimension checking if requested
+        if (!checkDimensions_)
+        {
+            dimensionSet::debug = 0;
+        }
+
+        bool stored = store
+        (
+            resultName_,
+            scale_*Foam::pow(x, n_) + offset_
+        );
+
+        // Reinstate dimension checking
+        if (!checkDimensions_)
+        {
+            dimensionSet::debug = dimensionSetDebug;
+        }
+
+        return stored;
+    }
+
+    return false;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::functionObjects::pow::pow
+(
+    const word& name,
+    const Time& runTime,
+    const dictionary& dict
+)
+:
+    fieldExpression(name, runTime, dict, typeName),
+    checkDimensions_(true),
+    n_(0.0),
+    scale_(1.0),
+    offset_(0.0)
+{
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::functionObjects::pow::read(const dictionary& dict)
+{
+    if (fvMeshFunctionObject::read(dict) && fieldExpression::read(dict))
+    {
+        checkDimensions_ = dict.getOrDefault<Switch>("checkDimensions", true);
+        n_ = dict.get<scalar>("n");
+        scale_ = dict.getOrDefault<scalar>("scale", 1.0);
+        offset_ = dict.getOrDefault<scalar>("offset", 0.0);
+
+        return true;
+    }
+
+    return false;
+}
+
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/pow/pow.H b/src/functionObjects/field/pow/pow.H
new file mode 100644
index 0000000000000000000000000000000000000000..8820a60b20873f3da442fd46c25d2e7d124c8bf0
--- /dev/null
+++ b/src/functionObjects/field/pow/pow.H
@@ -0,0 +1,197 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::functionObjects::pow
+
+Group
+    grpFieldFunctionObjects
+
+Description
+    Computes the power of an input \c volScalarField.
+
+    \f[
+        f = s f_0^n + t
+    \f]
+
+    where
+    \vartable
+        f     | Output volScalarField
+        f_0   | Input volScalarField
+        n     | Exponent
+        s     | Scaling factor
+        t     | Offset factor
+    \endvartable
+
+    \table
+      Operand       | Type           | Location
+      input         | volScalarField | $FOAM_CASE/\<time\>/\<inpField\>
+      output file   | -              | -
+      output field  | volScalarField | $FOAM_CASE/\<time\>/\<outField\>
+    \endtable
+
+Usage
+    Minimal example by using \c system/controlDict.functions:
+    \verbatim
+    pow1
+    {
+        // Mandatory entries (unmodifiable)
+        type            pow;
+        libs            (fieldFunctionObjects);
+
+        // Mandatory (inherited) entry (runtime modifiable)
+        field           <inpField>;
+
+        // Mandatory entry (runtime modifiable)
+        n               0.25;
+
+        // Optional entries (runtime modifiable)
+        checkDimensions false;
+        scale           1.0;
+        offset          0.0;
+
+        // Optional (inherited) entries
+        ...
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property     | Description                        | Type | Req'd | Dflt
+      type         | Type name: pow                     | word |  yes  | -
+      libs         | Library name: fieldFunctionObjects | word |  yes  | -
+      field        | Name of the operand field          | word |  yes  | -
+      n            | Exponent                           | scalar | yes | -
+      checkDimensions | Flag to check dimensions of the operand field <!--
+                                                    --> | bool   | no  | true
+      scale        | Scaling factor - \c s above        | scalar | no  | 1.0
+      offset       | Offset factor - \c t above         | scalar | no  | 0.0
+    \endtable
+
+    The inherited entries are elaborated in:
+     - \link functionObject.H \endlink
+     - \link fieldExpression.H \endlink
+
+    Minimal example by using the \c postProcess utility:
+    \verbatim
+        postProcess -func "pow(<inpField>)" -scale 1.0 -offset 0.0
+    \endverbatim
+
+See also
+    - Foam::functionObject
+    - Foam::functionObjects::fieldExpression
+    - Foam::functionObjects::fvMeshFunctionObject
+    - ExtendedCodeGuide::functionObjects::field::pow
+
+SourceFiles
+    pow.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef functionObjects_pow_H
+#define functionObjects_pow_H
+
+#include "fieldExpression.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class pow Declaration
+\*---------------------------------------------------------------------------*/
+
+class pow
+:
+    public fieldExpression
+{
+    // Private Data
+
+        //- Flag to check dimensions of the operand
+        Switch checkDimensions_;
+
+        //- Exponent
+        scalar n_;
+
+        //- Scaling factor
+        scalar scale_;
+
+        //- Offset factor
+        scalar offset_;
+
+
+    // Private Member Functions
+
+        //- Calculate the pow field and return true if successful
+        virtual bool calc();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("pow");
+
+
+    // Constructors
+
+        //- Construct from Time and dictionary
+        pow
+        (
+            const word& name,
+            const Time& runTime,
+            const dictionary& dict
+        );
+
+        //- No copy construct
+        pow(const pow&) = delete;
+
+        //- No copy assignment
+        void operator=(const pow&) = delete;
+
+
+    //- Destructor
+    virtual ~pow() = default;
+
+
+    // Member Functions
+
+        //- Read the randomise data
+        virtual bool read(const dictionary&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace functionObjects
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/distributionModels/RosinRammler/RosinRammler.C b/src/lagrangian/distributionModels/RosinRammler/RosinRammler.C
index 2db2244653cf90fd703e526748e59b086b376078..185d9c1570815a4664800e09b499a15b4c34487a 100644
--- a/src/lagrangian/distributionModels/RosinRammler/RosinRammler.C
+++ b/src/lagrangian/distributionModels/RosinRammler/RosinRammler.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2011-2020 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -77,10 +77,10 @@ Foam::distributionModels::RosinRammler::~RosinRammler()
 
 Foam::scalar Foam::distributionModels::RosinRammler::sample() const
 {
-    scalar K = 1.0 - exp(-pow((maxValue_ - minValue_)/d_, n_));
-    scalar y = rndGen_.sample01<scalar>();
-    scalar x = minValue_ + d_*::pow(-log(1.0 - y*K), 1.0/n_);
-    return x;
+    const scalar minValueByDPowN = pow(minValue_/d_, n_);
+    const scalar K = 1 - exp(- pow(maxValue_/d_, n_) + minValueByDPowN);
+    const scalar y = rndGen_.sample01<scalar>();
+    return d_*pow(minValueByDPowN - log(1 - K*y), 1/n_);
 }
 
 
diff --git a/src/lagrangian/distributionModels/RosinRammler/RosinRammler.H b/src/lagrangian/distributionModels/RosinRammler/RosinRammler.H
index 35c650987725ba76e4efabb64ef7bd6d61a795fc..ec4f94aa0fa0839e07a2d926464ff310a78701a6 100644
--- a/src/lagrangian/distributionModels/RosinRammler/RosinRammler.H
+++ b/src/lagrangian/distributionModels/RosinRammler/RosinRammler.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2011-2020 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -30,19 +30,18 @@ Description
     Rosin-Rammler distributionModel
 
    \f[
-       cumulative model =
-           (1.0 - exp( -(( x - d0)/d)^n )
-         / (1.0 - exp( -((d1 - d0)/d)^n )
+       CDF(x) =
+            (1 - exp(-(x/d)^n + (d_0/d)^n)
+           /(1 - exp(-(d_1/d)^n + (d_0/d)^n)
    \f]
 
-
 SourceFiles
     RosinRammler.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef distributionModels_RosinRammler_H
-#define distributionModels_RosinRammler_H
+#ifndef RosinRammler_H
+#define RosinRammler_H
 
 #include "distributionModel.H"
 
@@ -61,7 +60,7 @@ class RosinRammler
 :
     public distributionModel
 {
-    // Private data
+    // Private Data
 
         //- Distribution minimum
         scalar minValue_;
@@ -71,7 +70,10 @@ class RosinRammler
 
         // Model coefficients
 
+            //- Scale parameter
             scalar d_;
+
+            //- Shape parameter
             scalar n_;
 
 
diff --git a/src/lagrangian/spray/submodels/BreakupModel/SHF/SHF.C b/src/lagrangian/spray/submodels/BreakupModel/SHF/SHF.C
index 0cb8c00cbc5f936082ee6b174d13fd2081a7f880..a5e4524bd8a36e8612fc2ce27a8cc0406fc54afe 100644
--- a/src/lagrangian/spray/submodels/BreakupModel/SHF/SHF.C
+++ b/src/lagrangian/spray/submodels/BreakupModel/SHF/SHF.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011-2013 OpenFOAM Foundation
+    Copyright (C) 2011-2020 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -158,9 +158,21 @@ bool Foam::SHF<CloudType>::update
 
     scalar weGasCorr = weGas/(1.0 + weCorrCoeff_*ohnesorge);
 
-    // droplet deformation characteristic time
+    // update the droplet characteristic time
+    tc += dt;
 
-    scalar tChar = d/Urmag*sqrt(rho/rhoc);
+    // droplet deformation characteristic rate
+    scalar rChar = Urmag/d*sqrt(rhoc/rho);
+
+    // return if the characteristic deformation rate is too low for the
+    // following modelling to be calculable
+    if (tc*rChar < SMALL)
+    {
+        return false;
+    }
+
+    // droplet deformation characteristic time
+    scalar tChar = 1/rChar;
 
     scalar tFirst = cInit_*tChar;
 
@@ -173,9 +185,6 @@ bool Foam::SHF<CloudType>::update
     bool success = false;
 
 
-    // update the droplet characteristic time
-    tc += dt;
-
     if (weGas > weConst_)
     {
         if (weGas < weCrit1_)
diff --git a/src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/fixedMultiPhaseHeatFlux/fixedMultiPhaseHeatFluxFvPatchScalarField.C b/src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/fixedMultiPhaseHeatFlux/fixedMultiPhaseHeatFluxFvPatchScalarField.C
index 5ff761fbe42f7634bc921869f8323a79cbc15c6d..77722395ed019c770f92e7f841cc68ad045dd76d 100644
--- a/src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/fixedMultiPhaseHeatFlux/fixedMultiPhaseHeatFluxFvPatchScalarField.C
+++ b/src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/fixedMultiPhaseHeatFlux/fixedMultiPhaseHeatFluxFvPatchScalarField.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2015-2019 OpenFOAM Foundation
+    Copyright (C) 2015-2020 OpenFOAM Foundation
     Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -168,6 +168,31 @@ void Foam::fixedMultiPhaseHeatFluxFvPatchScalarField::updateCoeffs()
 }
 
 
+void Foam::fixedMultiPhaseHeatFluxFvPatchScalarField::autoMap
+(
+    const fvPatchFieldMapper& m
+)
+{
+    fixedValueFvPatchScalarField::autoMap(m);
+    m(q_);
+}
+
+
+void Foam::fixedMultiPhaseHeatFluxFvPatchScalarField::rmap
+(
+    const fvPatchScalarField& ptf,
+    const labelList& addr
+)
+{
+    fixedValueFvPatchScalarField::rmap(ptf, addr);
+
+    const fixedMultiPhaseHeatFluxFvPatchScalarField& mptf =
+        refCast<const fixedMultiPhaseHeatFluxFvPatchScalarField>(ptf);
+
+    q_.rmap(mptf.q_, addr);
+}
+
+
 void Foam::fixedMultiPhaseHeatFluxFvPatchScalarField::write(Ostream& os) const
 {
     fvPatchField<scalar>::write(os);
diff --git a/src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/fixedMultiPhaseHeatFlux/fixedMultiPhaseHeatFluxFvPatchScalarField.H b/src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/fixedMultiPhaseHeatFlux/fixedMultiPhaseHeatFluxFvPatchScalarField.H
index 9b059248765e80b0aa91499f155a4d7a4f417fbe..2d6e16eb5b58a8092951472362a9d04120ce3fe5 100644
--- a/src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/fixedMultiPhaseHeatFlux/fixedMultiPhaseHeatFluxFvPatchScalarField.H
+++ b/src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/fixedMultiPhaseHeatFlux/fixedMultiPhaseHeatFluxFvPatchScalarField.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2015-2018 OpenFOAM Foundation
+    Copyright (C) 2015-2020 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -144,6 +144,17 @@ public:
 
     // Member Functions
 
+        // Mapping functions
+
+            //- Map (and resize as needed) from self given a mapping object
+            //  Used to update fields following mesh topology change
+            virtual void autoMap(const fvPatchFieldMapper&);
+
+            //- Reverse map the given fvPatchField onto this fvPatchField
+            //  Used to reconstruct fields
+            virtual void rmap(const fvPatchScalarField&, const labelList&);
+
+
         // Evaluation Functions
 
             //- Update the coefficients associated with the patch field
diff --git a/src/sampling/meshToMesh0/calculateMeshToMesh0Addressing.C b/src/sampling/meshToMesh0/calculateMeshToMesh0Addressing.C
index e529a4e19621409a2a560acf2660c0aa56fe2143..6ca92bc00f001b7c24e7839faaa1e2e1c2b70843 100644
--- a/src/sampling/meshToMesh0/calculateMeshToMesh0Addressing.C
+++ b/src/sampling/meshToMesh0/calculateMeshToMesh0Addressing.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2011-2020 OpenFOAM Foundation
     Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -32,9 +32,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "meshToMesh0.H"
-#include "SubField.H"
 
-#include "indexedOctree.H"
 #include "treeDataCell.H"
 #include "treeDataFace.H"
 
@@ -273,6 +271,11 @@ void Foam::meshToMesh0::cellAddresses
             if (boundaryCell[curCell])
             {
                 cellAddressing_[toI] = oc.findInside(p);
+
+                if (cellAddressing_[toI] != -1)
+                {
+                    curCell = cellAddressing_[toI];
+                }
             }
             else
             {
@@ -325,6 +328,11 @@ void Foam::meshToMesh0::cellAddresses
                 {
                     // Still not found so use the octree
                     cellAddressing_[toI] = oc.findInside(p);
+
+                    if (cellAddressing_[toI] != -1)
+                    {
+                        curCell = cellAddressing_[toI];
+                    }
                 }
             }
         }
diff --git a/src/sampling/meshToMesh0/meshToMesh0.H b/src/sampling/meshToMesh0/meshToMesh0.H
index bb2ce9f784e9275632e7af9d83cd63766f8397f1..1d30a3fd63659b3e9afd8b45cfa2fa5ad48f6fbc 100644
--- a/src/sampling/meshToMesh0/meshToMesh0.H
+++ b/src/sampling/meshToMesh0/meshToMesh0.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2011-2020 OpenFOAM Foundation
     Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -28,7 +28,7 @@ Class
     Foam::meshToMesh0
 
 Description
-    mesh to mesh interpolation class.
+    Serial mesh to mesh interpolation class.
 
 Note
     This class is due to be deprecated in favour of meshToMesh0New
@@ -105,6 +105,9 @@ class meshToMesh0
 
     // Private Member Functions
 
+        //- Calculates mesh to mesh addressing pattern.
+        //  For each cell from one mesh find the closest cell centre
+        //  in the other mesh
         void calcAddressing();
 
         void cellAddresses
diff --git a/tutorials/verificationAndValidation/turbulentInflow/0.orig/U b/tutorials/verificationAndValidation/turbulentInflow/PCF/0.orig/U
similarity index 75%
rename from tutorials/verificationAndValidation/turbulentInflow/0.orig/U
rename to tutorials/verificationAndValidation/turbulentInflow/PCF/0.orig/U
index 07516edcc8002fb994c48722bdbc76ce5a6bc36b..42f41c742079d104bebca725228ef2d85a8e77e4 100644
--- a/tutorials/verificationAndValidation/turbulentInflow/0.orig/U
+++ b/tutorials/verificationAndValidation/turbulentInflow/PCF/0.orig/U
@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volVectorField;
-    location    "0";
     object      U;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -21,34 +20,28 @@ internalField   uniform (0 0 0);
 
 boundaryField
 {
-    bottomWall
+    "bottomWall|topWall"
     {
         type            fixedValue;
-        value           uniform (0 0 0);
+        value           $internalField;
     }
-    topWall
-    {
-        type            fixedValue;
-        value           uniform (0 0 0);
-    }
-    left
-    {
-        type            cyclic;
-    }
-    right
+
+    "left|right"
     {
         type            cyclic;
     }
+
     inlet
     {
-        value           uniform (0 0 0);
+        value           $internalField;
     }
     #include            "inlet/U"
+
     outlet
     {
         type            inletOutlet;
-        inletValue      uniform (0 0 0);
-        value           uniform (0 0 0);
+        inletValue      $internalField;
+        value           $internalField;
     }
 }
 
diff --git a/tutorials/verificationAndValidation/turbulentInflow/0.orig/inlet.digitalFilter/U b/tutorials/verificationAndValidation/turbulentInflow/PCF/0.orig/inlet.DFM/U
similarity index 72%
rename from tutorials/verificationAndValidation/turbulentInflow/0.orig/inlet.digitalFilter/U
rename to tutorials/verificationAndValidation/turbulentInflow/PCF/0.orig/inlet.DFM/U
index a4b0994dfa517202b7142c0b847a8d8f5a28dc60..85bf7d7c5375fcad43e903c03e300f5fd142083e 100644
--- a/tutorials/verificationAndValidation/turbulentInflow/0.orig/inlet.digitalFilter/U
+++ b/tutorials/verificationAndValidation/turbulentInflow/PCF/0.orig/inlet.DFM/U
@@ -16,12 +16,15 @@ FoamFile
 
 inlet
 {
-    type                turbulentDigitalFilterInlet;
-    variant             digitalFilter;
-    planeDivisions      ( 64 70 );
-    L                   ( 0.78035508 0.31085352 0.342261 0.1728125 0.171875
-                          0.22459375 0.172787596 0.171889998 0.224578995 );
-    patchNormalSpeed    20.133;
+    type    turbulentDigitalFilterInlet;
+    n       ( 64 70 );
+    L
+    (
+        0.78035508 0.31085352 0.342261 0.1728125 0.171875
+        0.22459375 0.172787596 0.171889998 0.224578995
+    );
+    Ubulk   20.133;
 }
 
+
 // ************************************************************************* //
diff --git a/tutorials/verificationAndValidation/turbulentInflow/0.orig/inlet.DFSEM/U b/tutorials/verificationAndValidation/turbulentInflow/PCF/0.orig/inlet.DFSEM/U
similarity index 99%
rename from tutorials/verificationAndValidation/turbulentInflow/0.orig/inlet.DFSEM/U
rename to tutorials/verificationAndValidation/turbulentInflow/PCF/0.orig/inlet.DFSEM/U
index 9de2182bc31e820f66eb8ed0b5b4e9db5d656396..e110001e29ff229a63e95a08442a0db7c700c3d4 100644
--- a/tutorials/verificationAndValidation/turbulentInflow/0.orig/inlet.DFSEM/U
+++ b/tutorials/verificationAndValidation/turbulentInflow/PCF/0.orig/inlet.DFSEM/U
@@ -22,4 +22,5 @@ inlet
     mapMethod       nearestCell;
 }
 
+
 // ************************************************************************* //
diff --git a/tutorials/verificationAndValidation/turbulentInflow/0.orig/inlet.reducedDigitalFilter/U b/tutorials/verificationAndValidation/turbulentInflow/PCF/0.orig/inlet.FSM/U
similarity index 72%
rename from tutorials/verificationAndValidation/turbulentInflow/0.orig/inlet.reducedDigitalFilter/U
rename to tutorials/verificationAndValidation/turbulentInflow/PCF/0.orig/inlet.FSM/U
index 49adf756c5d75d321a26c5118d5f7b10ca2bcc64..48fead67b76225a33bbc81e15f24421aff394ad4 100644
--- a/tutorials/verificationAndValidation/turbulentInflow/0.orig/inlet.reducedDigitalFilter/U
+++ b/tutorials/verificationAndValidation/turbulentInflow/PCF/0.orig/inlet.FSM/U
@@ -16,12 +16,16 @@ FoamFile
 
 inlet
 {
-    type                turbulentDigitalFilterInlet;
-    variant             reducedDigitalFilter;
-    planeDivisions      ( 64 70 );
-    L                   ( 0.78035508 0.31085352 0.342261 0.1728125 0.171875
-                          0.22459375 0.172787596 0.171889998 0.224578995 );
-    patchNormalSpeed    20.133;
+    type    turbulentDigitalFilterInlet;
+    fsm     true;
+    n       ( 64 70 );
+    L
+    (
+        0.78035508 0.31085352 0.342261 0.1728125 0.171875
+        0.22459375 0.172787596 0.171889998 0.224578995
+    );
+    Ubulk   20.133;
 }
 
+
 // ************************************************************************* //
diff --git a/tutorials/verificationAndValidation/turbulentInflow/0.orig/nut b/tutorials/verificationAndValidation/turbulentInflow/PCF/0.orig/nut
similarity index 80%
rename from tutorials/verificationAndValidation/turbulentInflow/0.orig/nut
rename to tutorials/verificationAndValidation/turbulentInflow/PCF/0.orig/nut
index 27ffbc17a8443b8111141f377b44c827ce5f8495..f914c43bbf5e7ede0f686e2bc4435841520d80da 100644
--- a/tutorials/verificationAndValidation/turbulentInflow/0.orig/nut
+++ b/tutorials/verificationAndValidation/turbulentInflow/PCF/0.orig/nut
@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    location    "0";
     object      nut;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -21,28 +20,17 @@ internalField   uniform 0;
 
 boundaryField
 {
-    bottomWall
+    "bottomWall|topWall"
     {
         type            zeroGradient;
     }
-    topWall
-    {
-        type            zeroGradient;
-    }
-    left
-    {
-        type            cyclic;
-    }
-    right
+
+    "left|right"
     {
         type            cyclic;
     }
-    inlet
-    {
-        type            calculated;
-        value           uniform 1e-08;
-    }
-    outlet
+
+    "inlet|outlet"
     {
         type            calculated;
         value           uniform 1e-08;
diff --git a/tutorials/verificationAndValidation/turbulentInflow/0.orig/p b/tutorials/verificationAndValidation/turbulentInflow/PCF/0.orig/p
similarity index 80%
rename from tutorials/verificationAndValidation/turbulentInflow/0.orig/p
rename to tutorials/verificationAndValidation/turbulentInflow/PCF/0.orig/p
index 5da3f8c16e10a321241fb25d0396d01a20ad6f7b..fe37c1a563e0fa1e67d7381a3d44749427900289 100644
--- a/tutorials/verificationAndValidation/turbulentInflow/0.orig/p
+++ b/tutorials/verificationAndValidation/turbulentInflow/PCF/0.orig/p
@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    location    "0";
     object      p;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -21,30 +20,20 @@ internalField   uniform 0;
 
 boundaryField
 {
-    bottomWall
+    "bottomWall|topWall|inlet"
     {
         type            zeroGradient;
     }
-    topWall
-    {
-        type            zeroGradient;
-    }
-    left
-    {
-        type            cyclic;
-    }
-    right
+
+    "left|right"
     {
         type            cyclic;
     }
-    inlet
-    {
-        type            zeroGradient;
-    }
+
     outlet
     {
         type            fixedValue;
-        value           uniform 0;
+        value           $internalField;
     }
 }
 
diff --git a/tutorials/verificationAndValidation/turbulentInflow/Allclean b/tutorials/verificationAndValidation/turbulentInflow/PCF/Allclean
similarity index 73%
rename from tutorials/verificationAndValidation/turbulentInflow/Allclean
rename to tutorials/verificationAndValidation/turbulentInflow/PCF/Allclean
index f66f87c3a8e4b5cb548e7c787db9a9bcb11574f9..3d330958d1a0d6bd3f29fa188906c953ce8af3f1 100755
--- a/tutorials/verificationAndValidation/turbulentInflow/Allclean
+++ b/tutorials/verificationAndValidation/turbulentInflow/PCF/Allclean
@@ -4,6 +4,11 @@ cd "${0%/*}" || exit                                # Run from this directory
 #------------------------------------------------------------------------------
 
 cleanCase0
-\rm -rf system/controlDict constant/boundaryData/inlet results
+
+rm -f system/controlDict
+rm -rf constant/boundaryData/inlet
+rm -rf results
+rm -f *.png
+rm -f constant/{R*,points*,UMean*}
 
 #------------------------------------------------------------------------------
diff --git a/tutorials/verificationAndValidation/turbulentInflow/Allrun b/tutorials/verificationAndValidation/turbulentInflow/PCF/Allrun
similarity index 55%
rename from tutorials/verificationAndValidation/turbulentInflow/Allrun
rename to tutorials/verificationAndValidation/turbulentInflow/PCF/Allrun
index 3a55dccafa02736629a9de02313de39745cbfb8d..4c57b41f0a35efef302fcd9de8beaf122b547b61 100755
--- a/tutorials/verificationAndValidation/turbulentInflow/Allrun
+++ b/tutorials/verificationAndValidation/turbulentInflow/PCF/Allrun
@@ -4,15 +4,6 @@ cd "${0%/*}" || exit                                # Run from this directory
 . ${WM_PROJECT_DIR:?}/bin/tools/CleanFunctions      # Tutorial clean functions
 #------------------------------------------------------------------------------
 
-endTime=10
-\cp system/controlDict.template system/controlDict
-if notTest "$@"
-then
-    endTime=85
-fi
-\sed -i "s|END_TIME|$endTime|g" system/controlDict
-
-
 # Collect data into the 'results' directory,
 # and clean the case for the next run
 #
@@ -20,14 +11,15 @@ fi
 # ----
 collectData(){
     model=$1
-    \echo "    Moving results into 'results/$model'"
-    results="results/$model"
-    \mkdir -p "$results"
+    runType=$2
+    echo "    Moving results into 'results/$model.$runType'"
+    results="results/$model.$runType"
+    mkdir -p "$results"
     timeDir=$(foamListTimes -latestTime)
-    \mv -f log* *.png postProcessing "$timeDir" "$results" 2>/dev/null
+    mv -f log* *.png postProcessing "$timeDir" "$results" 2>/dev/null
 
     cleanTimeDirectories
-    \rm -rf processor* > /dev/null 2>&1
+    rm -rf processor* > /dev/null 2>&1
 }
 
 
@@ -40,13 +32,13 @@ serialRun(){
     models=$*
     for model in $models
     do
-        \echo "    Running with the synthetic turbulence model: $model"
-        (\cd 0 && \ln -snf "inlet.$model" inlet)
-        (\cd constant/boundaryData && \ln -snf "inlet.$model" inlet)
+        echo "    Running with the synthetic turbulence model: $model"
+        (cd 0 && ln -snf "inlet.$model" inlet)
+        (cd constant/boundaryData && ln -snf "inlet.$model" inlet)
 
         runApplication -s "$model" $(getApplication)
         ./plot
-        collectData $model
+        collectData $model "serial"
     done
 }
 
@@ -60,36 +52,35 @@ parallelRun(){
     models=$*
     for model in $models
     do
-        \echo "    Running with the synthetic turbulence model: $model"
-        (\cd 0 && \ln -snf "inlet.$model" inlet)
-        (\cd constant/boundaryData && \ln -snf "inlet.$model" inlet)
+        echo "    Running with the synthetic turbulence model: $model"
+        (cd 0 && ln -snf "inlet.$model" inlet)
+        (cd constant/boundaryData && ln -snf "inlet.$model" inlet)
 
-        runApplication -s "$model" decomposePar
+        runApplication -s "$model" decomposePar -force
         runParallel -s "$model" $(getApplication)
+        runApplication -s "$model" reconstructPar -latestTime
         ./plot
-
-        collectData $model
+        collectData $model "parallel"
     done
 }
 
 
 #------------------------------------------------------------------------------
 
-# Synthetic inflow models
+# Prepare the numerical setup
+./Allrun.pre
+
+# Run with the synthetic turbulence models
 models="
-reducedDigitalFilter
-digitalFilter
+FSM
+DFM
 DFSEM
 "
 
-# Prepare the numerical setup
-runApplication blockMesh
-restore0Dir
-\rm -rf "results"
+parallelRun $models
 
-# Run with the synthetic turbulence models
 serialRun $models
-#parallelRun $models
 
 
 #------------------------------------------------------------------------------
+
diff --git a/tutorials/verificationAndValidation/turbulentInflow/PCF/Allrun.pre b/tutorials/verificationAndValidation/turbulentInflow/PCF/Allrun.pre
new file mode 100755
index 0000000000000000000000000000000000000000..2d3e50abf283808af53f87b5a86287a0f29c3505
--- /dev/null
+++ b/tutorials/verificationAndValidation/turbulentInflow/PCF/Allrun.pre
@@ -0,0 +1,20 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
+#------------------------------------------------------------------------------
+
+endTime=85
+if notTest "$@"
+then
+    endTime=10
+fi
+
+sed "s|END_TIME|$endTime|g" system/controlDict.template > system/controlDict
+
+restore0Dir
+
+runApplication blockMesh
+
+rm -rf results
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/verificationAndValidation/turbulentInflow/README b/tutorials/verificationAndValidation/turbulentInflow/PCF/README
similarity index 69%
rename from tutorials/verificationAndValidation/turbulentInflow/README
rename to tutorials/verificationAndValidation/turbulentInflow/PCF/README
index d17226d0b45458ef9c179b15bb7427349b8c6bf6..75345550869a95a8fb7acccd60d80448bd0e6929 100644
--- a/tutorials/verificationAndValidation/turbulentInflow/README
+++ b/tutorials/verificationAndValidation/turbulentInflow/PCF/README
@@ -1,12 +1,11 @@
-Synthetic turbulence inflow tests
-======================
+#------------------------------------------------------------------------------
 
 The following three synthetic turbulence inflow boundary conditions are
 examined through a single-cell-domain smooth-wall plane channel flow setup:
 
-- turbulentDFSEMInlet
-- turbulentDigitalFilterInlet variant=digitalFilter
-- turbulentDigitalFilterInlet variant=reducedDigitalFilter
+- turbulentDFSEMInlet (DFSEM)
+- turbulentDigitalFilterInlet (DFM)
+- turbulentDigitalFilterInlet with the forward-stepwise method (FSM)
 
 The input statistics are obtained from:
 
@@ -22,13 +21,12 @@ The data is available online from (Retrieved: 21-06-2019):
 
     https://turbulence.oden.utexas.edu/data/MKM/chan395/
 
-Serial executing (comment out 'parallelRun'):
-
-./Allrun
-
-Parallel (decompositionMethod=scotch) executing (comment out 'serialRun'):
+Executing:
 
 ./Allrun
 
 The script will run the test case, and collect the plots and samples into
 the 'results' directory.
+
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/verificationAndValidation/turbulentInflow/constant/boundaryData/inlet.DFSEM/0/R b/tutorials/verificationAndValidation/turbulentInflow/PCF/constant/boundaryData/inlet.DFM/0/R
similarity index 100%
rename from tutorials/verificationAndValidation/turbulentInflow/constant/boundaryData/inlet.DFSEM/0/R
rename to tutorials/verificationAndValidation/turbulentInflow/PCF/constant/boundaryData/inlet.DFM/0/R
diff --git a/tutorials/verificationAndValidation/turbulentInflow/constant/boundaryData/inlet.digitalFilter/0/UMean b/tutorials/verificationAndValidation/turbulentInflow/PCF/constant/boundaryData/inlet.DFM/0/UMean
similarity index 100%
rename from tutorials/verificationAndValidation/turbulentInflow/constant/boundaryData/inlet.digitalFilter/0/UMean
rename to tutorials/verificationAndValidation/turbulentInflow/PCF/constant/boundaryData/inlet.DFM/0/UMean
diff --git a/tutorials/verificationAndValidation/turbulentInflow/constant/boundaryData/inlet.DFSEM/points b/tutorials/verificationAndValidation/turbulentInflow/PCF/constant/boundaryData/inlet.DFM/points
similarity index 100%
rename from tutorials/verificationAndValidation/turbulentInflow/constant/boundaryData/inlet.DFSEM/points
rename to tutorials/verificationAndValidation/turbulentInflow/PCF/constant/boundaryData/inlet.DFM/points
diff --git a/tutorials/verificationAndValidation/turbulentInflow/constant/boundaryData/inlet.DFSEM/0/L b/tutorials/verificationAndValidation/turbulentInflow/PCF/constant/boundaryData/inlet.DFSEM/0/L
similarity index 100%
rename from tutorials/verificationAndValidation/turbulentInflow/constant/boundaryData/inlet.DFSEM/0/L
rename to tutorials/verificationAndValidation/turbulentInflow/PCF/constant/boundaryData/inlet.DFSEM/0/L
diff --git a/tutorials/verificationAndValidation/turbulentInflow/constant/boundaryData/inlet.digitalFilter/0/R b/tutorials/verificationAndValidation/turbulentInflow/PCF/constant/boundaryData/inlet.DFSEM/0/R
similarity index 100%
rename from tutorials/verificationAndValidation/turbulentInflow/constant/boundaryData/inlet.digitalFilter/0/R
rename to tutorials/verificationAndValidation/turbulentInflow/PCF/constant/boundaryData/inlet.DFSEM/0/R
diff --git a/tutorials/verificationAndValidation/turbulentInflow/constant/boundaryData/inlet.DFSEM/0/U b/tutorials/verificationAndValidation/turbulentInflow/PCF/constant/boundaryData/inlet.DFSEM/0/U
similarity index 100%
rename from tutorials/verificationAndValidation/turbulentInflow/constant/boundaryData/inlet.DFSEM/0/U
rename to tutorials/verificationAndValidation/turbulentInflow/PCF/constant/boundaryData/inlet.DFSEM/0/U
diff --git a/tutorials/verificationAndValidation/turbulentInflow/constant/boundaryData/inlet.digitalFilter/points b/tutorials/verificationAndValidation/turbulentInflow/PCF/constant/boundaryData/inlet.DFSEM/points
similarity index 100%
rename from tutorials/verificationAndValidation/turbulentInflow/constant/boundaryData/inlet.digitalFilter/points
rename to tutorials/verificationAndValidation/turbulentInflow/PCF/constant/boundaryData/inlet.DFSEM/points
diff --git a/tutorials/verificationAndValidation/turbulentInflow/constant/boundaryData/inlet.reducedDigitalFilter/0/R b/tutorials/verificationAndValidation/turbulentInflow/PCF/constant/boundaryData/inlet.FSM/0/R
similarity index 100%
rename from tutorials/verificationAndValidation/turbulentInflow/constant/boundaryData/inlet.reducedDigitalFilter/0/R
rename to tutorials/verificationAndValidation/turbulentInflow/PCF/constant/boundaryData/inlet.FSM/0/R
diff --git a/tutorials/verificationAndValidation/turbulentInflow/constant/boundaryData/inlet.reducedDigitalFilter/0/UMean b/tutorials/verificationAndValidation/turbulentInflow/PCF/constant/boundaryData/inlet.FSM/0/UMean
similarity index 100%
rename from tutorials/verificationAndValidation/turbulentInflow/constant/boundaryData/inlet.reducedDigitalFilter/0/UMean
rename to tutorials/verificationAndValidation/turbulentInflow/PCF/constant/boundaryData/inlet.FSM/0/UMean
diff --git a/tutorials/verificationAndValidation/turbulentInflow/constant/boundaryData/inlet.reducedDigitalFilter/points b/tutorials/verificationAndValidation/turbulentInflow/PCF/constant/boundaryData/inlet.FSM/points
similarity index 100%
rename from tutorials/verificationAndValidation/turbulentInflow/constant/boundaryData/inlet.reducedDigitalFilter/points
rename to tutorials/verificationAndValidation/turbulentInflow/PCF/constant/boundaryData/inlet.FSM/points
diff --git a/tutorials/verificationAndValidation/turbulentInflow/constant/transportProperties b/tutorials/verificationAndValidation/turbulentInflow/PCF/constant/transportProperties
similarity index 100%
rename from tutorials/verificationAndValidation/turbulentInflow/constant/transportProperties
rename to tutorials/verificationAndValidation/turbulentInflow/PCF/constant/transportProperties
diff --git a/tutorials/verificationAndValidation/turbulentInflow/constant/turbulenceProperties b/tutorials/verificationAndValidation/turbulentInflow/PCF/constant/turbulenceProperties
similarity index 100%
rename from tutorials/verificationAndValidation/turbulentInflow/constant/turbulenceProperties
rename to tutorials/verificationAndValidation/turbulentInflow/PCF/constant/turbulenceProperties
diff --git a/tutorials/verificationAndValidation/turbulentInflow/PCF/plot b/tutorials/verificationAndValidation/turbulentInflow/PCF/plot
new file mode 100755
index 0000000000000000000000000000000000000000..adf477bce334acb964cd78dba0e106de6a160c83
--- /dev/null
+++ b/tutorials/verificationAndValidation/turbulentInflow/PCF/plot
@@ -0,0 +1,135 @@
+#!/bin/sh
+cd "${0%/*}" || exit                        # Run from this directory
+#------------------------------------------------------------------------------
+
+plotCellR() {
+    timeDir=$1
+    echo "    Plotting the normal and Reynolds stresses on cell."
+
+    outName="stress-cell.png"
+    gnuplot<<PLT_CELL_R
+    set terminal pngcairo font "helvetica,20" size 1000, 800
+    set xrange [0:1]
+    set yrange [-1:8]
+    set grid
+    set key top right
+    set key samplen 2
+    set key spacing 0.75
+    set xlabel "Channel height from the bottomWall [m]"
+    set ylabel "<u_i^' u_i^'> [m2/s2]"
+    set offset .05, .05
+    set output "$outName"
+    set title "Normal and Reynolds stresses on cell"
+
+    input = "$timeDir/inletCell_UPrime2Mean.xy"
+    bench = "constant/pointsRdata"
+
+    plot \
+        input u 1:2 t "<u^' u^'>" w l lw 2 lc rgb "#009E73", \
+        input u 1:5 t "<v^' v^'>" w l lw 2 lc rgb "#F0E440", \
+        input u 1:7 t "<w^' w^'>" w l lw 2 lc rgb "#0072B2", \
+        input u 1:3 t "<u^' v^'>" w l lw 2 lc rgb "#D55E00", \
+        bench u 2:4 t "<u^' u^'>_{DNS}" w l lw 2 dt 2 lc rgb "#009E73", \
+        bench u 2:7 t "<v^' v^'>_{DNS}" w l lw 2 dt 2 lc rgb "#F0E440", \
+        bench u 2:9 t "<w^' w^'>_{DNS}" w l lw 2 dt 2 lc rgb "#0072B2", \
+        bench u 2:5 t "<u^' v^'>_{DNS}" w l lw 2 dt 2 lc rgb "#D55E00"
+PLT_CELL_R
+}
+
+
+plotPatchR() {
+    timeDir=$1
+    echo "    Plotting the normal and Reynolds stresses on inlet patch faces."
+
+    outName="stress-patch.png"
+    gnuplot<<PLT_PATCH_R
+    set terminal pngcairo font "helvetica,20" size 1000, 800
+    set xrange [0:1]
+    set yrange [-1:8]
+    set grid
+    set key top right
+    set key samplen 2
+    set key spacing 0.75
+    set xlabel "Channel height from the bottomWall [m]"
+    set ylabel "<u_i^' u_i^'> [m2/s2]"
+    set offset .05, .05
+    set output "$outName"
+    set title "Normal and Reynolds stresses on patch"
+
+    input = "$timeDir/inletPatch_UPrime2Mean.xy"
+    bench = "constant/pointsRdata"
+
+    plot \
+        input u 1:2 t "<u^' u^'>" w l lw 2 lc rgb "#009E73", \
+        input u 1:5 t "<v^' v^'>" w l lw 2 lc rgb "#F0E440", \
+        input u 1:7 t "<w^' w^'>" w l lw 2 lc rgb "#0072B2", \
+        input u 1:3 t "<u^' v^'>" w l lw 2 lc rgb "#D55E00", \
+        bench u 2:4 t "<u^' u^'>_{DNS}" w l lw 2 dt 2 lc rgb "#009E73", \
+        bench u 2:7 t "<v^' v^'>_{DNS}" w l lw 2 dt 2 lc rgb "#F0E440", \
+        bench u 2:9 t "<w^' w^'>_{DNS}" w l lw 2 dt 2 lc rgb "#0072B2", \
+        bench u 2:5 t "<u^' v^'>_{DNS}" w l lw 2 dt 2 lc rgb "#D55E00"
+PLT_PATCH_R
+}
+
+
+plotPatchUMean() {
+    timeDir=$1
+    echo "    Plotting the streamwise mean flow speed on inlet patch faces."
+
+    outName="u-patch.png"
+    gnuplot<<PLT_PATCH_UMEAN
+    set terminal pngcairo font "helvetica,20" size 1000, 800
+    set xrange [0:1]
+    set yrange [0:25]
+    set grid
+    set key top right
+    set key samplen 2
+    set key spacing 0.75
+    set xlabel "Channel height from the bottomWall [m]"
+    set ylabel "u [m/s]"
+    set offset .05, .05
+    set output "$outName"
+
+    input = "$timeDir/inletPatch_UMean.xy"
+    bench = "constant/pointsUMeanData"
+
+    plot \
+        input u 1:2 t "u" w l lw 2 lc rgb "#009E73", \
+        bench u 2:4 t "u_{DNS}" w l lw 2 dt 2 lc rgb "#009E73"
+PLT_PATCH_UMEAN
+}
+
+
+#------------------------------------------------------------------------------
+
+# Require gnuplot
+command -v gnuplot >/dev/null || {
+    echo "gnuplot not found - skipping graph creation" 1>&2
+    exit 1
+}
+
+# Prepare the benchmark data
+cp -f constant/boundaryData/inlet/0/R constant/R
+cp -f constant/boundaryData/inlet/points constant/points
+cp -f constant/boundaryData/inlet.DFM/0/UMean constant/UMean
+cat constant/R | tr -d '()' > constant/Rdata
+cat constant/points | tr -d '()' > constant/pointsData
+cat constant/UMean | tr -d '()' > constant/UMeanData
+paste constant/pointsData constant/Rdata > constant/pointsRdata
+paste constant/pointsData constant/UMeanData > constant/pointsUMeanData
+
+# The latestTime in postProcessing/sampling1
+timeDir=$(foamListTimes -case postProcessing/sampling1 -latestTime 2>/dev/null)
+[ -n "$timeDir" ] || {
+    echo "No postProcessing/sampling1 found - skipping graph creation" 1>&2
+    exit 1
+}
+
+timeDir="postProcessing/sampling1/$timeDir"
+
+plotCellR "$timeDir"
+plotPatchR "$timeDir"
+plotPatchUMean "$timeDir"
+
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/verificationAndValidation/turbulentInflow/system/blockMeshDict b/tutorials/verificationAndValidation/turbulentInflow/PCF/system/blockMeshDict
similarity index 99%
rename from tutorials/verificationAndValidation/turbulentInflow/system/blockMeshDict
rename to tutorials/verificationAndValidation/turbulentInflow/PCF/system/blockMeshDict
index fca41042cf8ad98e2d03cc6858a764f315fc3063..fc492a94bcd065ed5adea3d48424f1e76ac95bfe 100644
--- a/tutorials/verificationAndValidation/turbulentInflow/system/blockMeshDict
+++ b/tutorials/verificationAndValidation/turbulentInflow/PCF/system/blockMeshDict
@@ -88,4 +88,5 @@ mergePatchPairs
 (
 );
 
+
 // ************************************************************************* //
diff --git a/tutorials/verificationAndValidation/turbulentInflow/PCF/system/controlDict.template b/tutorials/verificationAndValidation/turbulentInflow/PCF/system/controlDict.template
new file mode 100644
index 0000000000000000000000000000000000000000..61b203660d9de898a60b69320e504d3b8132858e
--- /dev/null
+++ b/tutorials/verificationAndValidation/turbulentInflow/PCF/system/controlDict.template
@@ -0,0 +1,105 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     pimpleFoam;
+
+startFrom       latestTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         END_TIME;
+
+deltaT          4e-3;
+
+writeControl    timeStep;
+
+writeInterval   50;
+
+purgeWrite      3;
+
+writeFormat     ascii;
+
+writePrecision  8;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   8;
+
+runTimeModifiable false;
+
+adjustTimeStep  false;
+
+
+// Allow 10% of time for initialisation before sampling
+timeStart       #eval #{ 0.1 * ${/endTime} #};
+
+functions
+{
+    fieldAverage1
+    {
+        type                fieldAverage;
+        libs                (fieldFunctionObjects);
+        writeControl        writeTime;
+        timeStart           $/timeStart;
+
+        fields
+        (
+            U
+            {
+                mean        on;
+                prime2Mean  on;
+                base        time;
+            }
+        );
+    }
+
+    sampling1
+    {
+        type                sets;
+        libs                (sampling);
+        interpolationScheme cellPoint;
+        setFormat           raw;
+        writeControl        onEnd;
+        fields              (UPrime2Mean UMean);
+
+        sets
+        (
+            inletPatch
+            {
+                type        face;
+                axis        y;
+                start       (0.0 0 1.57);
+                end         (0.0 2 1.57);
+            }
+
+            inletCell
+            {
+                type        midPoint;
+                axis        y;
+                start       (0.062832 0 1.57);
+                end         (0.062832 2 1.57);
+            }
+        );
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/verificationAndValidation/turbulentInflow/system/decomposeParDict b/tutorials/verificationAndValidation/turbulentInflow/PCF/system/decomposeParDict
similarity index 89%
rename from tutorials/verificationAndValidation/turbulentInflow/system/decomposeParDict
rename to tutorials/verificationAndValidation/turbulentInflow/PCF/system/decomposeParDict
index f7450dd6845c7844d0675c1cbe943d3bae69f745..4d8edd3b6121fcc55bf1671c1861363901d0e386 100644
--- a/tutorials/verificationAndValidation/turbulentInflow/system/decomposeParDict
+++ b/tutorials/verificationAndValidation/turbulentInflow/PCF/system/decomposeParDict
@@ -10,13 +10,19 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    note        "mesh decomposition control dictionary";
     object      decomposeParDict;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-numberOfSubdomains  3;
+numberOfSubdomains  4;
+
+method              hierarchical;
+
+coeffs
+{
+    n    (1 2 2);
+}
 
-method          scotch;
 
 // ************************************************************************* //
+
diff --git a/tutorials/verificationAndValidation/turbulentInflow/system/fvSchemes b/tutorials/verificationAndValidation/turbulentInflow/PCF/system/fvSchemes
similarity index 99%
rename from tutorials/verificationAndValidation/turbulentInflow/system/fvSchemes
rename to tutorials/verificationAndValidation/turbulentInflow/PCF/system/fvSchemes
index 296000d302c1ddada09742cd9977e0fa1d625274..1828b72273d799be0c0dd57338a5750101ddc1f0 100644
--- a/tutorials/verificationAndValidation/turbulentInflow/system/fvSchemes
+++ b/tutorials/verificationAndValidation/turbulentInflow/PCF/system/fvSchemes
@@ -51,4 +51,5 @@ snGradSchemes
     default         corrected;
 }
 
+
 // ************************************************************************* //
diff --git a/tutorials/verificationAndValidation/turbulentInflow/system/fvSolution b/tutorials/verificationAndValidation/turbulentInflow/PCF/system/fvSolution
similarity index 100%
rename from tutorials/verificationAndValidation/turbulentInflow/system/fvSolution
rename to tutorials/verificationAndValidation/turbulentInflow/PCF/system/fvSolution
diff --git a/tutorials/verificationAndValidation/turbulentInflow/plot b/tutorials/verificationAndValidation/turbulentInflow/plot
deleted file mode 100755
index 3855ddd15cb6966f518cf0b50f643667ac4f99c0..0000000000000000000000000000000000000000
--- a/tutorials/verificationAndValidation/turbulentInflow/plot
+++ /dev/null
@@ -1,73 +0,0 @@
-#!/bin/sh
-cd "${0%/*}" || exit                        # Run from this directory
-#------------------------------------------------------------------------------
-
-plotStresses() {
-    timeDir=$1
-    \echo "    Plotting the normal and Reynolds stresses"
-
-    gnuplot<<PLT_STRESSES
-    set terminal pngcairo font "helvetica,20" size 1000, 800
-    set xrange [0:1]
-    set yrange [-1:8]
-    set grid
-    set key top right
-    set xlabel "Channel height from the bottomWall [m]"
-    set ylabel "<u_i^' u_i^'>"
-    set offset .05, .05
-
-    set style data linespoints
-    set linetype 1 lc rgb 'black' lw 2
-    set linetype 2 lc rgb 'red' lw 2
-    set linetype 3 lc rgb 'blue' lw 2
-    set linetype 4 lc rgb 'green' lw 2
-    set linetype 5 lc rgb 'black' pi -8 pt 4 ps 1.5
-    set linetype 6 lc rgb 'red' pi -8 pt 4 ps 1.5
-    set linetype 7 lc rgb 'blue' pi -8 pt 4 ps 1.5
-    set linetype 8 lc rgb 'green' pi -8 pt 4 ps 1.5
-
-    set title "Normal and Reynolds stresses on cell"
-    input="$timeDir/inletCell_UPrime2Mean.xy"
-    set output 'stress-cell.png'
-    plot \
-        input u 1:2 w lines t "<u^' u^'>" lt 1, \
-        input u 1:5 w lines t "<v^' v^'>" lt 2, \
-        input u 1:7 w lines t "<w^' w^'>" lt 3, \
-        input u 1:3 w lines t "<u^' v^'>" lt 4
-
-    set title "Normal and Reynolds stresses on patch"
-    input = "$timeDir/inletPatch_UPrime2Mean.xy"
-    set output 'stress-patch.png'
-    plot \
-        input u 1:2 w lines t "<u^' u^'>" lt 1, \
-        input u 1:5 w lines t "<v^' v^'>" lt 2, \
-        input u 1:7 w lines t "<w^' w^'>" lt 3, \
-        input u 1:3 w lines t "<u^' v^'>" lt 4
-
-PLT_STRESSES
-}
-
-
-#------------------------------------------------------------------------------
-
-# Require gnuplot
-command -v gnuplot >/dev/null || {
-    \echo "gnuplot not found - skipping graph creation" 1>&2
-    \exit 1
-}
-
-
-# The latestTime in postProcessing/inletSampling
-timeDir=$(foamListTimes -case postProcessing/inletSampling -latestTime 2>/dev/null)
-
-[ -n "$timeDir" ] || {
-    \echo "No postProcessing/inletSampling found - skipping graph creation" 1>&2
-    \exit 2
-}
-
-timeDir="postProcessing/inletSampling/$timeDir"
-
-plotStresses "$timeDir"
-
-
-#------------------------------------------------------------------------------