From 730233cd15d9e7b34987cdc3fa87feee54145589 Mon Sep 17 00:00:00 2001
From: Kutalmis Bercin <kutalmis.bercin@esi-group.com>
Date: Thu, 5 Mar 2020 11:37:45 +0000
Subject: [PATCH] ENH: add new FO Streaming-Total Dynamic Mode Decomposition
 (STDMD)

    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.

    Please refer to the header file documentation for further details.

  ENH: add new STDMD tutorial, pimpleFoam/laminar/cylinder2D
---
 src/functionObjects/field/Make/files          |    2 +
 src/functionObjects/field/Make/options        |    2 +
 src/functionObjects/field/STDMD/STDMD.C       | 1280 +++++++++++++++++
 src/functionObjects/field/STDMD/STDMD.H       |  644 +++++++++
 .../field/STDMD/STDMDTemplates.C              |  141 ++
 .../pimpleFoam/laminar/cylinder2D/0.orig/U    |   51 +
 .../pimpleFoam/laminar/cylinder2D/0.orig/p    |   51 +
 .../pimpleFoam/laminar/cylinder2D/Allclean    |   12 +
 .../pimpleFoam/laminar/cylinder2D/Allrun      |   25 +
 .../pimpleFoam/laminar/cylinder2D/Allrun.pre  |   20 +
 .../cylinder2D/constant/transportProperties   |   23 +
 .../cylinder2D/constant/turbulenceProperties  |   21 +
 .../cylinder2D/system/blockMeshDict.coarse    |   74 +
 .../cylinder2D/system/blockMeshDict.main      |  277 ++++
 .../cylinder2D/system/coarseMesh/fvSchemes    |   47 +
 .../cylinder2D/system/coarseMesh/fvSolution   |   27 +
 .../laminar/cylinder2D/system/controlDict     |  229 +++
 .../cylinder2D/system/decomposeParDict        |   28 +
 .../laminar/cylinder2D/system/fvSchemes       |   56 +
 .../laminar/cylinder2D/system/fvSolution      |   56 +
 .../laminar/cylinder2D/system/mirrorMeshDict  |   27 +
 .../cylinder2D/system/snappyHexMeshDict       |   83 ++
 22 files changed, 3176 insertions(+)
 create mode 100644 src/functionObjects/field/STDMD/STDMD.C
 create mode 100644 src/functionObjects/field/STDMD/STDMD.H
 create mode 100644 src/functionObjects/field/STDMD/STDMDTemplates.C
 create mode 100644 tutorials/incompressible/pimpleFoam/laminar/cylinder2D/0.orig/U
 create mode 100644 tutorials/incompressible/pimpleFoam/laminar/cylinder2D/0.orig/p
 create mode 100755 tutorials/incompressible/pimpleFoam/laminar/cylinder2D/Allclean
 create mode 100755 tutorials/incompressible/pimpleFoam/laminar/cylinder2D/Allrun
 create mode 100755 tutorials/incompressible/pimpleFoam/laminar/cylinder2D/Allrun.pre
 create mode 100644 tutorials/incompressible/pimpleFoam/laminar/cylinder2D/constant/transportProperties
 create mode 100644 tutorials/incompressible/pimpleFoam/laminar/cylinder2D/constant/turbulenceProperties
 create mode 100644 tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/blockMeshDict.coarse
 create mode 100644 tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/blockMeshDict.main
 create mode 100644 tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/coarseMesh/fvSchemes
 create mode 100644 tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/coarseMesh/fvSolution
 create mode 100644 tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/controlDict
 create mode 100644 tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/decomposeParDict
 create mode 100644 tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/fvSchemes
 create mode 100644 tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/fvSolution
 create mode 100644 tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/mirrorMeshDict
 create mode 100644 tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/snappyHexMeshDict

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;
+
+
+// ************************************************************************* //