diff --git a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/turbulentDigitalFilterInletFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/turbulentDigitalFilterInletFvPatchVectorField.C
index bcbf9caa89..2b95c58c82 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
-)
+Foam::List<Foam::scalar> Foam::turbulentDigitalFilterInletFvPatchVectorField::
+calcCoeffs1FSM() const
 {
-    #ifdef FULLDEBUG
-    Info<< "Starts: computeDFM()" << nl;
-    #endif
+    List<scalar> coeffs1FSM(pTraits<vector>::nComponents);
 
-    if (Pstream::master())
+    forAll(L_.x(), i)
     {
-        embedTwoPointCorrs();
-        rndShiftRefill();
+        coeffs1FSM[i] = exp(C1FSM_/L_[i]);
     }
 
-    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
+    return coeffs1FSM;
 }
 
 
 Foam::List<Foam::scalar> Foam::turbulentDigitalFilterInletFvPatchVectorField::
-computeConstList1FSM() const
+calcCoeffs2FSM() const
 {
-    List<scalar> constList1FSM(pTraits<vector>::nComponents);
+    List<scalar> coeffs2FSM(pTraits<vector>::nComponents);
 
     forAll(L_.x(), i)
     {
-        constList1FSM[i] = exp(const1FSM_/L_[i]);
+        coeffs2FSM[i] = Foam::sqrt(1.0 - exp(C2FSM_/L_[i]));
     }
 
-    return constList1FSM;
-}
-
-
-Foam::List<Foam::scalar> Foam::turbulentDigitalFilterInletFvPatchVectorField::
-computeConstList2FSM() const
-{
-    List<scalar> constList2FSM(pTraits<vector>::nComponents);
-
-    forAll(L_.x(), i)
-    {
-        constList2FSM[i] = Foam::sqrt(1.0 - exp(const2FSM_/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_
+    Ubulk_(dict.getCheck<scalar>("Ubulk", scalarMinMax::ge(ROOTSMALL))),
+    C1_
     (
-        dict.getCheck<scalar>
-        (
-            "patchNormalSpeed",
-            scalarMinMax::ge(tiny_)
-        )
+        dict.getOrDefault<scalar>("C1", -0.5*constant::mathematical::pi)
     ),
-    modelConst_
+    perturb_
     (
-        dict.getOrDefault<scalar>
-        (
-            "modelConst",
-           -0.5*constant::mathematical::pi
-        )
+        dict.getCheckOrDefault<scalar>("perturb", 1e-5, scalarMinMax::ge(SMALL))
     ),
-    perturb_(dict.getOrDefault<scalar>("perturb", 1e-5)),
-    initialFlowRate_(pTraits<scalar>::one),
-    rndGen_
-    (
-        isFixedSeed_
-      ? 1234
-      : time(0)
-    ),
-    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 bbb0aa54f0..3fe83e2371 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)
+    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
-    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
+        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
 
-    - Forward-stepwise method-based generator (FSM)
-
-    \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
-    \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
+        // 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 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;
+        // 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.
+    - \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.
 
-SeeAlso
-    turbulentDFSEMInletFvPatchVectorField.C
+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 to enable the forward-stepwise method
+        const bool fsm_;
 
-        //- 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 select correlation function form: Gaussian or exponential
+        const bool Gaussian_;
 
-        //- Flag: random-number generator seed is fixed (default=true) or
-        //- generated pseudo-randomly based on clock-time per simulation
-        const bool isFixedSeed_;
+        //- 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: 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 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: mass flow rate is corrected on turbulence plane (default=true)
-        const bool isCorrectedFlowRate_;
+        //- Flag to correct mass flow rate on turbulence plane
+        const bool correctFlowRate_;
 
-        //- Flag: interpolate R field (default=false)
-        bool interpolateR_;
+        //- Internal flag to read R from data files
+        bool interpR_;
 
-        //- Flag: interpolate UMean field (default=false)
-        bool interpolateUMean_;
-
-        //- 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 patch-normal bulk flow speed [m/s]
+        const scalar Ubulk_;
 
-        //- Characteristic (e.g. bulk) mean speed of flow in the patch normal
-        //- direction [m/s]
-        const scalar patchNormalSpeed_;
+        //- Model constant shaping autocorrelation function (KSJ:Eq. 14)
+        const scalar C1_;
 
-        //- Model constant shaping autocorr function [-] (default='-0.5*pi')
-        //  (Klein et al., 2003, Eq. 14)
-        const scalar modelConst_;
-
-        //- 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();
+        //- shifting from the back to the front, and add new plane to the back
+        void shiftRefill();
 
-        //- Map two-point correlated random-number sets on patch based on chosen
-        //- mapping method
-        void mapFilteredRandomBox(vectorField& U);
+        //- Map two-point correlated random-number sets
+        //- on patch based on chosen mapping method
+        void mapFilteredBox(vectorField& U);
 
-        //- Map R_ on patch
-        void embedOnePointCorrs(vectorField& U) const;
+        //- Embed one-point correlations, i.e. R, on patch
+        void onePointCorrs(vectorField& U) const;
 
-        //- Map UMean_ on patch
-        void embedMeanVelocity(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();
 
-        //- Correct mass/vol flow rate in (only) streamwise direction
-        void correctFlowRate(vectorField& U) const;
+        //- Compute coeffs1FSM_ once per simulation
+        List<scalar> calcCoeffs1FSM() const;
 
-        //- 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 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 03a17b787f..706d7cd720 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/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 07516edcc8..42f41c7420 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 a4b0994dfa..85bf7d7c53 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 9de2182bc3..e110001e29 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 49adf756c5..48fead67b7 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 27ffbc17a8..f914c43bbf 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
+
+    "left|right"
     {
         type            cyclic;
     }
-    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 5da3f8c16e..fe37c1a563 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
+
+    "left|right"
     {
         type            cyclic;
     }
-    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 f66f87c3a8..3d330958d1 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 3a55dccafa..4c57b41f0a 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 0000000000..2d3e50abf2
--- /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 d17226d0b4..7534555086 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 0000000000..adf477bce3
--- /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 fca41042cf..fc492a94bc 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/system/controlDict.template b/tutorials/verificationAndValidation/turbulentInflow/PCF/system/controlDict.template
similarity index 54%
rename from tutorials/verificationAndValidation/turbulentInflow/system/controlDict.template
rename to tutorials/verificationAndValidation/turbulentInflow/PCF/system/controlDict.template
index 91209ef2c3..61b203660d 100644
--- a/tutorials/verificationAndValidation/turbulentInflow/system/controlDict.template
+++ b/tutorials/verificationAndValidation/turbulentInflow/PCF/system/controlDict.template
@@ -29,7 +29,7 @@ deltaT          4e-3;
 
 writeControl    timeStep;
 
-writeInterval   1250;
+writeInterval   50;
 
 purgeWrite      3;
 
@@ -47,13 +47,59 @@ runTimeModifiable false;
 
 adjustTimeStep  false;
 
+
 // Allow 10% of time for initialisation before sampling
 timeStart       #eval #{ 0.1 * ${/endTime} #};
 
 functions
 {
-    #include "fieldAverage"
-    #include "sampling"
+    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 f7450dd684..4d8edd3b61 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 296000d302..1828b72273 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 3855ddd15c..0000000000
--- 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"
-
-
-#------------------------------------------------------------------------------
diff --git a/tutorials/verificationAndValidation/turbulentInflow/system/fieldAverage b/tutorials/verificationAndValidation/turbulentInflow/system/fieldAverage
deleted file mode 100644
index 2af894e391..0000000000
--- a/tutorials/verificationAndValidation/turbulentInflow/system/fieldAverage
+++ /dev/null
@@ -1,36 +0,0 @@
-/*--------------------------------*- 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      fieldAverage;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-fieldAverage1
-{
-    type                fieldAverage;
-    libs                (fieldFunctionObjects);
-    writeControl        writeTime;
-    timeStart           $/timeStart;
-
-    fields
-    (
-        U
-        {
-            mean        on;
-            prime2Mean  on;
-            base        time;
-        }
-    );
-}
-
-// ************************************************************************* //
diff --git a/tutorials/verificationAndValidation/turbulentInflow/system/sampling b/tutorials/verificationAndValidation/turbulentInflow/system/sampling
deleted file mode 100644
index ef34447e5a..0000000000
--- a/tutorials/verificationAndValidation/turbulentInflow/system/sampling
+++ /dev/null
@@ -1,48 +0,0 @@
-/*--------------------------------*- 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      sampling;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-inletSampling
-{
-    type                sets;
-    libs                (sampling);
-    writeControl        writeTime;
-    timeStart           $/timeStart;
-
-    interpolationScheme cellPoint;
-    setFormat           raw;
-    fields              (UPrime2Mean);
-
-    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);
-        }
-    );
-}
-
-// ************************************************************************* //