From 33ef139ac1ac539a556626e13fc49719866f984d Mon Sep 17 00:00:00 2001
From: Kutalmis Bercin <kutalmis.bercin@esi-group.com>
Date: Mon, 18 Feb 2019 20:12:48 +0000
Subject: [PATCH] ENH: Digital-Filter Based Synthetic Turbulence Generation
 Method for LES/DES Inflows

    Velocity boundary condition generating synthetic turbulence-alike
    time-series for LES and DES turbulent flow computations.

    To this end, two synthetic turbulence generators can be chosen:
    - Digital-filter method-based generator (DFM)

    \verbatim
    Klein, M., Sadiki, A., and Janicka, J.
        A digital filter based generation of inflow data for spatially
        developing direct numerical or large eddy simulations,
        Journal of Computational Physics (2003) 186(2):652-665.
        doi:10.1016/S0021-9991(03)00090-1
    \endverbatim

    - Forward-stepwise method-based generator (FSM)

    \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
    of target statistics (mostly mean flow, Reynolds stress tensor profiles and
    length-scale sets) are fused 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.

    Random number sets ---->-|
                             |
                         DFM or FSM ---> New stochastic time-series consisting
                             |           turbulence statistics
    Turbulence statistics ->-|

    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.
---
 src/finiteVolume/Make/files                   |    1 +
 ...lentDigitalFilterInletFvPatchVectorField.C | 1261 +++++++++++++++++
 ...lentDigitalFilterInletFvPatchVectorField.H |  569 ++++++++
 ...alFilterInletFvPatchVectorFieldTemplates.C |  127 ++
 4 files changed, 1958 insertions(+)
 create mode 100644 src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/turbulentDigitalFilterInletFvPatchVectorField.C
 create mode 100644 src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/turbulentDigitalFilterInletFvPatchVectorField.H
 create mode 100644 src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/turbulentDigitalFilterInletFvPatchVectorFieldTemplates.C

diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files
index 791ef91c32..81500b3e13 100644
--- a/src/finiteVolume/Make/files
+++ b/src/finiteVolume/Make/files
@@ -200,6 +200,7 @@ $(derivedFvPatchFields)/translatingWallVelocity/translatingWallVelocityFvPatchVe
 $(derivedFvPatchFields)/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C
 $(derivedFvPatchFields)/turbulentDFSEMInlet/eddy/eddy.C
 $(derivedFvPatchFields)/turbulentDFSEMInlet/eddy/eddyIO.C
+$(derivedFvPatchFields)/turbulentDigitalFilterInlet/turbulentDigitalFilterInletFvPatchVectorField.C
 $(derivedFvPatchFields)/turbulentInlet/turbulentInletFvPatchFields.C
 $(derivedFvPatchFields)/turbulentIntensityKineticEnergyInlet/turbulentIntensityKineticEnergyInletFvPatchScalarField.C
 $(derivedFvPatchFields)/uniformDensityHydrostaticPressure/uniformDensityHydrostaticPressureFvPatchScalarField.C
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/turbulentDigitalFilterInletFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/turbulentDigitalFilterInletFvPatchVectorField.C
new file mode 100644
index 0000000000..c1a2b841ab
--- /dev/null
+++ b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/turbulentDigitalFilterInletFvPatchVectorField.C
@@ -0,0 +1,1261 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "turbulentDigitalFilterInletFvPatchVectorField.H"
+#include "volFields.H"
+#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&
+Foam::turbulentDigitalFilterInletFvPatchVectorField::patchMapper() const
+{
+    // Initialise interpolation (2D planar interpolation by triangulation)
+    if (mapperPtr_.empty())
+    {
+        // Reread values and interpolate
+        fileName samplePointsFile
+        (
+            this->db().time().path()
+           /this->db().time().caseConstant()
+           /"boundaryData"
+           /this->patch().name()
+           /"points"
+        );
+
+        pointField samplePoints((IFstream(samplePointsFile)()));
+
+        // tbd: run-time selection
+        bool nearestOnly =
+        (
+           !mapMethod_.empty()
+         && mapMethod_ != "planarInterpolation"
+        );
+
+        // Allocate the interpolator
+        mapperPtr_.reset
+        (
+            new pointToPointPlanarInterpolation
+            (
+                samplePoints,
+                this->patch().patch().faceCentres(),
+                perturb_,
+                nearestOnly
+            )
+        );
+    }
+
+    return *mapperPtr_;
+}
+
+
+Foam::List<Foam::Pair<Foam::label>>
+Foam::turbulentDigitalFilterInletFvPatchVectorField::patchIndexPairs()
+{
+    // Get patch normal direction into the domain
+    const vector nf(computePatchNormal());
+
+    // Find the second local coordinate direction
+    direction minCmpt = -1;
+    scalar minMag = VGREAT;
+    for (direction cmpt = 0; cmpt < pTraits<vector>::nComponents; ++cmpt)
+    {
+        scalar s = mag(nf[cmpt]);
+        if (s < minMag)
+        {
+            minMag = s;
+            minCmpt = cmpt;
+        }
+    }
+
+    // Create the second local coordinate direction
+    vector e2(Zero);
+    e2[minCmpt] = 1.0;
+
+    // Remove normal component
+    e2 -= (nf&e2)*nf;
+
+    // Create the local coordinate system
+    coordSystem::cartesian cs
+    (
+        Zero,   // origin
+        nf,     // normal
+        e2      // 0-axis
+    );
+
+    // Convert patch points into local coordinate system
+    const pointField localPos
+    (
+        cs.localPosition
+        (
+            pointField
+            (
+                patch().patch().points(),
+                patch().patch().meshPoints()
+            )
+        )
+    );
+
+    const boundBox bb(localPos);
+
+    // Normalise to (i, j) coordinates
+    const Vector2D<label> n(planeDivisions_.first(), planeDivisions_.second());
+    invDelta_[0] = n[0]/bb.span()[0];
+    invDelta_[1] = n[1]/bb.span()[1];
+    const point localMinPt(bb.min());
+
+    // Compute virtual-actual patch index pairs
+    List<Pair<label>> indexPairs(this->size(), Pair<label>(Zero, Zero));
+
+    // Virtual turbulence plane indices
+    label j = 0;
+    label k = 0;
+
+    forAll(*this, facei)
+    {
+        const scalar& centre0 = localPos[facei][0];
+        const scalar& centre1 = localPos[facei][1];
+
+        j = label((centre0 - localMinPt[0])*invDelta_[0]);
+        k = label((centre1 - localMinPt[1])*invDelta_[1]);
+
+        indexPairs[facei] = Pair<label>(facei, k*n[0] + j);
+    }
+
+    return indexPairs;
+}
+
+
+void Foam::turbulentDigitalFilterInletFvPatchVectorField::
+checkRTensorRealisable() const
+{
+    const vectorField& faceCentres = this->patch().patch().faceCentres();
+
+    forAll(R_, facei)
+    {
+        if (R_[facei].xx() <= 0)
+        {
+            FatalErrorInFunction
+                << "Reynolds stress tensor component Rxx cannot be negative"
+                   "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);
+        }
+
+        scalar term0 = R_[facei].xx()*R_[facei].yy() - sqr(R_[facei].xy());
+
+        if (term0 <= 0)
+        {
+            FatalErrorInFunction
+                << "Reynolds stress tensor component group, Rxx*Ryy - Rxy^2"
+                << "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;
+
+        if (term3 < 0)
+        {
+            FatalErrorInFunction
+                << "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);
+        }
+    }
+
+    #ifdef FULLDEBUG
+    Info<< "Ends: checkRTensorRealisable()."
+        << " Reynolds tensor (on patch) is consistent." << nl;
+    #endif
+}
+
+
+Foam::symmTensorField Foam::turbulentDigitalFilterInletFvPatchVectorField::
+computeLundWuSquires() const
+{
+    checkRTensorRealisable();
+
+    symmTensorField LundWuSquires(symmTensorField(R_.size()));
+
+    forAll(R_, facei)
+    {
+        const symmTensor& R = R_[facei];
+        symmTensor& lws = LundWuSquires[facei];
+
+        // (Klein et al., 2003, Eq. 5)
+        lws.xx() = Foam::sqrt(R.xx());
+        lws.xy() = R.xy()/lws.xx();
+        lws.xz() = R.xz()/lws.xx();
+        lws.yy() = Foam::sqrt(R.yy() - sqr(lws.xy()));
+        lws.yz() = (R.yz() - lws.xy()*lws.xz())/lws.yy();
+        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();
+}
+
+
+Foam::scalar Foam::turbulentDigitalFilterInletFvPatchVectorField::
+computeInitialFlowRate() const
+{
+    const vector patchNormal(computePatchNormal());
+    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
+(
+    const tensor& L
+) const
+{
+    tensor Ls(L);
+    convertToTimeScale(Ls);
+    const scalar invDeltaT = 1.0/db().time().deltaTValue();
+
+    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
+{
+    label initValue = 0;
+    label rangeModifier = 0;
+
+    if (variant_ == variantType::FORWARD_STEPWISE)
+    {
+        // Initialise with 1 since x-dir possess 1 node with this variant
+        initValue = pTraits<label>::nComponents;
+        rangeModifier = pTraits<vector>::nComponents;
+    }
+
+    List<label> lenRandomBox(pTraits<tensor>::nComponents, initValue);
+    Vector<label> lenGrid
+    (
+        pTraits<label>::nComponents,
+        planeDivisions_.first(),
+        planeDivisions_.second()
+    );
+
+    // Initialise: Main convenience factor, lenRandomBox_
+    for
+    (
+        const label& i
+      : labelRange(rangeModifier, pTraits<tensor>::nComponents - rangeModifier)
+    )
+    {
+        // Slicing index
+        const label sliceI = label(i/pTraits<vector>::nComponents);
+
+        // Refer to 'computeFilterCoeffs()'
+        const label n = ceil(L_[i]);
+        const label twiceN = 4*n;       
+
+        // Initialise: Random-number set sizes
+        lenRandomBox[i] = lenGrid[sliceI] + twiceN;
+    }
+
+    return lenRandomBox;
+}
+
+
+Foam::List<Foam::label> Foam::turbulentDigitalFilterInletFvPatchVectorField::
+initBoxFactors2D() 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];
+    }
+
+    return boxFactors2D;
+}
+
+
+Foam::List<Foam::label> Foam::turbulentDigitalFilterInletFvPatchVectorField::
+initBoxFactors3D() const
+{
+    List<label> boxFactors3D(pTraits<vector>::nComponents);
+
+    for (direction dir = 0; dir < pTraits<vector>::nComponents; ++dir)
+    {
+        boxFactors3D[dir] = randomBoxFactors2D_[dir]*lenRandomBox_[dir];
+    }
+
+    return boxFactors3D;
+}
+
+
+Foam::List<Foam::label> Foam::turbulentDigitalFilterInletFvPatchVectorField::
+initBoxPlaneFactors() 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);
+    }
+
+    return planeFactors;
+}
+
+
+Foam::List<Foam::List<Foam::scalar>>
+Foam::turbulentDigitalFilterInletFvPatchVectorField::fillRandomBox()
+{
+    List<List<scalar>> randomBox(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]);
+    }
+
+    #ifdef FULLDEBUG
+    Info<< "Ends: fillRandomBox()."
+        << "Random-number box is created." << nl;
+    #endif
+
+    return randomBox;
+}
+
+
+Foam::List<Foam::List<Foam::scalar>>
+Foam::turbulentDigitalFilterInletFvPatchVectorField::computeFilterCoeffs() const
+{
+    List<List<scalar>> filterCoeffs
+    (
+        pTraits<tensor>::nComponents,
+        List<scalar>(1, 1.0)
+    );
+
+    label initValue = 0;
+
+    if(variant_ == variantType::FORWARD_STEPWISE)
+    {
+        initValue = pTraits<vector>::nComponents;
+    }
+
+    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)
+        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)
+        const label twiceN = 4*n;
+
+        // Initialise filter-coeffs containers with full filter-width size
+        // Extra elem is mid-zero within [-N, N]
+        filterCoeffs[dir] = List<scalar>(twiceN + 1, Zero);
+
+        // 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)
+        std::iota
+        (
+            filterCoeffs[dir].begin(),
+            filterCoeffs[dir].end(),
+            initElem
+        );
+
+        // Compute filter-coeffs
+        // (Klein et al., 2003, Eq. 14 (Gaussian))
+        // (Bercin et al., 2018, Fig. 9 (Exponential))
+        List<scalar> fTemp(filterCoeffs[dir]);
+        scalar fSum = 0;
+        const scalar nSqr = n*n;
+
+        if (isGaussian_)
+        {
+            fTemp = sqr(exp(modelConst_*sqr(fTemp)/nSqr));
+            fSum = Foam::sqrt(sum(fTemp));
+
+            filterCoeffs[dir] =
+                exp(modelConst_*sqr(filterCoeffs[dir])/nSqr)/fSum;
+        }
+        else
+        {
+            fTemp = sqr(exp(modelConst_*mag(fTemp)/n));
+            fSum = Foam::sqrt(sum(fTemp));
+
+            filterCoeffs[dir] = exp(modelConst_*mag(filterCoeffs[dir])/n)/fSum;
+        }
+    }
+
+    #ifdef FULLDEBUG
+    Info<< "Ends: computeFilterCoeffs()."
+        << " Filter coefficients are computed." << nl;
+    #endif
+
+    return filterCoeffs;
+}
+
+
+void Foam::turbulentDigitalFilterInletFvPatchVectorField::rndShiftRefill()
+{
+    forAll(randomBox_, dir)
+    {
+        List<scalar>& r = randomBox_[dir];
+
+        // Shift forward from the back to the front / First Out
+        inplaceRotateList(r, randomBoxFactors2D_[dir]);
+
+        // Refill the back with a new random-number set / First In
+        for (label i = 0; i < randomBoxFactors2D_[dir]; ++i)
+        {
+            r[i] = rndGen_.GaussNormal<scalar>();
+        }
+    }
+}
+
+
+void Foam::turbulentDigitalFilterInletFvPatchVectorField::mapFilteredRandomBox
+(
+    vectorField& U
+)
+{
+    for (const auto& x : indexPairs_)
+    {
+        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];
+    }
+}
+
+
+void Foam::turbulentDigitalFilterInletFvPatchVectorField::embedOnePointCorrs
+(
+    vectorField& U
+) const
+{
+    forAll(LundWuSquires_, facei)
+    {
+        vector& Us = U[facei];
+        const symmTensor& lws = LundWuSquires_[facei];
+
+        // (Klein et al. 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();
+    }
+}
+
+
+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()
+{
+    for (direction dir = 0; dir < pTraits<vector>::nComponents; ++dir)
+    {
+        List<scalar>& in = randomBox_[dir];
+        List<scalar>& out = filteredRandomBox_[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 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 validSlice2 = sz2 - (szfilter2 - label(1));
+        const label validSlice3 = sz3 - (szfilter3 - label(1));
+
+        // Convolution summation - Along 1st direction
+        scalarField tmp(sz123);
+        label filterCentre = label(szfilter2/label(2));
+        label endIndex = sz2 - filterCentre;
+        label i0 = 0;
+        label i1 = 0;
+        label i2 = 0;
+
+        for (label i = 0; i < sz1; ++i)
+        {
+            for (label j = 0; j < sz3; ++j)
+            {
+                i1 += filterCentre;
+
+                for (label k = filterCentre; k < endIndex; ++k)
+                {
+                    tmp[i1] = 0.0;
+                    label q = 0;
+
+                    for (label p = szfilter2 - 1; p >= 0; --p, ++q)
+                    {
+                        tmp[i1] += in[i0 + q]*filter2[p];
+                    }
+                    ++i0;
+                    ++i1;
+
+                }
+                i0 += (filterCentre + filterCentre);
+                i1 += filterCentre;
+            }
+        }
+
+        // Convolution summation - Along 2nd direction
+        scalarField tmp2(tmp);
+        filterCentre = label(szfilter3/label(2));
+        endIndex = sz3 - filterCentre;
+        i1 = 0;
+        i2 = 0;
+
+        for (label i = 0; i < sz1; ++i)
+        {
+           label sl = i*sz23;
+
+            for (label j = 0; j < sz2; ++j)
+            {
+                i1 = j + sl;
+                i2 = i1;
+
+                for (label k = 0; k < endIndex - filterCentre; ++k)
+                {
+                    tmp[i1] = 0.0;
+                    label q = 0;
+
+                    for (label p = szfilter3 - 1; p >= 0; --p, ++q)
+                    {
+                        tmp[i1] += tmp2[i2 + q*sz2]*filter3[p];
+                    }
+                    i1 += sz2;
+                    i2 += sz2;
+
+                }
+                i1 += (sz2 + filterCentre);
+                i2 += (sz2 + filterCentre);
+            }
+        }
+
+        // Convolution summation - Along 3rd direction
+        filterCentre = label(szfilter1/label(2));
+        endIndex = sz1 - filterCentre;
+        i1 = (szfilter2 - label(1))/label(2);
+        i2 = (szfilter2 - label(1))/label(2);
+        label i3 = 0;
+
+        for (label i = 0; i < validSlice3; ++i)
+        {
+            for (label j = 0; j < validSlice2; ++j)
+            {
+                scalar sum = 0.0;
+                i1 = i2 + j;
+
+                for (label k = szfilter1 - 1; k >= 0; --k)
+                {
+                    sum += tmp[i1]*filter1[k];
+                    i1 += sz23;
+                }
+                out[i3] = sum;
+                ++i3;
+
+            }
+            i2 += sz2;
+        }
+    }
+}
+
+
+void Foam::turbulentDigitalFilterInletFvPatchVectorField::computeDFM
+(
+    vectorField& U
+)
+{
+    #ifdef FULLDEBUG
+    Info<< "Starts: computeDFM()" << nl;
+    #endif
+
+    if (Pstream::master())
+    {
+        embedTwoPointCorrs();
+        rndShiftRefill();
+    }
+
+    Pstream::scatter(filteredRandomBox_);
+
+    mapFilteredRandomBox(U);
+
+    embedOnePointCorrs(U);
+
+    embedMeanVelocity(U);
+
+    if (isCorrectedFlowRate_)
+    {
+        correctFlowRate(U);
+    }
+
+    #ifdef FULLDEBUG
+    Info<< "Ends: computeDFM()" << nl;
+    #endif
+}
+
+
+void Foam::turbulentDigitalFilterInletFvPatchVectorField::computeReducedDFM
+(
+    vectorField& U
+)
+{
+    #ifdef FULLDEBUG
+    Info<< "Starts: computeReducedDFM()" << nl;
+    #endif
+
+    if (Pstream::master())
+    {
+        embedTwoPointCorrs();
+        rndShiftRefill();
+    }
+
+    Pstream::scatter(filteredRandomBox_);
+
+    mapFilteredRandomBox(U);
+
+    computeFSM(U);
+
+    embedOnePointCorrs(U);
+
+    embedMeanVelocity(U);
+
+    if (isCorrectedFlowRate_)
+    {
+        correctFlowRate(U);
+    }
+
+    #ifdef FULLDEBUG
+    Info<< "Ends: computeReducedDFM()" << nl;
+    #endif
+}
+
+
+Foam::List<Foam::scalar> Foam::turbulentDigitalFilterInletFvPatchVectorField::
+computeConstList1FSM() const
+{
+    List<scalar> constList1FSM(pTraits<vector>::nComponents);
+
+    forAll(L_.x(), i)
+    {
+        constList1FSM[i] = exp(const1FSM_/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;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::turbulentDigitalFilterInletFvPatchVectorField::
+turbulentDigitalFilterInletFvPatchVectorField
+(
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF
+)
+:
+    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),
+    mapMethod_("nearestCell"),
+    curTimeIndex_(-1),
+    tiny_(1e-8),
+    patchNormalSpeed_(Zero),
+    modelConst_(-0.5*constant::mathematical::pi),
+    perturb_(1e-5),
+    initialFlowRate_(pTraits<scalar>::one),
+    rndGen_(1234),
+    planeDivisions_(Zero, Zero),
+    invDelta_(Zero),
+    indexPairs_(Zero),
+    R_(Zero),
+    LundWuSquires_(Zero),
+    UMean_(Zero),
+    Lbak_(Zero),
+    L_(Zero),
+    const1FSM_(Zero),
+    const2FSM_(Zero),
+    constList1FSM_(Zero),
+    constList2FSM_(Zero),
+    lenRandomBox_(Zero),
+    randomBoxFactors2D_(Zero),
+    randomBoxFactors3D_(Zero),
+    iNextToLastPlane_(Zero),
+    randomBox_(Zero),
+    filterCoeffs_(Zero),
+    filteredRandomBox_(Zero),
+    U0_(Zero),
+    computeVariant(nullptr)
+{}
+
+
+Foam::turbulentDigitalFilterInletFvPatchVectorField::
+turbulentDigitalFilterInletFvPatchVectorField
+(
+    const turbulentDigitalFilterInletFvPatchVectorField& ptf,
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    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_),
+    mapMethod_(ptf.mapMethod_),
+    curTimeIndex_(-1),
+    tiny_(ptf.tiny_),
+    patchNormalSpeed_(ptf.patchNormalSpeed_),
+    modelConst_(ptf.modelConst_),
+    perturb_(ptf.perturb_),
+    initialFlowRate_(ptf.initialFlowRate_),
+    rndGen_(ptf.rndGen_),
+    planeDivisions_(ptf.planeDivisions_),
+    invDelta_(ptf.invDelta_),
+    indexPairs_(ptf.indexPairs_),
+    R_(ptf.R_),
+    LundWuSquires_(ptf.LundWuSquires_),
+    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_),
+    iNextToLastPlane_(ptf.iNextToLastPlane_),
+    randomBox_(ptf.randomBox_),
+    filterCoeffs_(ptf.filterCoeffs_),
+    filteredRandomBox_(ptf.filteredRandomBox_),
+    U0_(ptf.U0_),
+    computeVariant(ptf.computeVariant)
+{}
+
+
+Foam::turbulentDigitalFilterInletFvPatchVectorField::
+turbulentDigitalFilterInletFvPatchVectorField
+(
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    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)),
+    mapMethod_(dict.getOrDefault<word>("mapMethod", "nearestCell")),
+    curTimeIndex_(-1),
+    tiny_(dict.getOrDefault<scalar>("threshold", 1e-8)),
+    patchNormalSpeed_
+    (
+        dict.getCheck<scalar>
+        (
+            "patchNormalSpeed",
+            scalarMinMax::ge(tiny_)
+        )
+    ),
+    modelConst_
+    (
+        dict.getOrDefault<scalar>
+        (
+            "modelConst",
+           -0.5*constant::mathematical::pi
+        )
+    ),
+    perturb_(dict.getOrDefault<scalar>("perturb", 1e-5)),
+    initialFlowRate_(pTraits<scalar>::one),
+    rndGen_
+    (
+        isFixedSeed_
+      ? 1234
+      : time(0)
+    ),
+    planeDivisions_
+    (
+        dict.getCheck<Tuple2<label, label>>
+        (
+            "planeDivisions",
+            [&](const Tuple2<label, label>& len)
+            {
+                return tiny_ < min(len.first(), len.second()) ? true : false;
+            }
+        )
+    ),
+    invDelta_(),
+    indexPairs_(patchIndexPairs()),
+    R_(interpolateOrRead<symmTensor>("R", dict, interpolateR_)),
+    LundWuSquires_(computeLundWuSquires()),
+    UMean_(interpolateOrRead<vector>("UMean", dict, interpolateUMean_)),
+    Lbak_
+    (
+        dict.getCheck<tensor>
+        (
+            "L",
+            [&](const tensor& l){return tiny_ < cmptMin(l) ? true : false;}
+        )
+    ),
+    L_(convertScalesToGridUnits(Lbak_)),
+    const1FSM_
+    (
+        (variant_ == variantType::FORWARD_STEPWISE)
+      ? dict.getOrDefault<scalar>
+        (
+            "const1FSM",
+           -0.25*constant::mathematical::pi
+        )
+      : Zero
+    ),
+    const2FSM_
+    (
+        (variant_ == variantType::FORWARD_STEPWISE)
+      ? dict.getOrDefault<scalar>
+        (
+            "const2FSM",
+           -0.5*constant::mathematical::pi
+        )
+      : Zero
+    ),
+    constList1FSM_
+    (
+        (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())
+      ? dict.getOrDefault<List<List<scalar>>>
+        (
+            "randomBox",
+            fillRandomBox() // First time-step
+        )
+      :
+        fillRandomBox()
+    ),
+    filterCoeffs_(computeFilterCoeffs()),
+    filteredRandomBox_
+    (
+        pTraits<vector>::nComponents,
+        List<scalar>(planeDivisions_.first()*planeDivisions_.second(), Zero)
+    ),
+    U0_
+    (
+        (variant_ == variantType::FORWARD_STEPWISE)
+      ? generateRandomSet<vectorField, vector>(patch().size())
+      : vectorField()
+    ),
+    computeVariant
+    (
+        (variant_ == variantType::FORWARD_STEPWISE)
+      ? &turbulentDigitalFilterInletFvPatchVectorField::computeReducedDFM
+      : &turbulentDigitalFilterInletFvPatchVectorField::computeDFM
+    )
+{
+    if (isCorrectedFlowRate_)
+    {
+        initialFlowRate_ = computeInitialFlowRate();
+    }
+
+    // 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;
+    }
+
+    #ifdef FULLDEBUG
+    Info<< "Ends: Resource acquisition/initialisation for the"
+        << " synthetic turbulence boundary condition." << nl;
+    #endif
+}
+
+
+Foam::turbulentDigitalFilterInletFvPatchVectorField::
+turbulentDigitalFilterInletFvPatchVectorField
+(
+    const turbulentDigitalFilterInletFvPatchVectorField& ptf
+)
+:
+    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_),
+    mapMethod_(ptf.mapMethod_),
+    curTimeIndex_(ptf.curTimeIndex_),
+    tiny_(ptf.tiny_),
+    patchNormalSpeed_(ptf.patchNormalSpeed_),
+    modelConst_(ptf.modelConst_),
+    perturb_(ptf.perturb_),
+    initialFlowRate_(ptf.initialFlowRate_),
+    rndGen_(ptf.rndGen_),
+    planeDivisions_(ptf.planeDivisions_),
+    invDelta_(ptf.invDelta_),
+    indexPairs_(ptf.indexPairs_),
+    R_(ptf.R_),
+    LundWuSquires_(ptf.LundWuSquires_),
+    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_),
+    iNextToLastPlane_(ptf.iNextToLastPlane_),
+    randomBox_(ptf.randomBox_),
+    filterCoeffs_(ptf.filterCoeffs_),
+    filteredRandomBox_(ptf.filteredRandomBox_),
+    U0_(ptf.U0_),
+    computeVariant(ptf.computeVariant)
+{}
+
+
+Foam::turbulentDigitalFilterInletFvPatchVectorField::
+turbulentDigitalFilterInletFvPatchVectorField
+(
+    const turbulentDigitalFilterInletFvPatchVectorField& ptf,
+    const DimensionedField<vector, volMesh>& iF
+)
+:
+    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_),
+    mapMethod_(ptf.mapMethod_),
+    curTimeIndex_(-1),
+    tiny_(ptf.tiny_),
+    patchNormalSpeed_(ptf.patchNormalSpeed_),
+    modelConst_(ptf.modelConst_),
+    perturb_(ptf.perturb_),
+    initialFlowRate_(ptf.initialFlowRate_),
+    rndGen_(ptf.rndGen_),
+    planeDivisions_(ptf.planeDivisions_),
+    invDelta_(ptf.invDelta_),
+    indexPairs_(ptf.indexPairs_),
+    R_(ptf.R_),
+    LundWuSquires_(ptf.LundWuSquires_),
+    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_),
+    iNextToLastPlane_(ptf.iNextToLastPlane_),
+    randomBox_(ptf.randomBox_),
+    filterCoeffs_(ptf.filterCoeffs_),
+    filteredRandomBox_(ptf.filteredRandomBox_),
+    U0_(ptf.U0_),
+    computeVariant(ptf.computeVariant)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::turbulentDigitalFilterInletFvPatchVectorField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    if (curTimeIndex_ != db().time().timeIndex())
+    {
+        vectorField& U = *this;
+
+        computeVariant(this, U);
+
+        curTimeIndex_ = db().time().timeIndex();
+    }
+
+    fixedValueFvPatchVectorField::updateCoeffs();
+}
+
+
+void Foam::turbulentDigitalFilterInletFvPatchVectorField::write
+(
+    Ostream& os
+) const
+{
+    fixedValueFvPatchField<vector>::write(os),
+    os.writeEntry("variant", variantNames[variant_]);
+    os.writeEntry("planeDivisions", planeDivisions_);
+    os.writeEntry("L", Lbak_);
+
+    if (!interpolateR_)
+    {
+        R_.writeEntry("R", os);
+    }
+
+    if (!interpolateUMean_)
+    {
+        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_);
+
+    if (!mapMethod_.empty())
+    {
+        os.writeEntryIfDifferent<word>
+        (
+            "mapMethod",
+            "nearestCell",
+            mapMethod_
+        );
+    }
+
+    os.writeEntry("threshold", tiny_);
+    os.writeEntry("patchNormalSpeed", patchNormalSpeed_);
+    os.writeEntry("modelConst", modelConst_);
+    os.writeEntry("perturb", perturb_);
+
+    if (variant_ == variantType::FORWARD_STEPWISE)
+    {
+        os.writeEntry("const1FSM", const1FSM_);
+        os.writeEntry("const2FSM", const2FSM_);
+    }
+
+    if (isContinuous_ && Pstream::master())
+    {
+        os.writeEntry("randomBox", randomBox_);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+   makePatchTypeField
+   (
+       fvPatchVectorField,
+       turbulentDigitalFilterInletFvPatchVectorField
+   );
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/turbulentDigitalFilterInletFvPatchVectorField.H b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/turbulentDigitalFilterInletFvPatchVectorField.H
new file mode 100644
index 0000000000..2ec5ef8c7e
--- /dev/null
+++ b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/turbulentDigitalFilterInletFvPatchVectorField.H
@@ -0,0 +1,569 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::turbulentDigitalFilterFvPatchVectorField
+
+Group
+    grpInletBoundaryConditions
+
+Description
+    Velocity boundary condition generating synthetic turbulence-alike
+    time-series for LES and DES turbulent flow computations.
+
+    To this end, two synthetic turbulence generators can be chosen:
+
+    - Digital-filter method-based generator (DFM)
+
+    \verbatim
+    Klein, M., Sadiki, A., and Janicka, J.
+        A digital filter based generation of inflow data for spatially 
+        developing direct numerical or large eddy simulations,
+        Journal of Computational Physics (2003) 186(2):652-665.
+        doi:10.1016/S0021-9991(03)00090-1
+    \endverbatim
+
+    - Forward-stepwise method-based generator (FSM)
+
+    \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
+    of target statistics (mostly mean flow, Reynolds stress tensor profiles and
+    length-scale sets) are fused 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.
+
+    Random number sets ---->-|
+                             |
+                         DFM or FSM ---> New stochastic time-series consisting
+                             |           turbulence statistics
+    Turbulence statistics ->-|
+
+    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.
+
+    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.
+
+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:
+    \verbatim
+    <patchName>
+    {
+        type                turbulentDigitalFilter;
+        variant             digitalFilter;          // reducedDigitalFilter;
+        planeDivisions      (<planeDivisionsHeight> <planeDivisionsWidth>);
+        L               (<Lxu> <Lxv> <Lxw> <Lyu> <Lyv> <Lyw> <Lzu> <Lzv> <Lzw>);
+        R                   (<Rxx> <Rxy> <Rxz> <Ryy> <Ryz> <Rzz>);
+        patchNormalSpeed    <characteristic flow speed>;
+        value               uniform (0 0 0);        // mandatory placeholder
+
+        // Optional entries with default input
+        isGaussian          true;           // false    // always false for FSM
+        isFixedSeed         true;           // false
+        isContinuous        false;          // true
+        isCorrectedFlowRate true;           // false
+        interpolateR        false;          // placeholder
+        interpolateUMean    false;          // placeholder
+        isInsideMesh        false;          // placeholder
+        isTaylorHypot       true;           // placeholder
+        mapMethod           nearestCell;    // planarInterpolation
+        threshold           1e-8;
+        modelConst          -1.5707;        //-0.5*PI;
+        perturb             1e-5;
+
+        // Optional entries for only FSM with default input
+        const1FSM           -0.7854         //-0.25*PI;
+        const2FSM           -1.5707;        //-0.5*PI;
+    }
+    \endverbatim 
+
+    Among the dictionary entries, two entries can be input as patch profiles:
+
+    \verbatim
+     - Reynolds stress tensor, R
+     - Mean velocity, UMean
+    \endverbatim
+
+    Profile data and corresponding coordinates are then input in the following
+    directories:
+
+    \verbatim
+     - $FOAM_CASE/constant/boundaryData/\<patchName\>/points
+     - $FOAM_CASE/constant/boundaryData/\<patchName\>/0/\{R|UMean\}
+    \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.
+
+    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.
+
+    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.
+
+Note
+    - \c mapMethod \c planarInterpolation option requires point coordinates
+    which can form a plane, thus input point coordinates varying only in a
+    single direction will trigger error.
+    - \c adjustTimeStep = true option is not fully supported at the moment.
+
+SeeAlso
+    turbulentDFSEMInletFvPatchVectorField.C
+
+SourceFiles
+    turbulentDigitalFilterInletFvPatchVectorField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef turbulentDigitalFilterInletFvPatchVectorField_H
+#define turbulentDigitalFilterInletFvPatchVectorField_H
+
+#include "fixedValueFvPatchFields.H"
+#include "Random.H"
+#include <functional>
+#include "fieldTypes.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class pointToPointPlanarInterpolation;
+
+/*---------------------------------------------------------------------------*\
+        Class turbulentDigitalFilterInletFvPatchVectorField Declaration
+\*---------------------------------------------------------------------------*/
+
+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)
+        mutable autoPtr<pointToPointPlanarInterpolation> mapperPtr_;
+
+        //- Selected option for the synthetic turbulence generator variant
+        const enum variantType variant_;
+
+        //- Flag: correlation function form is Gaussian or Exponential
+        //  for variantType::DIGITAL_FILTER, default=Gaussian
+        //  for variantType::FORWARD_STEPWISE, default=Exponential (only option)
+        const bool isGaussian_;
+
+        //- Flag: random-number generator seed is fixed (default=true) or
+        //- generated pseudo-randomly based on clock-time per simulation
+        const bool isFixedSeed_;
+
+        //- 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: mass flow rate is corrected on turbulence plane (default=true)
+        const bool isCorrectedFlowRate_;
+
+        //- Flag: interpolate R field (default=false)
+        bool interpolateR_;
+
+        //- 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_;
+
+        //- Method for interpolation between a patch and turbulence plane
+        //- (default=nearestCell)
+        //  Options:
+        //   - nearestCell: one-to-one direct map, no interpolation
+        //   - planarInterpolation: bilinear interpolation
+        const word mapMethod_;
+
+        //- Current time index
+        label curTimeIndex_;
+
+        //- Threshold to avoid unintentional 'tiny' input scalars (default=1e-8)
+        const scalar tiny_;
+
+        //- Characteristic (e.g. bulk) mean speed of flow in the patch normal
+        //- direction [m/s]
+        const scalar patchNormalSpeed_;
+
+        //- 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)
+        const scalar perturb_;
+
+        //- Initial (first time-step) mass/vol flow rate [m^3/s]
+        scalar initialFlowRate_;
+
+        //- Random number generator
+        Random rndGen_;
+
+        //- Number of nodes on turbulence plane (e2 e3) [-]
+        const Tuple2<label, label> planeDivisions_;
+
+        //- Turbulence plane mesh size (reversed) (e2 e3) [1/m]
+        Vector2D<scalar> invDelta_;
+
+        //- Nearest cell mapping: Index pairs between patch and turbulence plane
+        const List<Pair<label>> indexPairs_;
+
+        //- Reynolds stress tensor profile (xx xy xz yy yz zz) in global
+        //- coordinates [m^2/s^2]
+        symmTensorField R_;
+
+        //- Lund-Wu-Squires transformation (Cholesky decomp.) [m/s]
+        //- Mapped onto actual mesh patch rather than turbulence plane
+        //  (Klein et al., 2003, Eq. 5)
+        symmTensorField LundWuSquires_;
+
+        //- Mean inlet velocity profile 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]
+        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_;
+
+        //- One of the two model coefficients in FSM
+        //  (Xie-Castro, 2008, the argument of the second exp func in Eq. 14)
+        const scalar const2FSM_;
+
+        //- One of the two exponential functions in FSM
+        //  (Xie-Castro, 2008, the first exponential function in Eq. 14)
+        const List<scalar> constList1FSM_;
+
+        //- One of the two exponential functions in FSM
+        //  (Xie-Castro, 2008, the first exponential function in Eq. 14)
+        const List<scalar> constList2FSM_;
+
+        //- Number of nodes in random-number box [node]
+        //- Random-number sets within box are filtered with filterCoeffs_
+        //- (e1u, e1v, e1w, e2u, e2v, e2w, e3u, e3v, e3w)
+        const List<label> lenRandomBox_;
+
+        //- Convenience factors for 2-D random-number box [-]
+        const List<label> randomBoxFactors2D_;
+
+        //- Convenience factors for 3-D randomNum box [-]
+        const List<label> randomBoxFactors3D_;
+
+        //- 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_;
+
+        //- 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 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
+
+        //- Return a reference to the patch mapper object
+        const pointToPointPlanarInterpolation& patchMapper() const;
+
+        //- Helper function to interpolate values from the boundary data or
+        //- read from dictionary
+        template<class Type>
+        tmp<Field<Type>> interpolateOrRead
+        (
+            const word& fieldName,
+            const dictionary& dict,
+            bool& interpolateField
+        ) const;
+
+        //- Helper function to interpolate values from the boundary data
+        template<class Type>
+        tmp<Field<Type>> interpolateBoundaryData
+        (
+            const word& fieldName
+        ) const;
+
+        //- Generate random-number sets obeying the standard normal distribution
+        template<class Form, class Type>
+        Form generateRandomSet(const label len);
+
+        //- Compute nearest cell index-pairs between turbulence plane and patch
+        List<Pair<label>> patchIndexPairs();
+
+        //- Check R_ (mapped on actual mesh) for mathematical domain errors
+        void checkRTensorRealisable() const;
+
+        //- Compute Lund-Wu-Squires transformation
+        //  (Klein et al., 2003, Eq. 5)
+        symmTensorField computeLundWuSquires() const;
+
+        //- Compute patch-normal into the domain
+        vector computePatchNormal() 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;
+
+        //- Resource allocation functions for the convenience factors
+        List<label> initLenRandomBox() const;
+        List<label> initBoxFactors2D() const;
+        List<label> initBoxFactors3D() const;
+        List<label> initBoxPlaneFactors() const;
+
+        //- Compute various convenience factors for random-number box
+        List<List<scalar>> fillRandomBox();
+
+        //- Compute filter coeffs once per simulation
+        //  (Klein et al., 2003, Eq. 14)
+        List<List<scalar>> computeFilterCoeffs() const;
+
+        //- Discard current time-step random-box plane (closest to patch) by
+        //- shifting from the back to the front, and dd new plane to the back
+        void rndShiftRefill();
+
+        //- Map two-point correlated random-number sets on patch based on chosen
+        //- mapping method
+        void mapFilteredRandomBox(vectorField& U);
+
+        //- Map R_ on patch
+        void embedOnePointCorrs(vectorField& U) const;
+
+        //- Map UMean_ on patch
+        void embedMeanVelocity(vectorField& U) const;
+
+        //- Correct mass/vol flow rate in (only) streamwise direction
+        void correctFlowRate(vectorField& U) 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);
+
+
+public:
+
+   //- Runtime type information
+   TypeName("turbulentDigitalFilterInlet");
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        turbulentDigitalFilterInletFvPatchVectorField
+        (
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        turbulentDigitalFilterInletFvPatchVectorField
+        (
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //- turbulentDigitalFilterInletFvPatchVectorField onto a new patch
+        turbulentDigitalFilterInletFvPatchVectorField
+        (
+            const turbulentDigitalFilterInletFvPatchVectorField&,
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        turbulentDigitalFilterInletFvPatchVectorField
+        (
+            const turbulentDigitalFilterInletFvPatchVectorField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchVectorField> clone() const
+        {
+            return tmp<fvPatchVectorField>
+            (
+                new turbulentDigitalFilterInletFvPatchVectorField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        turbulentDigitalFilterInletFvPatchVectorField
+        (
+            const turbulentDigitalFilterInletFvPatchVectorField&,
+            const DimensionedField<vector, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchVectorField> clone
+        (
+            const DimensionedField<vector, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchVectorField>
+            (
+                new turbulentDigitalFilterInletFvPatchVectorField(*this, iF)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~turbulentDigitalFilterInletFvPatchVectorField() = default;
+
+
+    // Member Functions
+
+        // Evaluation functions
+
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs();
+
+
+       //- Write
+       virtual void write(Ostream&) const;
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "turbulentDigitalFilterInletFvPatchVectorFieldTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/turbulentDigitalFilterInletFvPatchVectorFieldTemplates.C b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/turbulentDigitalFilterInletFvPatchVectorFieldTemplates.C
new file mode 100644
index 0000000000..f43b2a2d9f
--- /dev/null
+++ b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/turbulentDigitalFilterInletFvPatchVectorFieldTemplates.C
@@ -0,0 +1,127 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 - 2019 OpenCFD Ltd.
+     \\/     M anipulation  | 
+-------------------------------------------------------------------------------
+                              Copyright (C) 2015 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "pointToPointPlanarInterpolation.H"
+#include "Time.H"
+#include "IFstream.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>>
+Foam::turbulentDigitalFilterInletFvPatchVectorField::interpolateOrRead
+(
+    const word& fieldName,
+    const dictionary& dict,
+    bool& interpolateField
+) const
+{
+    if (dict.found(fieldName))
+    {
+        tmp<Field<Type>> tFld
+        (
+            new Field<Type>
+            (
+                fieldName,
+                dict,
+                this->patch().size()
+            )
+        );
+
+        interpolateField = false;
+        return tFld;
+    }
+    else
+    {
+        interpolateField = true;
+        return interpolateBoundaryData<Type>(fieldName);
+    }
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>>
+Foam::turbulentDigitalFilterInletFvPatchVectorField::interpolateBoundaryData
+(
+    const word& fieldName
+) const
+{
+    const word& patchName = this->patch().name();
+
+    fileName valsFile
+    (
+        fileHandler().filePath
+        (
+            fileName
+            (
+                this->db().time().path()
+               /this->db().time().caseConstant()
+               /"boundaryData"
+               /patchName
+               /"0"
+               /fieldName
+            )
+        )
+    );
+
+    autoPtr<ISstream> isPtr
+    (
+        fileHandler().NewIFstream
+        (
+            valsFile
+        )
+    );
+
+    Field<Type> vals(isPtr());
+
+    Info<< "Turbulent DFM/FSM patch " << patchName
+        << ": Interpolating field " << fieldName
+        << " from " << valsFile << endl;
+
+    return patchMapper().interpolate(vals);
+}
+
+
+template<class Form, class Type>
+Form Foam::turbulentDigitalFilterInletFvPatchVectorField::generateRandomSet
+(
+    const label len
+)
+{
+    Form randomSet(len);
+
+    std::generate
+    (
+        randomSet.begin(),
+        randomSet.end(),
+        [&]{return rndGen_.GaussNormal<Type>();}
+    );
+
+    return randomSet;
+}
+
+// ************************************************************************* //