diff --git a/src/functionObjects/field/Make/files b/src/functionObjects/field/Make/files
index 649176d1f4..580d9a3340 100644
--- a/src/functionObjects/field/Make/files
+++ b/src/functionObjects/field/Make/files
@@ -118,4 +118,6 @@ stabilityBlendingFactor/stabilityBlendingFactor.C
 
 interfaceHeight/interfaceHeight.C
 
+STDMD/STDMD.C
+
 LIB = $(FOAM_LIBBIN)/libfieldFunctionObjects
diff --git a/src/functionObjects/field/Make/options b/src/functionObjects/field/Make/options
index c24acfb632..5d3d5f11d6 100644
--- a/src/functionObjects/field/Make/options
+++ b/src/functionObjects/field/Make/options
@@ -1,4 +1,5 @@
 EXE_INC = \
+    -I$(LIB_SRC)/OpenFOAM/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/fileFormats/lnInclude \
     -I$(LIB_SRC)/surfMesh/lnInclude \
@@ -22,6 +23,7 @@ EXE_INC = \
     -I$(LIB_SRC)/phaseSystemModels/reactingEulerFoam/phaseSystems/lnInclude
 
 LIB_LIBS = \
+    -lOpenFOAM \
     -lfiniteVolume \
     -lfileFormats \
     -lsurfMesh \
diff --git a/src/functionObjects/field/STDMD/STDMD.C b/src/functionObjects/field/STDMD/STDMD.C
new file mode 100644
index 0000000000..d550a4cefb
--- /dev/null
+++ b/src/functionObjects/field/STDMD/STDMD.C
@@ -0,0 +1,1280 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "STDMD.H"
+#include "addToRunTimeSelectionTable.H"
+
+using namespace Foam::constant::mathematical;
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+    defineTypeNameAndDebug(STDMD, 0);
+    addToRunTimeSelectionTable(functionObject, STDMD, dictionary);
+}
+}
+
+
+const Foam::Enum
+<
+    Foam::functionObjects::STDMD::modeSorterType
+>
+Foam::functionObjects::STDMD::modeSorterTypeNames
+({
+    { modeSorterType::KIEWAT , "kiewat" },
+    { modeSorterType::KOU_ZHANG , "kouZhang" },
+    { modeSorterType::FIRST_SNAPSHOT, "firstSnap" }
+});
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::scalar Foam::functionObjects::STDMD::parnorm
+(
+    const RMatrix& colVector
+) const
+{
+    #ifdef FULLDEBUG
+    // Check if the given RectangularMatrix is effectively a column vector
+    if (colVector.n() != 1)
+    {
+        FatalErrorInFunction
+            << "  # Input matrix is not a column vector. #"
+            << abort(FatalError);
+    }
+    #endif
+    const bool noSqrt = true;
+    scalar result(colVector.columnNorm(0, noSqrt));
+    reduce(result, sumOp<scalar>());
+
+    // Heuristic addition to avoid very small or zero norm
+    return max(SMALL, Foam::sqrt(result));
+}
+
+
+void Foam::functionObjects::STDMD::snapshot()
+{
+    bool processed = false;
+    processed = processed || getSnapshotType<scalar>();
+    processed = processed || getSnapshotType<vector>();
+    processed = processed || getSnapshotType<sphericalTensor>();
+    processed = processed || getSnapshotType<symmTensor>();
+    processed = processed || getSnapshotType<tensor>();
+
+    if (!processed)
+    {
+        FatalErrorInFunction
+            << "  # Unknown type of input field during snapshot loading = "
+            << fieldName_ << " #" << nl
+            << "  # Do you execute required functionObjects "
+            << "before executing STDMD, e.g. mapFields?"
+            << abort(FatalError);
+    }
+}
+
+
+void Foam::functionObjects::STDMD::init()
+{
+    bool processed = false;
+    processed = processed || getComps<scalar>();
+    processed = processed || getComps<vector>();
+    processed = processed || getComps<sphericalTensor>();
+    processed = processed || getComps<symmTensor>();
+    processed = processed || getComps<tensor>();
+
+    if (!processed)
+    {
+        FatalErrorInFunction
+            << "  # Unknown type of input field during initialisation = "
+            << fieldName_ << " #" << nl
+            << "  # Do you execute required functionObjects "
+            << "before executing STDMD, e.g. mapFields?"
+            << abort(FatalError);
+    }
+
+    nSnap_ = nComps_*mesh_.nCells();
+
+    if (nSnap_ <= 0)
+    {
+        FatalErrorInFunction
+            << "  # Zero-size input field = " << fieldName_ << " #"
+            << abort(FatalError);
+    }
+
+    currIndex_ = 0;
+    zNorm_ = 0;
+    ezNorm_ = 0;
+
+    z_ = RMatrix(2*nSnap_, 1, Zero);
+    ez_ = z_;
+    X1_ = RMatrix(nSnap_, 1, Zero);
+    Qz_ = z_;
+    Gz_ = SMatrix(1);
+
+    RxInv_.clear();
+    Ap_.clear();
+    oEVecs_.clear();
+    oEVals_.clear();
+    oAmps_.clear();
+    oFreqs_.clear();
+    iFreqs_.clear();
+    oMags_.clear();
+    iMags_.clear();
+
+    initialised_ = true;
+}
+
+
+void Foam::functionObjects::STDMD::initBasis()
+{
+    std::copy(z_.cbegin(), z_.cbegin() + nSnap_, X1_.begin());
+
+    zNorm_ = parnorm(z_);
+    Qz_ = z_/zNorm_;
+    Gz_(0,0) = sqr(zNorm_);
+}
+
+
+void Foam::functionObjects::STDMD::GramSchmidt()
+{
+    ez_ = z_;
+    RMatrix dz(Qz_.n(), 1, Zero);
+
+    for (label i = 0; i < nGramSchmidt_; ++i)
+    {
+        dz = Qz_ & ez_;
+        reduce(dz, sumOp<RMatrix>());
+        ez_ -= Qz_*dz;
+    }
+}
+
+
+void Foam::functionObjects::STDMD::expandBasis()
+{
+    Log<< tab << "# " << name() << ":"
+        << " Expanding orthonormal basis for field = " << fieldName_ << " #"
+        << endl;
+
+    // Stack a column (ez_/ezNorm_) to Qz_
+    Qz_.resize(Qz_.m(), Qz_.n() + 1);
+    Qz_.subColumn(Qz_.n() - 1) = ez_/ezNorm_;
+
+    // Stack a row (Zero) and column (Zero) to Gz_
+    Gz_.resize(Gz_.m() + 1);
+}
+
+
+void Foam::functionObjects::STDMD::updateGz()
+{
+    RMatrix zTilde(Qz_ & z_);
+    reduce(zTilde, sumOp<RMatrix>());
+
+    const SMatrix zTildes(zTilde^zTilde);
+
+    Gz_ += zTildes;
+}
+
+
+void Foam::functionObjects::STDMD::compressBasis()
+{
+    Log<< tab << "# " << name() << ":"
+        << " Compressing orthonormal basis for field = " << fieldName_ << " #"
+        << endl;
+
+    RMatrix qz;
+
+    if (Pstream::master())
+    {
+        const bool symmetric = true;
+        const EigenMatrix<scalar> EM(Gz_, symmetric);
+        const SquareMatrix<scalar>& EVecs = EM.EVecs();
+        DiagonalMatrix<scalar> EVals(EM.EValsRe());
+
+        if (testEigen_)
+        {
+            testEigenvalues(Gz_, EVals);
+            testEigenvectors(Gz_, EVals, EVecs);
+        }
+
+        // Sort eigenvalues in descending order, and track indices
+        const auto descend = [&](scalar a, scalar b){ return a > b; };
+        const List<label> permut(EVals.sortPermutation(descend));
+        EVals.applyPermutation(permut);
+        EVals.resize(EVals.size() - 1);
+
+        // Update Gz_
+        Gz_ = SMatrix(maxRank_, Zero);
+        Gz_.diag(EVals);
+
+        qz.resize(Qz_.n(), maxRank_);
+        for (label i = 0; i < maxRank_; ++i)
+        {
+            qz.subColumn(i) = EVecs.subColumn(permut[i]);
+        }
+    }
+    Pstream::scatter(Gz_);
+    Pstream::scatter(qz);
+
+    // Update Qz_
+    Qz_ = Qz_*qz;
+}
+
+
+void Foam::functionObjects::STDMD::calcAp()
+{
+    Log<< tab << "# " << name() << ": Computing Ap matrix #" << endl;
+
+    Log<< tab << "# " << name() << ": Computing local Rx #" << endl;
+    RMatrix Rx;
+    {
+        QRMatrix<RMatrix> QRM
+        (
+            QzUH_,
+            QRMatrix<RMatrix>::outputTypes::REDUCED_R,
+            QRMatrix<RMatrix>::colPivoting::FALSE
+        );
+        Rx = QRM.R();
+    }
+    Rx.round();
+
+    RMatrix A1; // Convenience objects A1, A2, A3 to compute Ap
+    if (Pstream::parRun())
+    {
+        // Parallel direct tall-skinny QR decomposition
+        // (BGD:Fig. 5) & (DGHL:Fig. 2)
+        // Tests revealed that the distribution of Qz_ does not affect
+        // the final outcome of TSQR decomposition up to sign
+
+        Log<< tab << "# " << name() << ": Gathering all local Rx #" << endl;
+        List<RMatrix> RxList(Pstream::nProcs());
+        RxList[Pstream::myProcNo()] = Rx;
+        Pstream::gatherList(RxList);
+        Pstream::scatterList(RxList);
+
+        // Create a global object Rx
+        if (Pstream::master())
+        {
+            // Determine the row size of the global Rx
+            label mRowSum = 0;
+            forAll(RxList, i)
+            {
+                mRowSum += RxList[i].m();
+            }
+
+            Rx.resize(mRowSum, Rx.n());
+
+            Log<< tab << "# " << name()
+                << ": Populating the global Rx #" << endl;
+            label m = 0;
+            for (label i = 0; i < Pstream::nProcs(); ++i)
+            {
+                label mRows = RxList[i].m();
+
+                Rx.subMatrix(m, 0, mRows) = RxList[i];
+
+                m += mRows;
+            }
+
+            Log<< tab << "# " << name()
+                << ": Computing the parallel-direct tall-skinny QR decomp. #"
+                << endl;
+            QRMatrix<RMatrix> QRM
+            (
+                Rx,
+                QRMatrix<RMatrix>::outputTypes::REDUCED_R,
+                QRMatrix<RMatrix>::storeMethods::IN_PLACE,
+                QRMatrix<RMatrix>::colPivoting::FALSE
+            );
+            Rx.round();
+
+            Log<< tab << "# " << name()
+                << ": Computing Moore-Penrose pseudo-inverse of Rx #"
+                << endl;
+            RxInv_ = pinv(Rx);
+
+            Log<< tab << "# " << name() << ": Computing Gx #" << endl;
+            const RMatrix Gx(Rx*(Gz_^Rx));
+
+            Log<< tab << "# " << name()
+                << ": Computing Moore-Penrose pseudo-inverse of Gx #"
+                << endl;
+            const RMatrix GxInv(pinv(Gx));
+
+            Log<< tab << "# " << name() << ": Computing A1 #" << endl;
+            A1 = RxInv_*GxInv;
+        }
+        Pstream::scatter(RxInv_);
+        Pstream::scatter(A1);
+
+        Log<< tab << "# " << name() << ": Computing A2 #" << endl;
+        SMatrix A2(QzUH_ & QzUH_);
+        reduce(A2, sumOp<SMatrix>());
+
+        Log<< tab << "# " << name() << ": Computing A3 #" << endl;
+        SMatrix A3(QzUH_ & QzLH_);
+        reduce(A3, sumOp<SMatrix>());
+
+        Log<< tab << "# " << name() << ": Computing Ap #" << endl;
+        // by optimized matrix chain multiplication
+        // obtained by dynamic programming memoisation
+        Ap_ = SMatrix(RxInv_ & ((A3*(Gz_*A2))*A1));
+    }
+    else
+    {
+        {
+            Log<< tab << "# " << name()
+                << ": Computing Moore-Penrose pseudo-inverse of Rx #"
+                << endl;
+            RxInv_ = pinv(Rx);
+
+            Log<< tab << "# " << name() << ": Computing Gx #" << endl;
+            const RMatrix Gx(Rx*(Gz_^Rx));
+
+            Log<< tab << "# " << name()
+                << ": Computing Moore-Penrose pseudo-inverse of Gx #"
+                << endl;
+            const RMatrix GxInv(pinv(Gx));
+
+            Log<< tab << "# " << name() << ": Computing A1 #" << endl;
+            A1 = RxInv_*GxInv;
+        }
+
+        Log<< tab << "# " << name() << ": Computing A2 #" << endl;
+        SMatrix A2(QzUH_ & QzUH_);
+
+        Log<< tab << "# " << name() << ": Computing A3 #" << endl;
+        SMatrix A3(QzUH_ & QzLH_);
+
+        Log<< tab << "# " << name() << ": Computing Ap #" << endl;
+        // by optimized matrix chain multiplication
+        // obtained by dynamic programming memoisation
+        Ap_ = SMatrix(RxInv_ & ((A3*(Gz_*A2))*A1));
+    }
+}
+
+
+void Foam::functionObjects::STDMD::calcEigen()
+{
+    Log<< tab << "# " << name() << ": Computing eigendecomposition #" << endl;
+
+    // Compute eigenvalues, and clip eigenvalues with values < minMagEVal_
+    if (Pstream::master())
+    {
+        // Replace Ap by eigenvalues (in-place eigendecomposition)
+        const EigenMatrix<scalar> EM(Ap_);
+        const DiagonalMatrix<scalar>& evalsRe = EM.EValsRe();
+        const DiagonalMatrix<scalar>& evalsIm = EM.EValsIm();
+
+        oEVals_.resize(evalsRe.size());
+        oEVecs_ = RectangularMatrix<complex>(EM.complexEVecs());
+
+        // Convert scalar eigenvalue pairs to complex eigenvalues
+        label i = 0;
+        for (auto& eval : oEVals_)
+        {
+            eval = complex(evalsRe[i], evalsIm[i]);
+            ++i;
+        }
+
+        if (testEigen_)
+        {
+            testEigenvalues(Ap_, evalsRe);
+            testEigenvectors(Ap_, oEVals_, oEVecs_);
+        }
+    }
+}
+
+
+bool Foam::functionObjects::STDMD::close
+(
+    const scalar s1,
+    const scalar s2,
+    const scalar absTol,
+    const scalar relTol
+) const
+{
+    if (s1 == s2) return true;
+
+    return (mag(s1 - s2) <= max(relTol*max(mag(s1), mag(s2)), absTol));
+}
+
+
+void Foam::functionObjects::STDMD::testEigenvalues
+(
+    const SquareMatrix<scalar>& A,
+    const DiagonalMatrix<scalar>& EValsRe
+) const
+{
+    const scalar trace = A.trace();
+    // Imaginary part of complex conjugates cancel each other
+    const scalar EValsSum = sum(EValsRe);
+
+    const bool clse = close(EValsSum - trace, 0, absTol_, relTol_);
+
+    if (clse)
+    {
+        Info<< tab
+            << "  ## PASS: Eigenvalues ##" << endl;
+    }
+    else
+    {
+        Info<< tab
+            << "  ## INCONCLUSIVE: Eigenvalues ##" << nl
+            << "  # sum(eigenvalues) = " << EValsSum << nl
+            << "  # trace(A) = " << trace << nl
+            << "  # sum(eigenvalues) - trace(A) ~ 0 ?= "
+            << EValsSum - trace << nl
+            << "  #######################"
+            << endl;
+
+        if (dumpEigen_)
+        {
+            Info<< tab
+                << "  ## Operands ##" << nl
+                << "  # eigenvalues:" << nl << EValsRe << nl
+                << "  # input matrix A:" << nl << A << nl
+                << "  ##############"
+                << endl;
+        }
+    }
+}
+
+
+void Foam::functionObjects::STDMD::testEigenvectors
+(
+    const SquareMatrix<scalar>& A,
+    const DiagonalMatrix<scalar>& EValsRe,
+    const SquareMatrix<scalar>& EVecs
+) const
+{
+    unsigned nInconclusive = 0;
+
+    for (label i = 0; i < A.n(); ++i)
+    {
+        const RectangularMatrix<scalar> EVec(EVecs.subColumn(i));
+        const scalar EVal = EValsRe[i];
+        const RectangularMatrix<scalar> leftSide(A*EVec);
+        const RectangularMatrix<scalar> rightSide(EVal*EVec);
+        const scalar rslt = (leftSide - rightSide).norm();
+
+        const bool clse = close(rslt, 0, absTol_, relTol_);
+
+        if (!clse)
+        {
+            Info<< tab
+                << "  ## INCONCLUSIVE: Eigenvector ##" << nl
+                << "  # (A & EVec - EVal*EVec).norm() ~ 0 ?= " << rslt << nl
+                << "  ##################################"
+                << endl;
+
+            if (dumpEigen_)
+            {
+                Info<< tab
+                    << "  ## Operands ##" << nl
+                    << "  # eigenvalue:" << nl << EVal << nl
+                    << "  # input matrix A:" << nl << A << nl
+                    << "  # eigenvector:" << nl << EVec << nl
+                    << "  # (A & EVec):" << nl << leftSide << nl
+                    << "  # (EVal*EVec):" << nl << rightSide << nl
+                    << "  ##############"
+                    << endl;
+            }
+
+            ++nInconclusive;
+        }
+    }
+
+    if (!nInconclusive)
+    {
+        Info<< tab
+            << "  ## PASS: Eigenvectors ##" << endl;
+    }
+}
+
+
+void Foam::functionObjects::STDMD::testEigenvectors
+(
+    const SquareMatrix<scalar>& A,
+    const List<complex>& EVals,
+    const RectangularMatrix<complex>& EVecs
+) const
+{
+    SquareMatrix<complex> B(A.m());
+    auto convertToComplex = [&](const scalar& val) { return complex(val); };
+    std::transform
+    (
+        A.cbegin(),
+        A.cend(),
+        B.begin(),
+        convertToComplex
+    );
+
+    unsigned nInconclusive = 0;
+
+    for (label i = 0; i < B.n(); ++i)
+    {
+        const RectangularMatrix<complex> EVec(EVecs.subColumn(i));
+        const complex EVal = EVals[i];
+        const RectangularMatrix<complex> leftSide(B*EVec);
+        const RectangularMatrix<complex> rightSide(EVal*EVec);
+        const scalar rslt = (leftSide - rightSide).norm();
+
+        const bool clse = close(rslt, 0, absTol_, relTol_);
+
+        if (!clse)
+        {
+            Info<< tab
+                << "  ## INCONCLUSIVE: Eigenvector ##" << nl
+                << "  # (A & EVec - EVal*EVec).norm() ~ 0 ?= " << rslt << nl
+                << "  ##################################"
+                << endl;
+
+            if (dumpEigen_)
+            {
+                Info<< tab
+                    << "  ## Operands ##" << nl
+                    << "  # eigenvalue:" << nl << EVal << nl
+                    << "  # input matrix A:" << nl << A << nl
+                    << "  # eigenvector:" << nl << EVec << nl
+                    << "  # (A & EVec):" << nl << leftSide << nl
+                    << "  # (EVal*EVec):" << nl << rightSide << nl
+                    << "  ##############"
+                    << endl;
+            }
+
+            ++nInconclusive;
+        }
+    }
+
+    if (!nInconclusive)
+    {
+        Info<< tab
+            << "  ## PASS: Eigenvectors ##" << endl;
+    }
+}
+
+
+void Foam::functionObjects::STDMD::filterEVals()
+{
+    Log<< tab << "# " << name() << ": Filtering eigenvalues #" << endl;
+
+    if (Pstream::master())
+    {
+        List<complex> cpEVals(oEVals_.size());
+        auto it =
+            std::copy_if
+            (
+                oEVals_.cbegin(),
+                oEVals_.cend(),
+                cpEVals.begin(),
+                [&](const complex& x){ return minMagEVal_ < mag(x); }
+            );
+        cpEVals.resize(std::distance(cpEVals.begin(), it));
+
+        if (cpEVals.size() == 0)
+        {
+            WarningInFunction
+                << "No eigenvalue with mag(eigenvalue) larger than "
+                << "minMagEVal_ = " << minMagEVal_ << " was found."
+                << endl;
+        }
+        else
+        {
+            oEVals_ = cpEVals;
+        }
+    }
+    Pstream::scatter(oEVals_);
+    Pstream::scatter(oEVecs_);
+}
+
+
+void Foam::functionObjects::STDMD::calcFreqs()
+{
+    Log<< tab << "# " << name() << ": Computing frequencies #" << endl;
+
+    if (Pstream::master())
+    {
+        oFreqs_.resize(oEVals_.size());
+
+        // Frequency equation (K:Eq. 81)
+        auto fEq =
+            [&](const complex& EVal)
+            {
+                return Foam::log(max(EVal, complex(SMALL))).imag()/(twoPi*dt_);
+            };
+
+        // Compute frequencies
+        std::transform(oEVals_.cbegin(), oEVals_.cend(), oFreqs_.begin(), fEq);
+    }
+}
+
+
+void Foam::functionObjects::STDMD::calcFreqI()
+{
+    Log<< tab << "# " << name() << ": Computing frequency indices #" << endl;
+
+    if (Pstream::master())
+    {
+        // Initialise iterator by the first search
+        auto margin = [&](const scalar& x){ return (fMin_ <= x && x < fMax_); };
+
+        auto it = std::find_if(oFreqs_.cbegin(), oFreqs_.cend(), margin);
+
+        while (it != oFreqs_.end())
+        {
+            iFreqs_.append(std::distance(oFreqs_.cbegin(), it));
+            it = std::find_if(std::next(it), oFreqs_.cend(), margin);
+        }
+    }
+    Pstream::scatter(oFreqs_);
+    Pstream::scatter(iFreqs_);
+}
+
+
+void Foam::functionObjects::STDMD::calcAmps()
+{
+    Log<< tab << "# " << name() << ": Computing amplitudes #" << endl;
+
+    RMatrix temp((RxInv_.T()^QzUH_)*X1_);
+    reduce(temp, sumOp<RMatrix>());
+
+    if (Pstream::master())
+    {
+        oAmps_.resize(temp.m());
+        const RCMatrix pinvEVecs(pinv(oEVecs_));
+
+        // oAmps_ = pinvEVecs*temp;
+        for (label i = 0; i < oAmps_.size(); ++i)
+        {
+            for (label k = 0; k < temp.m(); ++k)
+            {
+                oAmps_[i] += pinvEVecs(i, k)*temp(k, 0);
+            }
+        }
+    }
+    Pstream::scatter(oAmps_);
+}
+
+
+void Foam::functionObjects::STDMD::calcMags()
+{
+    Log<< tab << "# " << name() << ": Computing magnitudes #" << endl;
+
+    if (Pstream::master())
+    {
+        oMags_.resize(oAmps_.size());
+
+        Log<< tab << "# " << name() << ": Sorting modes with ";
+
+        switch (modeSorter_)
+        {
+            case modeSorterType::FIRST_SNAPSHOT:
+            {
+                Log<< "method of first snapshot #" << endl;
+
+                std::transform
+                (
+                    oAmps_.cbegin(),
+                    oAmps_.cend(),
+                    oMags_.begin(),
+                    [&](const complex& val){ return mag(val); }
+                );
+                break;
+            }
+
+            case modeSorterType::KIEWAT:
+            {
+                Log<< "modified weighted amplitude scaling method #" << endl;
+
+                //  Eigendecomposition returns eigenvectors with
+                //  the unity norm, hence modeNorm = 1
+                const scalar modeNorm = 1;
+                const scalar pr = 1;
+                List<scalar> w(currIndex_);
+                std::iota(w.begin(), w.end(), 1);
+                w = sin(twoPi/currIndex_*(w - 1 - 0.25*currIndex_))*pr + pr;
+
+                label i = 0;
+                for (const auto& amp : oAmps_)
+                {
+                    oMags_[i] = sorter(w, amp, oEVals_[i], modeNorm);
+                    ++i;
+                }
+                break;
+            }
+
+            case modeSorterType::KOU_ZHANG:
+            {
+                Log<< "weighted amplitude scaling method #" << endl;
+
+                const scalar modeNorm = 1;
+                const List<scalar> w(currIndex_, 1.0);
+                label i = 0;
+                for (const auto& amp : oAmps_)
+                {
+                    oMags_[i] = sorter(w, amp, oEVals_[i], modeNorm);
+                    ++i;
+                }
+                break;
+            }
+
+            default:
+                break;
+        }
+    }
+}
+
+
+void Foam::functionObjects::STDMD::calcMagI()
+{
+    Log<< tab << "# " << name() << ": Computing magnitude indices #" << endl;
+
+    if (Pstream::master())
+    {
+        iMags_ = iFreqs_;
+
+        auto descend =
+            [&](const label i1, const label i2)
+            {
+                return !(oMags_[i1] < oMags_[i2]);
+            };
+
+        std::sort(iMags_.begin(), iMags_.end(), descend);
+    }
+    Pstream::scatter(oMags_);
+    Pstream::scatter(iMags_);
+}
+
+
+void Foam::functionObjects::STDMD::calcModes()
+{
+    Log<< tab << "# " << name() << ": Computing modes #" << endl;
+
+    // Resize the number of output variables to nModes_, if need be
+    if (nModes_ < iMags_.size())
+    {
+        iMags_.resize(nModes_);
+    }
+
+    // Compute and write out modes one by one
+    label oModeI = 0;
+    const RMatrix temp(QzUH_*RxInv_);
+
+    for (const auto& i : iMags_)
+    {
+        List<complex> mode(temp.m(), Zero);
+
+        // mode = temp*oEVecs_.subColumn(i);
+        for (label p = 0; p < mode.size(); ++p)
+        {
+            for (label q = 0; q < oEVecs_.m(); ++q)
+            {
+                mode[p] += temp(p, q)*oEVecs_(q, i);
+            }
+        }
+
+        volVectorField modeReal
+        (
+            IOobject
+            (
+                "modeReal" + std::to_string(oModeI) + fieldName_ + "_" + name(),
+                time_.timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedVector(dimless, Zero)
+        );
+
+        volVectorField modeImag
+        (
+            IOobject
+            (
+                "modeImag" + std::to_string(oModeI) + fieldName_ + "_" + name(),
+                time_.timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedVector(dimless, Zero)
+        );
+
+        label s = 0;
+        for (label k = 0; k < nComps_; ++k)
+        {
+            for (label j = 0; j < modeReal.size(); ++j)
+            {
+                const complex& m = mode[s];
+                modeReal[j].replace(k, m.real());
+                modeImag[j].replace(k, m.imag());
+                ++s;
+            }
+        }
+
+        modeReal.write();
+        modeImag.write();
+
+        ++oModeI;
+    }
+}
+
+
+Foam::scalar Foam::functionObjects::STDMD::sorter
+(
+    const List<scalar>& weight,
+    const complex& amplitude,
+    const complex& eval,
+    const scalar modeNorm
+) const
+{
+    // Omit eigenvalues with very large or very small magnitudes
+    if (!(mag(eval) < GREAT && mag(eval) > VSMALL))
+    {
+        Info<< "  Returning zero magnitude for mag(eval) = " << mag(eval)
+            << endl;
+
+        return 0.0;
+    }
+
+    // Omit eigenvalue-STDMD step combinations that pose a risk of overflow
+    if (mag(eval)*currIndex_ > sortLimiter_)
+    {
+        Info<< "  Returning zero magnitude for"
+            << " mag(eval) = " << mag(eval)
+            << " currIndex = " << currIndex_
+            << " sortLimiter = " << sortLimiter_
+            << endl;
+
+        return 0.0;
+    }
+
+    scalar magnitude = 0;
+
+    for (label j = 0; j < currIndex_; ++j)
+    {
+        magnitude += weight[j]*modeNorm*mag(amplitude*pow(eval, j + 1));
+    }
+
+    return magnitude;
+}
+
+
+void Foam::functionObjects::STDMD::writeFileHeader(Ostream& os) const
+{
+    writeHeader(os, "STDMD output");
+    writeCommented(os, "Frequency");
+    writeTabbed(os, "Magnitude");
+    writeTabbed(os, "Amplitude (real)");
+    writeTabbed(os, "Amplitude (imag)");
+    writeTabbed(os, "Eigenvalue (real)");
+    writeTabbed(os, "Eigenvalue (imag)");
+    os  << endl;
+}
+
+
+void Foam::functionObjects::STDMD::filterOutput()
+{
+    Log<< tab << "# " << name() << ": Filtering text output #" << endl;
+
+    if (Pstream::master())
+    {
+        // Filter objects according to iMags
+        filterIndexed(oEVals_, iMags_);
+        filterIndexed(oEVecs_, iMags_);
+        filterIndexed(oFreqs_, iMags_);
+        filterIndexed(oAmps_, iMags_);
+        filterIndexed(oMags_, iMags_);
+
+        // Clip objects if need be (assuming objects have the same len)
+        if (nModes_ < oFreqs_.size())
+        {
+            oEVals_.resize(nModes_);
+            oEVecs_.resize(oEVecs_.m(), nModes_);
+            oFreqs_.resize(nModes_);
+            oAmps_.resize(nModes_);
+            oMags_.resize(nModes_);
+        }
+    }
+}
+
+
+void Foam::functionObjects::STDMD::writeOutput(OFstream& os) const
+{
+    Log<< tab << "# " << name() << ": Writing text output #" << endl;
+
+    writeFileHeader(os);
+
+    for (const auto& i : labelRange(0, oFreqs_.size()))
+    {
+        os  << oFreqs_[i] << tab
+            << oMags_[i] << tab
+            << oAmps_[i].real() << tab
+            << oAmps_[i].imag() << tab
+            << oEVals_[i].real() << tab
+            << oEVals_[i].imag() << endl;
+    }
+}
+
+
+void Foam::functionObjects::STDMD::calcOutput()
+{
+    Log<< tab << "# " << name() << ":"
+        << " Starts output processing for field = " << fieldName_ << " #"
+        << endl;
+
+    // Move upper and lower halves of Qz_ to new containers
+    QzUH_ = Qz_.subMatrix(0, 0, nSnap_);
+    QzLH_ = Qz_.subMatrix(nSnap_, 0, nSnap_);
+    Qz_.clear();
+
+    calcAp();
+
+    QzLH_.clear();
+    Gz_.clear();
+
+    calcEigen();
+
+    filterEVals();
+
+    Ap_.clear();
+
+    calcFreqs();
+
+    calcFreqI();
+
+    calcAmps();
+
+    calcMags();
+
+    calcMagI();
+
+    calcModes();
+
+    // Write out unfiltered data
+    if (Pstream::master() && writeToFile_)
+    {
+        autoPtr<OFstream> osPtr =
+            createFile("uSTDMD" + fieldName_, mesh_.time().timeOutputValue());
+        OFstream& os = osPtr.ref();
+
+        writeOutput(os);
+    }
+
+    filterOutput();
+
+    // Write out filtered data
+    if (Pstream::master() && writeToFile_)
+    {
+        autoPtr<OFstream> osPtr =
+            createFile("STDMD" + fieldName_, mesh_.time().timeOutputValue());
+        OFstream& os = osPtr.ref();
+
+        writeOutput(os);
+    }
+
+    Log<< tab << "# " << name() << ":"
+        << " Ends output processing for field = " << fieldName_ << " #"
+        << endl;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::functionObjects::STDMD::STDMD
+(
+    const word& name,
+    const Time& runTime,
+    const dictionary& dict
+)
+:
+    fvMeshFunctionObject(name, runTime, dict),
+    writeFile(mesh_, name, typeName, dict, false),
+    modeSorter_
+    (
+        modeSorterTypeNames.getOrDefault
+        (
+            "modeSorter",
+            dict,
+            modeSorterType::FIRST_SNAPSHOT
+        )
+    ),
+    fieldName_(dict.get<word>("field")),
+    initialised_(false),
+    testEigen_(false),
+    dumpEigen_(false),
+    nModes_
+    (
+        dict.getCheckOrDefault<label>
+        (
+            "nModes",
+            pTraits<label>::max,
+            labelMinMax::ge(1)
+        )
+    ),
+    maxRank_
+    (
+        dict.getCheckOrDefault<label>
+        (
+            "maxRank",
+            pTraits<label>::max,
+            labelMinMax::ge(1)
+        )
+    ),
+    nGramSchmidt_
+    (
+        dict.getCheckOrDefault<label>
+        (
+            "nGramSchmidt",
+            5,
+            labelMinMax::ge(1)
+        )
+    ),
+    fMin_
+    (
+        dict.getCheckOrDefault<label>
+        (
+            "fMin",
+            0,
+            labelMinMax::ge(0)
+        )
+    ),
+    fMax_
+    (
+        dict.getCheckOrDefault<label>
+        (
+            "fMax",
+            pTraits<label>::max,
+            labelMinMax::ge(fMin_ + 1)
+        )
+    ),
+    nComps_(Zero),
+    nSnap_(Zero),
+    currIndex_(Zero),
+    sortLimiter_(500),
+    absTol_(1e-2),
+    relTol_(1e-3),
+    minBasis_(Zero),
+    dt_(Zero),
+    minMagEVal_(Zero),
+    zNorm_(Zero),
+    ezNorm_(Zero),
+    z_(),
+    ez_(),
+    X1_(),
+    Qz_(),
+    Gz_(),
+    QzUH_(),
+    QzLH_(),
+    RxInv_(),
+    Ap_(),
+    oEVecs_(),
+    oEVals_(),
+    oAmps_(),
+    oFreqs_(),
+    iFreqs_(),
+    oMags_(),
+    iMags_()
+{
+    Info<< nl
+        << "    ### This STDMD release is the beta release ###" << nl
+        << "    Therefore small-to-medium changes in input/output interfaces "
+        << "and internal structures should be expected in the next versions."
+        << endl;
+
+    // Check if varying or fixed time-step computation
+    if (runTime.isAdjustTimeStep())
+    {
+        WarningInFunction
+            << "  # STDMD: Viable only for fixed time-step computations. #"
+            << endl;
+    }
+
+    // Check if mesh topology changing
+    if (mesh_.topoChanging())
+    {
+        FatalErrorInFunction
+            << "  # STDMD: Available only for non-changing mesh topology. #"
+            << abort(FatalError);
+    }
+
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::functionObjects::STDMD::read(const dictionary& dict)
+{
+    fvMeshFunctionObject::read(dict);
+
+    testEigen_ = dict.getOrDefault<bool>("testEigen", false);
+
+    if (testEigen_)
+    {
+        dumpEigen_ = dict.getOrDefault<bool>("dumpEigen", false);
+    }
+
+    sortLimiter_ =
+        dict.getCheckOrDefault<scalar>
+        (
+            "sortLimiter",
+            500,
+            scalarMinMax::ge(1)
+        );
+
+    absTol_ =
+        dict.getCheckOrDefault<scalar>
+        (
+            "absTol",
+            1e-2,
+            scalarMinMax::ge(0)
+        );
+
+    relTol_ =
+        dict.getCheckOrDefault<scalar>
+        (
+            "relTol",
+            1e-3,
+            scalarMinMax::ge(0)
+        );
+
+    minBasis_ =
+        dict.getCheckOrDefault<scalar>
+        (
+            "minBasis",
+            1e-8,
+            scalarMinMax::ge(0)
+        );
+
+    dt_ =
+        dict.getCheckOrDefault<scalar>
+        (
+            "stdmdInterval",
+            (
+                dict.getCheck<label>("executeInterval", labelMinMax::ge(1))
+               *mesh_.time().deltaT().value()
+            ),
+            scalarMinMax::ge(0)
+        );
+
+    minMagEVal_ =
+        dict.getCheckOrDefault<scalar>
+        (
+            "minEVal",
+            1e-8,
+            scalarMinMax::ge(0)
+        );
+
+    writeToFile_ = dict.getOrDefault<bool>("writeToFile", true);
+
+    Info<< nl << name() << ":" << nl
+        << tab << "# Input is read for field = "
+        << fieldName_  << " #" << nl << endl;
+
+    return true;
+}
+
+
+bool Foam::functionObjects::STDMD::execute()
+{
+    Log << type() << " " << name() << " execute:" << endl;
+
+    // STDMD incremental orthonormal basis update (K:Fig. 15)
+    snapshot();
+
+    if (currIndex_ == 1)
+    {
+        initBasis();
+    }
+
+    if (1 < currIndex_)
+    {
+        GramSchmidt();
+
+        // Check basis for z_ and, if necessary, expand Qz_ and Gz_
+        zNorm_ = parnorm(z_);
+        ezNorm_ = parnorm(ez_);
+
+        if (minBasis_ < ezNorm_/zNorm_)
+        {
+            expandBasis();
+        }
+
+        updateGz();
+
+        // Compress the orthonormal basis if required
+        if (maxRank_ < Qz_.n())
+        {
+            compressBasis();
+        }
+    }
+    ++currIndex_;
+
+    Log<< tab << "# " << name() << ":"
+        << " Execution index = " << currIndex_
+        << " for field = " << fieldName_ << " #"
+        << endl;
+
+    return true;
+}
+
+
+bool Foam::functionObjects::STDMD::write()
+{
+    Log << type() << " " << name() << " write:" << endl;
+
+    if (currIndex_ < 2)
+    {
+        WarningInFunction
+            << "  # STDMD needs at least three snapshots to produce output #"
+            << nl
+            << "  # Only " << currIndex_ + 1 << " snapshots are available #"
+            << nl
+            << "  # Skipping STDMD output calculation and write #"
+            << endl;
+
+        return false;
+    }
+
+    // STDMD mode evaluation (K:Fig. 16)
+    calcOutput();
+
+    // Restart STDMD incremental orthonormal basis update
+    initialised_ = false;
+
+    mesh_.time().printExecutionTime(Info);
+
+    return true;
+}
+
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/STDMD/STDMD.H b/src/functionObjects/field/STDMD/STDMD.H
new file mode 100644
index 0000000000..8133004baf
--- /dev/null
+++ b/src/functionObjects/field/STDMD/STDMD.H
@@ -0,0 +1,644 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::functionObjects::STDMD
+
+Group
+    grpFieldFunctionObjects
+
+Description
+    (Beta Release) STDMD (i.e. Streaming Total Dynamic Mode Decomposition) is
+    a variant of a data-driven dimensionality reduction method.
+
+    STDMD is being used as a mathematical post-processing tool to compute
+    a set of dominant modes out of a given flow (or dataset) each of which is
+    associated with a constant frequency and decay rate, so that dynamic
+    features of a given flow may become interpretable, and tractable.
+    Among other Dynamic Mode Decomposition (DMD) variants, STDMD is presumed
+    to provide the general DMD method capabilities alongside economised and
+    feasible memory and CPU usage.
+
+    The code implementation corresponds to Figs. 15-16 of the first citation
+    below, more broadly to Section 2.4.
+
+    References:
+    \verbatim
+    STDMD and mode-sorting algorithms (tags:K, HRDC, KZ, HWR):
+        Kiewat, M. (2019).
+        Streaming modal decomposition approaches for vehicle aerodynamics.
+        PhD thesis. Munich: Technical University of Munich.
+        URL:mediatum.ub.tum.de/doc/1482652/1482652.pdf
+
+        Hemati, M. S., Rowley, C. W., Deem, E. A., & Cattafesta, L. N. (2017).
+        De-biasing the dynamic mode decomposition for applied Koopman
+        spectral analysis of noisy datasets.
+        Theoretical and Computational Fluid Dynamics, 31(4), 349-368.
+        DOI:10.1007/s00162-017-0432-2
+
+        Kou, J., & Zhang, W. (2017).
+        An improved criterion to select dominant modes from dynamic mode
+        decomposition.
+        European Journal of Mechanics-B/Fluids, 62, 109-129.
+        DOI:10.1016/j.euromechflu.2016.11.015
+
+        Hemati, M. S., Williams, M. O., & Rowley, C. W. (2014).
+        Dynamic mode decomposition for large and streaming datasets.
+        Physics of Fluids, 26(11), 111701.
+        DOI:10.1063/1.4901016
+
+    Parallel classical Gram-Schmidt process (tag:Ka):
+        Katagiri, T. (2003).
+        Performance evaluation of parallel Gram-Schmidt re-orthogonalization
+        methods.
+        In: Palma J. M. L. M., Sousa A. A., Dongarra J., Hernández V. (eds)
+        High Performance Computing for Computational Science — VECPAR 2002.
+        Lecture Notes in Computer Science, vol 2565, p. 302-314.
+        Berlin, Heidelberg: Springer.
+        DOI:10.1007/3-540-36569-9_19
+
+    Parallel direct tall-skinny QR decomposition (tags:BGD, DGHL):
+        Benson, A. R., Gleich, D. F., & Demmel, J. (2013).
+        Direct QR factorizations for tall-and-skinny matrices in MapReduce
+        architectures.
+        2013 IEEE International Conference on Big Data.
+        DOI:10.1109/bigdata.2013.6691583
+
+        Demmel, J., Grigori, L., Hoemmen, M., & Langou, J. (2012).
+        Communication-optimal parallel and sequential QR and LU
+        factorizations.
+        SIAM Journal on Scientific Computing, 34(1), A206-A239.
+        DOI:10.1137/080731992
+
+    DMD properties:
+        Brunton S. L. (2018).
+        Dynamic mode decomposition overview.
+        Seattle, Washington: University of Washington.
+        youtu.be/sQvrK8AGCAo (Retrieved:24-04-20)
+    \endverbatim
+
+    Operands:
+    \table
+      Operand      | Type              | Location
+      input        | {vol,surface}\<Type\>Field <!--
+               --> | $FOAM_CASE/\<time\>/\<inpField\>
+      output file  | dat                        <!--
+               --> | $FOAM_CASE/postProcessing/\<FO\>/\<time\>/\<file\>(s)
+      output field | volScalarField(s) | $FOAM_CASE/\<time\>/\<outField\>(s)
+    \endtable
+
+    where Type={Scalar,SphericalTensor,SymmTensor,Tensor}.
+
+    Output fields:
+    \verbatim
+      modeReal<modeIndex><field> | Real part of a mode field
+      modeImag<modeIndex><field> | Imaginary part of a mode field
+    \endverbatim
+
+    Output files:
+    \verbatim
+      uSTDMD.dat                 | Unfiltered STDMD output
+      STDMD.dat                  | Filtered STDMD output
+    \endverbatim
+
+    wherein for each mode, the following quantities are output into files:
+    \verbatim
+      Frequency
+      Magnitude
+      Amplitude (real part)
+      Amplitude (imaginary part)
+      Eigenvalue (real part)
+      Eigenvalue (imaginary part)
+    \endverbatim
+
+Note
+    Operations on boundary fields, e.g. \c wallShearStress, are currently not
+    available.
+
+Usage
+    Minimal example by using \c system/controlDict.functions:
+    \verbatim
+    STDMD1
+    {
+        // Mandatory entries (unmodifiable)
+        type                STDMD;
+        libs                (fieldFunctionObjects);
+        field               <inpField>;
+
+        // Conditional mandatory entries (unmodifiable)
+        // either of the below
+        stdmdInterval       5.5;
+        executeInterval     10;
+
+        // Optional entries (unmodifiable)
+        modeSorter          kiewat;
+        nModes              50;
+        maxRank             50;
+        nGramSchmidt        5;
+        fMin                0;
+        fMax                1000000000;
+
+        // Optional entries (run-time modifiable, yet not recommended)
+        testEigen           false;
+        dumpEigen           false;
+        minBasis            0.00000001;
+        minEVal             0.00000001;
+        absTol              0.001;
+        relTol              0.0001;
+        sortLimiter         500;
+
+        // Optional (inherited) entries
+        ...
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property     | Description                        | Type | Req'd | Dflt
+      type         | Type name: STDMD                   | word |  yes  | -
+      libs         | Library name: fieldFunctionObjects | word |  yes  | -
+      field        | Name of the operand field          | word |  yes  | -
+      stdmdInterval | STDMD time-step size [s]   | scalar| conditional | <!--
+                  --> executeInterval*(current time-step of the simulation)
+      modeSorter   | Mode-sorting algorithm variant     | word | no  | firstSnap
+      nModes | Number of output modes in input freq range | label | no | GREAT
+      maxRank  | Max columns in accumulated matrices      | label | no | GREAT
+      nGramSchmidt | Number of Gram-Schmidt iterations    | label | no | 5
+      fMin         | Min (non-negative) output frequency  | label | no | 0
+      fMax         | Max output frequency                 | label | no | GREAT
+      testEigen    | Flag to verify eigendecompositions   | bool  | no | false
+      dumpEigen    | Flag to log operands of eigendecomps | bool  | no | false
+      minBasis | Orthogonal basis expansion threshold     | scalar| no | 1e-8
+      minEVal  | Min EVal for below EVals are omitted     | scalar| no | 1e-8
+      absTol   | Min abs tol in eigendecomposition tests  | scalar| no | 1e-4
+      relTol   | Relative tol in eigendecomposition tests | scalar| no | 1e-6
+      sortLimiter | Maximum allowable magnitude for <!--
+               -->  mag(eigenvalue)*(number of STDMD steps) to be used in <!--
+               -->  modeSorter={kiewat,kouZhang} to avoid overflow errors <!--
+               -->                                        | scalar | no | 500.0
+    \endtable
+
+    Options for the \c modeSorter entry:
+    \verbatim
+      kiewat    | Modified weighted amplitude scaling method
+      kouZhang  | Weighted amplitude scaling method
+      firstSnap | Method of first snapshot amplitude magnitude
+    \endverbatim
+
+    The inherited entries are elaborated in:
+     - \link functionObject.H \endlink
+     - \link writeFile.H \endlink
+
+    Usage by \c postProcess utility is not available.
+
+Note
+    - To specify the STDMD time-step size (not necessarily equal to the
+    time step of the simulation), entries of either \c stdmdInterval or
+    \c executeInterval must be available (highest to lowest precedence). While
+    \c stdmdInterval allows users to directly specify the STDMD time-step size
+    in seconds, in absence of \c stdmdInterval, for convenience,
+    \c executeInterval allows users to compute the STDMD time-step internally
+    by multiplying itself with the current time-step size of the simulation.
+    - Limitation: Restart is currently not available since intermediate writing
+    of STDMD matrices are not supported.
+    - Limitation: Non-physical input (e.g. full-zero fields) may upset STDMD.
+    - Warning: DMD is an active research area at the time of writing; therefore,
+    there would be cases whereat oddities might be encountered.
+    - Warning: This STDMD release is the \b beta release; therefore,
+    small-to-medium changes in input/output interfaces and internal structures
+    should be expected in the next versions.
+    - Warning: In \c modeSorter algorithms of \c kiewat and \c kouZhang, there
+    exists a function of \f$power(x,y)\f$ where \c x is the magnitude of an
+    eigenvalue, and \c y is the number of STDMD snapshots. This function poses
+    a risk of overflow errors since, for example, if \c x ends up above 1.5 with
+    500 snapshots or more, this function automatically throws the floating point
+    error meaning overflow. Therefore, the domain-range of this function is
+    heuristically constrained by the optional entry \c sortLimiter which skips
+    the evaluation of this function for a given
+    mag(eigenvalue)*(no. of STDMD steps), i.e. x*y, whose magnitude is larger
+    than \c sortLimiter.
+
+See also
+    - Foam::functionObjects::fvMeshFunctionObject
+    - Foam::functionObjects::writeFile
+    - ExtendedCodeGuide::functionObjects::field::STDMD
+
+SourceFiles
+    STDMD.C
+    STDMDTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef functionObjects_STDMD_H
+#define functionObjects_STDMD_H
+
+#include "fvMeshFunctionObject.H"
+#include "writeFile.H"
+#include "EigenMatrix.H"
+#include "QRMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+
+/*---------------------------------------------------------------------------*\
+                              Class STDMD Declaration
+\*---------------------------------------------------------------------------*/
+
+class STDMD
+:
+    public fvMeshFunctionObject,
+    public writeFile
+{
+
+    typedef SquareMatrix<scalar> SMatrix;
+    typedef SquareMatrix<complex> SCMatrix;
+    typedef RectangularMatrix<scalar> RMatrix;
+    typedef RectangularMatrix<complex> RCMatrix;
+
+
+    // Private Enumerations
+
+        //- Options for the mode-sorting algorithm
+        enum modeSorterType
+        {
+            KIEWAT,          //!< "Modified weighted amplitude scaling method"
+            KOU_ZHANG,       //!< "Weighted amplitude scaling method"
+            FIRST_SNAPSHOT   //!< "Method of first snapshot amplitude magnitude"
+        };
+
+        //- Names for modeSorterType
+        static const Enum<modeSorterType> modeSorterTypeNames;
+
+
+    // Private Data
+
+        //- Mode-sorting algorithm (default=modeSorterType::FIRST_SNAPSHOT)
+        const enum modeSorterType modeSorter_;
+
+        //- Name of the operand volume or surface field
+        const word fieldName_;
+
+        //- Flag: First execution-step initialisation
+        bool initialised_;
+
+        //- Flag: Verify eigendecompositions by using theoretical relations
+        bool testEigen_;
+
+        //- Flag: Write operands of eigendecompositions to log
+        //  To activate it, testEigen=true
+        bool dumpEigen_;
+
+        //- Number of output modes within input frequency range
+        //- starting from the most energetic mode
+        const label nModes_;
+
+        //- Maximum allowable matrix column for Qz_ and Gz_
+        //  Qz_ is assumed to always have full-rank, thus Qz_.n() = rank
+        const label maxRank_;
+
+        //- Number of maximum iterations of the classical Gram-Schmidt procedure
+        const label nGramSchmidt_;
+
+        //- Min frequency: Output only entries corresponding to freqs >= fMin
+        const label fMin_;
+
+        //- Max frequency: Output only entries corresponding to freqs <= fMax
+        const label fMax_;
+
+        //- Number of base components of input field, e.g. 3 for vector
+        label nComps_;
+
+        //- Number of elements in a snapshot
+        //  A snapshot is an input dataset to be processed per execution-step
+        label nSnap_;
+
+        //- Current execution-step index of STDMD,
+        //- not necessarily that of the simulation
+        label currIndex_;
+
+        //- Maximum allowable magnitude for mag(eigenvalue)*(number of
+        //- STDMD steps) to be used in modeSorter={kiewat,kouZhang} to avoid
+        //- overflow errors
+        scalar sortLimiter_;
+
+        //- Min absolute tolerance being used in eigendecomposition tests
+        scalar absTol_;
+
+        //- Relative tolerance being used in eigendecomposition tests
+        scalar relTol_;
+
+        //- Min value to execute orthogonal basis expansion of Qz_ and Gz_
+        scalar minBasis_;
+
+        //- STDMD time-step size that equals to
+        //- (executeInterval of STDMD)*(deltaT of simulation) [s]
+        scalar dt_;
+
+        //- Min EVal magnitude threshold below where EVals are omitted
+        scalar minMagEVal_;
+
+        //- L2-norm of column vector z_
+        scalar zNorm_;
+
+        //- L2-norm of column vector ez_
+        scalar ezNorm_;
+
+        //- Augmented snapshot matrix (effectively a column vector) (K:Eq. 60)
+        //  Upper half z_ = current-time snapshot slot
+        //  Lower half z_ = previous-time snapshot slot
+        RMatrix z_;
+
+        //- Working copy of the augmented snapshot matrix z_
+        //- being used in the classical Gram-Schmidt process
+        RMatrix ez_;
+
+        //- First-processed snapshot required by the mode-sorting
+        //- algorithms at the final output computations (K:p. 43)
+        RMatrix X1_;
+
+        //- Accumulated-in-time unitary similarity matrix produced by the
+        //- orthonormal decomposition of the augmented snapshot matrix z_
+        //  (K:Eq. 60)
+        RMatrix Qz_;
+
+        //- Accumulated-in-time (squared) upper triangular matrix produced by
+        //- the orthonormal decomposition of the augmented snapshot matrix z_
+        //  (K:Eq. 64)
+        SMatrix Gz_;
+
+        //- Upper half of Qz_
+        RMatrix QzUH_;
+
+        //- Lower half of Qz_
+        RMatrix QzLH_;
+
+        //- Moore-Penrose pseudo-inverse of R produced by
+        //- the QR decomposition of the last time-step QzUH_
+        RMatrix RxInv_;
+
+        //- Projected STDMD operator (K:Eq. 78)
+        SMatrix Ap_;
+
+        //- Output eigenvectors
+        RCMatrix oEVecs_;
+
+        //- Output eigenvalues
+        List<complex> oEVals_;
+
+        //- Output amplitudes
+        List<complex> oAmps_;
+
+        //- Output (non-negative) frequencies
+        List<scalar> oFreqs_;
+
+        //- Indices of oFreqs_ where freqs are
+        //- non-negative and within [fMin_, fMax_]
+        DynamicList<label> iFreqs_;
+
+        //- Output (descending) magnitudes of (complex) amplitudes
+        List<scalar> oMags_;
+
+        //- Indices of oMags_
+        List<label> iMags_;
+
+
+    // Private Member Functions
+
+        // Process
+
+        //- Move previous-time snapshot into previous-time slot in z_
+        //- and copy new current-time snapshot into current-time slot in z_
+        template<class GeoFieldType>
+        bool getSnapshot();
+
+        //- Get the input field type to be processed by snapshot()
+        template<class Type>
+        bool getSnapshotType();
+
+        //- Get the number of base components of input field
+        template<class GeoFieldType>
+        bool getComps();
+
+        //- Return (parallel) L2-norm of a given column vector
+        scalar parnorm(const RMatrix& colVector) const;
+
+        //- Move the current-time snapshot into the previous-time snapshot in z_
+        //- and copy the new field into the current-time snapshot
+        void snapshot();
+
+        //- Initialise all working matrices at the first execution-step
+        void init();
+
+        //- Initialise Qz_, Gz_ (both require the first two snapshots) and X1_
+        void initBasis();
+
+        //- Execute (parallel) classical Gram-Schmidt
+        //- process to orthonormalise ez_ (Ka:Fig. 5)
+        void GramSchmidt();
+
+        //- Expand orthonormal bases Qz_ and Gz_ by stacking a column
+        //- (ez_/ezNorm_) to Qz_, and a row (Zero) and column (Zero)
+        //- to Gz_ if (minBasis_ < ezNorm_/zNorm_)
+        void expandBasis();
+
+        //- Update Gz_ before the potential orthonormal basis compression
+        void updateGz();
+
+        //- Compress orthonormal basis for Qz_ and Gz_ if (maxRank_ < Qz_.n())
+        void compressBasis();
+
+
+        // Postprocess
+
+        //- Return a new List containing elems of List at indices
+        template<class Type>
+        void filterIndexed
+        (
+            List<Type>& lst,
+            const UList<label>& indices
+        );
+
+        //- Return a new Matrix containing columns of Matrix at indices
+        template<class MatrixType>
+        void filterIndexed
+        (
+            MatrixType& lst,
+            const UList<label>& indices
+        );
+
+        //- Compute global Ap (K:Eq. 78)
+        void calcAp();
+
+        //- Compute eigenvalues and eigenvectors
+        void calcEigen();
+
+        //- Weak-type floating-point comparison
+        //  bit.ly/2Trdbgs (Eq. 2), and PEP-485
+        bool close
+        (
+            const scalar s1,
+            const scalar s2,
+            const scalar absTol = 0,   //<! comparisons near zero
+            const scalar relTol = 1e-8 //<! e.g. vals the same within 8 decimals
+        ) const;
+
+        //- Test real/complex eigenvalues by using
+        //- the theoretical relation: (sum(eigenvalues) - trace(A) ~ 0)
+        void testEigenvalues
+        (
+            const SquareMatrix<scalar>& A,
+            const DiagonalMatrix<scalar>& EValsRe
+        ) const;
+
+        //- Test real eigenvectors by using the theoretical relation:
+        //- ((A & EVec - EVal*EVec).norm() ~ 0)
+        void testEigenvectors
+        (
+            const SquareMatrix<scalar>& A,
+            const DiagonalMatrix<scalar>& EValsRe,
+            const SquareMatrix<scalar>& EVecs
+        ) const;
+
+        //- Test complex eigenvectors by using the theoretical relation:
+        //- ((A & EVec - EVal*EVec).norm() ~ 0)
+        void testEigenvectors
+        (
+            const SquareMatrix<scalar>& A,
+            const List<complex>& EVals,
+            const RectangularMatrix<complex>& EVecs
+        ) const;
+
+        //- Remove any eigenvalues whose magnitude is smaller than
+        //- minMagEVal_ while keeping the order of elements the same
+        void filterEVals();
+
+        //- Compute and filter frequencies (K:Eq. 81)
+        void calcFreqs();
+
+        //- Compute frequency indices
+        //  Locate indices where oFreqs_ are
+        //  in [fMin_, fMax_], and store them in iFreqs_ indices
+        void calcFreqI();
+
+        //- Compute amplitudes
+        void calcAmps();
+
+        //- Compute magnitudes
+        //  Includes advanced sorting methods using
+        //  the chosen weighted amplitude scaling method
+        void calcMags();
+
+        //- Compute the ordered magnitude indices
+        void calcMagI();
+
+        //- Compute modes
+        void calcModes();
+
+        //- Eigenvalue weighted amplitude scaling (KZ:Eq. 33)
+        //- Modified eigenvalue weighted amplitude scaling (K)
+        scalar sorter
+        (
+            const List<scalar>& weight,
+            const complex& amplitude,
+            const complex& eval,
+            const scalar modeNorm
+        ) const;
+
+        //- Output file header information
+        virtual void writeFileHeader(Ostream& os) const;
+
+        //- Filter objects according to iFreqs_ and iMags_
+        void filterOutput();
+
+        //- Write unfiltered/filtered data
+        void writeOutput(OFstream& os) const;
+
+        //- Compute STDMD output
+        void calcOutput();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("STDMD");
+
+
+    // Constructors
+
+        //- No default construct
+        STDMD() = delete;
+
+        //- Construct from Time and dictionary
+        STDMD
+        (
+            const word& name,
+            const Time& runTime,
+            const dictionary& dict
+        );
+
+        //- No copy construct
+        STDMD(const STDMD&) = delete;
+
+        //- No copy assignment
+        void operator=(const STDMD&) = delete;
+
+
+    //- Destructor
+    virtual ~STDMD() = default;
+
+
+    // Member Functions
+
+        //- Read STDMD settings
+        virtual bool read(const dictionary&);
+
+        //- Execute STDMD
+        virtual bool execute();
+
+        //- Write STDMD output
+        virtual bool write();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace functionObjects
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "STDMDTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/STDMD/STDMDTemplates.C b/src/functionObjects/field/STDMD/STDMDTemplates.C
new file mode 100644
index 0000000000..e32bb2ddc3
--- /dev/null
+++ b/src/functionObjects/field/STDMD/STDMDTemplates.C
@@ -0,0 +1,141 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "volFields.H"
+#include "surfaceFields.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class GeoFieldType>
+bool Foam::functionObjects::STDMD::getSnapshot()
+{
+    if (!initialised_)
+    {
+        init();
+    }
+
+    // Move previous-time snapshot into previous-time slot in z_
+    // Effectively moves the lower half of z_ to its upper half
+    std::rotate(z_.begin(), z_.begin() + nSnap_, z_.end());
+
+    // Copy new current-time snapshot into current-time slot in z_
+    // Effectively copies the new field elements into the lower half of z_
+    const GeoFieldType& Field = lookupObject<GeoFieldType>(fieldName_);
+    const label nField = Field.size();
+
+    for (direction dir = 0; dir < nComps_; ++dir)
+    {
+        z_.subColumn(0, nSnap_ + dir*nField, nField) = Field.component(dir);
+    }
+
+    return true;
+}
+
+
+template<class Type>
+bool Foam::functionObjects::STDMD::getSnapshotType()
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
+    typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
+
+    if (foundObject<VolFieldType>(fieldName_))
+    {
+        return getSnapshot<VolFieldType>();
+    }
+    else if (foundObject<SurfaceFieldType>(fieldName_))
+    {
+        return getSnapshot<SurfaceFieldType>();
+    }
+
+    return false;
+}
+
+
+template<class Type>
+bool Foam::functionObjects::STDMD::getComps()
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
+    typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
+
+    if (foundObject<VolFieldType>(fieldName_))
+    {
+        nComps_ = pTraits<typename VolFieldType::value_type>::nComponents;
+        return true;
+    }
+    else if (foundObject<SurfaceFieldType>(fieldName_))
+    {
+        nComps_ = pTraits<typename SurfaceFieldType::value_type>::nComponents;
+        return true;
+    }
+
+    return false;
+}
+
+
+template<class Type>
+void Foam::functionObjects::STDMD::filterIndexed
+(
+    List<Type>& lst,
+    const UList<label>& indices
+)
+{
+    // Elems within [a, b]
+    List<Type> lstWithin(indices.size());
+
+    // Copy if frequency of elem is within [a, b]
+    label j = 0;
+    for (const auto& i : indices)
+    {
+        lstWithin[j] = lst[i];
+        ++j;
+    }
+    lst.transfer(lstWithin);
+}
+
+
+template<class MatrixType>
+void Foam::functionObjects::STDMD::filterIndexed
+(
+    MatrixType& mat,
+    const UList<label>& indices
+)
+{
+    // Elems within [a, b]
+    MatrixType matWithin(labelPair(mat.m(), indices.size()));
+
+    // Copy if frequency of elem is within [a, b]
+    label j = 0;
+    for (const auto& i : indices)
+    {
+        matWithin.subColumn(j) = mat.subColumn(i);
+        ++j;
+    }
+    mat.transfer(matWithin);
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/0.orig/U b/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/0.orig/U
new file mode 100644
index 0000000000..6e2a24d9fa
--- /dev/null
+++ b/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/0.orig/U
@@ -0,0 +1,51 @@
+/*--------------------------------*- 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       volVectorField;
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0.012984 0 0);
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+
+    outlet
+    {
+        type            zeroGradient;
+    }
+
+    topAndBottom
+    {
+        type            symmetry;
+    }
+
+    "left|right"
+    {
+        type            empty;
+    }
+
+    cylinder
+    {
+        type            noSlip;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/0.orig/p b/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/0.orig/p
new file mode 100644
index 0000000000..1b6348bf3d
--- /dev/null
+++ b/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/0.orig/p
@@ -0,0 +1,51 @@
+/*--------------------------------*- 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       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            zeroGradient;
+    }
+
+    outlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+
+    topAndBottom
+    {
+        type            symmetry;
+    }
+
+    "left|right"
+    {
+        type            empty;
+    }
+
+    cylinder
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/Allclean b/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/Allclean
new file mode 100755
index 0000000000..41459e3ce5
--- /dev/null
+++ b/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/Allclean
@@ -0,0 +1,12 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/CleanFunctions      # Tutorial clean functions
+#------------------------------------------------------------------------------
+
+cleanCase0
+
+rm -rf dynamicCode
+rm -rf constant/coarseMesh
+rm -f system/coarseMesh/decomposeParDict
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/Allrun b/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/Allrun
new file mode 100755
index 0000000000..0307797569
--- /dev/null
+++ b/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/Allrun
@@ -0,0 +1,25 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
+#------------------------------------------------------------------------------
+
+./Allrun.pre
+
+runApplication decomposePar
+
+runParallel renumberMesh -overwrite
+
+cp -f system/decomposeParDict system/coarseMesh
+
+runApplication -s coarseMesh decomposePar -region coarseMesh
+
+runParallel -s coarseMesh renumberMesh -overwrite -region coarseMesh
+
+runParallel $(getApplication)
+
+runParallel -s main redistributePar -reconstruct -latestTime
+
+runParallel -s coarseMesh redistributePar -reconstruct \
+    -region coarseMesh -time '50,100,200'
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/Allrun.pre b/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/Allrun.pre
new file mode 100755
index 0000000000..bc21f038db
--- /dev/null
+++ b/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/Allrun.pre
@@ -0,0 +1,20 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
+#------------------------------------------------------------------------------
+
+runApplication -s coarseMesh blockMesh -dict system/blockMeshDict.coarse
+
+runApplication snappyHexMesh -overwrite
+
+mkdir -p constant/coarseMesh
+
+mv -f constant/polyMesh constant/coarseMesh
+
+runApplication -s main blockMesh -dict system/blockMeshDict.main
+
+runApplication mirrorMesh -overwrite
+
+restore0Dir
+
+#------------------------------------------------------------------------------
\ No newline at end of file
diff --git a/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/constant/transportProperties b/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/constant/transportProperties
new file mode 100644
index 0000000000..ae22a9bc6d
--- /dev/null
+++ b/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/constant/transportProperties
@@ -0,0 +1,23 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+transportModel  Newtonian;
+
+nu              1.558e-5;
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/constant/turbulenceProperties b/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/constant/turbulenceProperties
new file mode 100644
index 0000000000..e0b42c8d7d
--- /dev/null
+++ b/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/constant/turbulenceProperties
@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType laminar;
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/blockMeshDict.coarse b/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/blockMeshDict.coarse
new file mode 100644
index 0000000000..794d0060d2
--- /dev/null
+++ b/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/blockMeshDict.coarse
@@ -0,0 +1,74 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+scale    1;
+
+x0       -0.16;
+x1       0.8;
+y0       -0.24;
+y1       0.24;
+z0       -0.00375;
+z1       0.00375;
+
+nx       90;
+ny       60;
+nz       1;
+
+vertices
+(
+    ($x0 $y0 $z0)
+    ($x1 $y0 $z0)
+    ($x1 $y1 $z0)
+    ($x0 $y1 $z0)
+    ($x0 $y0 $z1)
+    ($x1 $y0 $z1)
+    ($x1 $y1 $z1)
+    ($x0 $y1 $z1)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) ($nx $ny $nz) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    allPatches
+    {
+        type     empty;
+        inGroups (allPatches);
+        faces
+        (
+            (3 7 6 2)
+            (1 5 4 0)
+            (0 4 7 3)
+            (2 6 5 1)
+            (0 3 2 1)
+            (4 5 6 7)
+        );
+    }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //
+
diff --git a/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/blockMeshDict.main b/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/blockMeshDict.main
new file mode 100644
index 0000000000..1ea36be199
--- /dev/null
+++ b/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/blockMeshDict.main
@@ -0,0 +1,277 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+scale     1;
+
+R         0.06; // Cylinder radius
+x0        #eval{ 0.5*$R };
+x2        #eval{ 2.38732*$R };
+x3        #eval{ 10.0*$R };
+xOutlet   #eval{ 18.6667*$R };
+xInlet    #eval{ -10.125*$R };
+
+RsinPi8   #eval{ $R*sin(0.125*pi()) };
+RsinPi8n  #eval{ -$RsinPi8 };
+RcosPi8   #eval{ $R*cos(0.125*pi()) };
+RcosPi8n  #eval{ -$RcosPi8 };
+RsinPi4   #eval{ $R*sin(0.25*pi()) };
+
+x2sinPi8  #eval{ $x2*sin(0.125*pi()) };
+x2sinPi8n #eval{ -$x2sinPi8 };
+x2cosPi8  #eval{ $x2*cos(0.125*pi()) };
+x2cosPi8n #eval{ -$x2cosPi8 };
+x2sinPi4  #eval{ $x2*sin(0.25*pi()) };
+
+z0        -0.0075;
+z1        0.0075;
+nz        1;
+
+
+vertices #codeStream
+{
+    codeInclude
+    #{
+        #include "pointField.H"
+    #};
+
+    code
+    #{
+        pointField points(19);
+        points[0] = point($R, 0, $z0);
+        points[1] = point($x2, 0, $z0);
+        points[2] = point($x3, 0, $z0);
+        points[3] = point($x3, $x2sinPi4, $z0);
+        points[4] = point($x2sinPi4, $x2sinPi4, $z0);
+        points[5] = point($RsinPi4, $RsinPi4, $z0);
+        points[6] = point($x3, $x3, $z0);
+        points[7] = point($x2sinPi4, $x3, $z0);
+
+        // Mirror +x points to -x side
+        points[11] = point(-points[0].x(), points[0].y(), points[0].z());
+        points[12] = point(-points[1].x(), points[1].y(), points[1].z());
+        points[13] = point(-points[2].x(), points[2].y(), points[2].z());
+        points[14] = point(-points[3].x(), points[3].y(), points[3].z());
+        points[15] = point(-points[4].x(), points[4].y(), points[4].z());
+        points[16] = point(-points[5].x(), points[5].y(), points[5].z());
+        points[17] = point(-points[6].x(), points[6].y(), points[6].z());
+        points[18] = point(-points[7].x(), points[7].y(), points[7].z());
+
+        // Points on the y-z-plane
+        points[8] = point(0, $x3, $z0);
+        points[9] = point(0, $x2, $z0);
+        points[10] = point(0, $R, $z0);
+
+        // Mirror -z points to +z side
+        label sz = points.size();
+        points.setSize(2*sz);
+        for (label i = 0; i < sz; ++i)
+        {
+            const point& pt = points[i];
+            points[i + sz] = point(pt.x(), pt.y(), $z1);
+        }
+
+        // Add an inner cylinder
+        sz = points.size();
+        label nAdd = 6;
+        points.setSize(sz + nAdd);
+
+        // Points within the inner cylinder
+        points[sz] = point(0, 0, $z0);
+        points[sz + 1] = point(0, $x0, $z0);
+        points[sz + 2] = point($x0, $x0, $z0);
+        points[sz + 3] = point($x0, 0, $z0);
+
+        // Mirror points from +x side to -x side
+        points[sz + 4] =
+            point(-points[sz + 2].x(), points[sz + 2].y(), points[sz + 2].z());
+        points[sz + 5] =
+            point(-points[sz + 3].x(), points[sz + 3].y(), points[sz + 3].z());
+
+        // Mirror -z points to +z side
+        sz = points.size();
+        points.setSize(sz + nAdd);
+        for (label i = 0; i < nAdd; ++i)
+        {
+            const point& pt = points[i+sz-nAdd];
+            points[i+sz] = point(pt.x(), pt.y(), $z1);
+        }
+
+        // Add downstream and upstream blocks
+        sz = points.size();
+        nAdd = 6;
+        points.setSize(sz + nAdd);
+
+        // Points on outlet
+        points[sz] = point($xOutlet, 0, $z0);
+        points[sz + 1] = point($xOutlet, $x3, $z0);
+        points[sz + 4] = point($xOutlet, $x2sinPi4, $z0);
+
+        // Points on inlet
+        points[sz + 2] = point($xInlet, 0, $z0);
+        points[sz + 3] = point($xInlet, $x3, $z0);
+        points[sz + 5] = point($xInlet, $x2sinPi4, $z0);
+
+        // Mirror -z points to +z side
+        sz = points.size();
+        points.setSize(sz + nAdd);
+        for (label i = 0; i < nAdd; ++i)
+        {
+            const point& pt = points[i + sz - nAdd];
+            points[i + sz] = point(pt.x(), pt.y(), $z1);
+        }
+
+        os  << points;
+    #};
+};
+
+
+blocks
+(
+    hex ( 5  4  9 10 24 23 28 29) (16 15 $nz) simpleGrading (1.6 1 1)
+    hex ( 0  1  4  5 19 20 23 24) (16 15 $nz) simpleGrading (1.6 1 1)
+    hex ( 1  2  3  4 20 21 22 23) (61 15 $nz) simpleGrading (1 1 1)
+    hex ( 4  3  6  7 23 22 25 26) (61 61 $nz) simpleGrading (1 1 1)
+    hex ( 9  4  7  8 28 23 26 27) (15 61 $nz) simpleGrading (1 1 1)
+    hex (15 16 10  9 34 35 29 28) (16 15 $nz) simpleGrading (0.625 1 1)
+    hex (12 11 16 15 31 30 35 34) (16 15 $nz) simpleGrading (0.625 1 1)
+    hex (13 12 15 14 32 31 34 33) (61 15 $nz) simpleGrading (1 1 1)
+    hex (14 15 18 17 33 34 37 36) (61 61 $nz) simpleGrading (1 1 1)
+    hex (15  9  8 18 34 28 27 37) (15 61 $nz) simpleGrading (1 1 1)
+    hex (2 50 54 3 21 56 60 22)  (69 15 $nz) simpleGrading (1 1 1) // downstream
+    hex (3 54 51 6 22 60 57 25)  (69 61 $nz) simpleGrading (1 1 1)
+    hex (52 13 14 55 58 32 33 61)  (1 15 $nz) simpleGrading (1 1 1) // upstream
+    hex (55 14 17 53 61 33 36 59)  (1 61 $nz) simpleGrading (1 1 1)
+);
+
+edges
+(
+    arc  0  5 ($RcosPi8   $RsinPi8  $z0)
+    arc  5 10 ($RsinPi8   $RcosPi8  $z0)
+    arc  1  4 ($x2cosPi8  $x2sinPi8 $z0)
+    arc  4  9 ($x2sinPi8  $x2cosPi8 $z0)
+    arc 19 24 ($RcosPi8   $RsinPi8  $z1)
+    arc 24 29 ($RsinPi8   $RcosPi8  $z1)
+    arc 20 23 ($x2cosPi8  $x2sinPi8 $z1)
+    arc 23 28 ($x2sinPi8  $x2cosPi8 $z1)
+    arc 11 16 ($RcosPi8n  $RsinPi8  $z0)
+    arc 16 10 ($RsinPi8n  $RcosPi8  $z0)
+    arc 12 15 ($x2cosPi8n $x2sinPi8 $z0)
+    arc 15  9 ($x2sinPi8n $x2cosPi8 $z0)
+    arc 30 35 ($RcosPi8n  $RsinPi8  $z1)
+    arc 35 29 ($RsinPi8n  $RcosPi8  $z1)
+    arc 31 34 ($x2cosPi8n $x2sinPi8 $z1)
+    arc 34 28 ($x2sinPi8n $x2cosPi8 $z1)
+);
+
+boundary
+(
+    inlet
+    {
+        type patch;
+        faces
+        (
+            (53 59 61 55)
+            (55 61 58 52)
+        );
+    }
+    outlet
+    {
+        type patch;
+        faces
+        (
+            (60 57 51 54)
+            (56 60 54 50)
+        );
+    }
+    topAndBottom
+    {
+        type symmetry;
+        faces
+        (
+            (59 53 17 36)
+            (36 17 18 37)
+            (37 18  8 27)
+            (27  8  7 26)
+            (26  7  6 25)
+            (25  6 51 57)
+            (52 58 32 13)
+            (13 32 31 12)
+            (12 31 30 11)
+            (0 19 20 1)
+            (1 20 21 2)
+            (2 21 56 50)
+        );
+    }
+    left
+    {
+        type empty;
+        faces
+        (
+            (55 52 13 14)
+            (53 55 14 17)
+            (13 12 15 14)
+            (14 15 18 17)
+            (15 9 8 18)
+            (12 11 16 15)
+            (15 16 10 9)
+            (9 10 5 4)
+            (5 0 1 4)
+            (9 4 7 8)
+            (4 3 6 7)
+            (1 2 3 4)
+            (3 2 50 54)
+            (6 3 54 51)
+        );
+    }
+    right
+    {
+        type empty;
+        faces
+        (
+            (59 36 33 61)
+            (61 33 32 58)
+            (36 37 34 33)
+            (33 34 31 32)
+            (37 27 28 34)
+            (27 26 23 28)
+            (26 25 22 23)
+            (23 22 21 20)
+            (25 57 60 22)
+            (22 60 56 21)
+            (24 23 20 19)
+            (28 23 24 29)
+            (34 28 29 35)
+            (31 34 35 30)
+        );
+    }
+    cylinder
+    {
+        type wall;
+        faces
+        (
+            (29 24 5 10)
+            (24 19 0 5)
+            (30 35 16 11)
+            (35 29 10 16)
+        );
+    }
+);
+
+mergePatchPairs
+();
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/coarseMesh/fvSchemes b/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/coarseMesh/fvSchemes
new file mode 100644
index 0000000000..172e4c34ad
--- /dev/null
+++ b/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/coarseMesh/fvSchemes
@@ -0,0 +1,47 @@
+/*--------------------------------*- 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      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+}
+
+gradSchemes
+{
+}
+
+divSchemes
+{
+}
+
+laplacianSchemes
+{
+}
+
+interpolationSchemes
+{
+}
+
+snGradSchemes
+{
+}
+
+fluxRequired
+{
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/coarseMesh/fvSolution b/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/coarseMesh/fvSolution
new file mode 100644
index 0000000000..133f1e3bf7
--- /dev/null
+++ b/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/coarseMesh/fvSolution
@@ -0,0 +1,27 @@
+/*--------------------------------*- 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      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+}
+
+PIMPLE
+{
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/controlDict b/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/controlDict
new file mode 100644
index 0000000000..5ef1852478
--- /dev/null
+++ b/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/controlDict
@@ -0,0 +1,229 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     pimpleFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         200;
+
+deltaT          0.05;
+
+writeControl    runTime;
+
+writeInterval   50;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  16;
+
+timeFormat      general;
+
+timePrecision   8;
+
+runTimeModifiable false;
+
+
+functions
+{
+    mapFields1
+    {
+        type            mapFields;
+        libs            (fieldFunctionObjects);
+        mapRegion       coarseMesh;
+        mapMethod       cellVolumeWeight;
+        consistent      no;
+        patchMap        ();
+        cuttingPatches  ();
+        fields          ( U p );
+
+        timeStart       10;
+        timeEnd         2000;
+
+        executeControl  timeStep;
+        executeInterval 10;
+
+        writeControl    timeStep;
+        writeInterval   50;
+    }
+
+    STDMD1U
+    {
+        // Mandatory entries
+        type                STDMD;
+        libs                (fieldFunctionObjects);
+        field               U;
+
+        // Optional entries
+        modeSorter          firstSnap;
+
+        // Optional (inherited) entries
+        region              coarseMesh;  // mapFields must be executed before.
+        timeStart           10;
+        executeControl      timeStep;
+        executeInterval     10;
+        writeControl        onEnd;
+    }
+
+    STDMD1p
+    {
+        // Mandatory entries
+        type                STDMD;
+        libs                (fieldFunctionObjects);
+        field               p;
+
+        // Optional entries
+        modeSorter          firstSnap;
+
+        // Optional (inherited) entries
+        region              coarseMesh;
+        timeStart           10;
+        executeControl      timeStep;
+        executeInterval     10;
+        writeControl        onEnd;
+    }
+
+    STDMD2U
+    {
+        // Mandatory entries
+        type                STDMD;
+        libs                (fieldFunctionObjects);
+        field               U;
+
+        // Optional entries
+        modeSorter          firstSnap;
+        nModes              2;
+        maxRank             2;
+        nGramSchmidt        5;
+        fMin                0;
+        fMax                1000000000;
+        testEigen           true;
+        dumpEigen           true;
+        minBasis            0.00000001;
+        minEVal             0.00000001;
+        absTol              0.001;
+        relTol              0.0001;
+
+        // Optional (inherited) entries
+        timeStart           50;
+        timeEnd             2000;
+        executeControl      timeStep;
+        executeInterval     10;
+        writeControl        onEnd;
+    }
+
+    STDMD3U
+    {
+        // Mandatory entries
+        type                STDMD;
+        libs                (fieldFunctionObjects);
+        field               U;
+
+        // Optional entries
+        modeSorter          firstSnap;
+
+        // Optional (inherited) entries
+        region              coarseMesh;
+        timeStart           10;
+        timeEnd             150;
+        executeControl      timeStep;
+        executeInterval     10;
+        writeControl        runTime;
+        writeInterval       50;
+    }
+
+    STDMD4U
+    {
+        // Mandatory entries
+        type                STDMD;
+        libs                (fieldFunctionObjects);
+        field               U;
+
+        // Optional entries
+        modeSorter          kiewat;
+
+        // Optional (inherited) entries
+        region              coarseMesh;
+        timeStart           10;
+        executeControl      timeStep;
+        executeInterval     10;
+        writeControl        onEnd;
+    }
+
+    STDMD4p
+    {
+        // Mandatory entries
+        type                STDMD;
+        libs                (fieldFunctionObjects);
+        field               p;
+
+        // Optional entries
+        modeSorter          kiewat;
+
+        // Optional (inherited) entries
+        region              coarseMesh;
+        timeStart           10;
+        executeControl      timeStep;
+        executeInterval     10;
+        writeControl        onEnd;
+    }
+
+    STDMD5U
+    {
+        // Mandatory entries
+        type                STDMD;
+        libs                (fieldFunctionObjects);
+        field               U;
+
+        // Optional entries
+        modeSorter          kouZhang;
+
+        // Optional (inherited) entries
+        region              coarseMesh;
+        timeStart           10;
+        executeControl      timeStep;
+        executeInterval     10;
+        writeControl        onEnd;
+    }
+
+    STDMD5p
+    {
+        // Mandatory entries
+        type                STDMD;
+        libs                (fieldFunctionObjects);
+        field               p;
+
+        // Optional entries
+        modeSorter          kouZhang;
+
+        // Optional (inherited) entries
+        region              coarseMesh;
+        timeStart           10;
+        executeControl      timeStep;
+        executeInterval     10;
+        writeControl        onEnd;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/decomposeParDict b/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/decomposeParDict
new file mode 100644
index 0000000000..5c92caa799
--- /dev/null
+++ b/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/decomposeParDict
@@ -0,0 +1,28 @@
+/*--------------------------------*- 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      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains  8;
+
+method              hierarchical;
+
+coeffs
+{
+    n               (8 1 1);
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/fvSchemes b/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/fvSchemes
new file mode 100644
index 0000000000..45794b927c
--- /dev/null
+++ b/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/fvSchemes
@@ -0,0 +1,56 @@
+/*--------------------------------*- 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      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    default         Gauss linear;
+    div((nuEff*dev(T(grad(U))))) Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+fluxRequired
+{
+    default         no;
+    p               ;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/fvSolution b/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/fvSolution
new file mode 100644
index 0000000000..e76425bfb1
--- /dev/null
+++ b/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/fvSolution
@@ -0,0 +1,56 @@
+/*--------------------------------*- 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      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    p
+    {
+        solver           GAMG;
+        smoother         DICGaussSeidel;
+        tolerance        1e-6;
+        relTol           0.01;
+    }
+
+    pFinal
+    {
+        $p;
+        relTol          0;
+    }
+
+    U
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-05;
+        relTol          0.1;
+    }
+
+    UFinal
+    {
+        $U;
+        relTol          0;
+    }
+}
+
+PIMPLE
+{
+    nNonOrthogonalCorrectors 0;
+    nCorrectors         2;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/mirrorMeshDict b/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/mirrorMeshDict
new file mode 100644
index 0000000000..7b48e324c7
--- /dev/null
+++ b/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/mirrorMeshDict
@@ -0,0 +1,27 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      mirrorMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+planeType           pointAndNormal;
+
+pointAndNormalDict
+{
+    point   (0 0 0);
+    normal  (0 -1 0);
+}
+
+planeTolerance      1e-3;
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/snappyHexMeshDict b/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/snappyHexMeshDict
new file mode 100644
index 0000000000..268920eab9
--- /dev/null
+++ b/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/snappyHexMeshDict
@@ -0,0 +1,83 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      snappyHexMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+castellatedMesh true;
+snap            false;
+addLayers       false;
+
+geometry
+{
+    cylinder
+    {
+        type    searchableCylinder;
+        point1  (0 0 -0.00375);
+        point2  (0 0 0.00375);
+        radius  0.06;
+    }
+}
+
+
+castellatedMeshControls
+{
+    maxLocalCells       100000;
+    maxGlobalCells      100000000;
+    minRefinementCells  10;
+    maxLoadUnbalance    0.10;
+    nCellsBetweenLevels 1;
+    resolveFeatureAngle 30;
+    allowFreeStandingZoneFaces true;
+    features
+    ();
+    refinementSurfaces
+    {
+        cylinder
+        {
+            level (0 0);
+            patchInfo
+            {
+                type     empty;
+                inGroups (allPatches);
+            }
+        }
+    };
+    refinementRegions{};
+    locationInMesh (-0.1 -0.2 0);
+}
+
+
+snapControls
+{
+    nSmoothPatch    3;
+    tolerance       2.0;
+    nSolveIter      100;
+    nRelaxIter      5;
+}
+
+
+addLayersControls
+{}
+
+
+meshQualityControls
+{
+    #includeEtc "caseDicts/meshQualityDict"
+}
+
+debug 0;
+mergeTolerance 1e-6;
+
+
+// ************************************************************************* //