From 2a8e1c08653b5698de28404996a3e958fc409d2e Mon Sep 17 00:00:00 2001
From: Kutalmis Bercin <kutalmis.bercin@esi-group.com>
Date: Mon, 12 Apr 2021 12:00:47 +0100
Subject: [PATCH] ENH: DMD: refactor dynamic mode decomposition FO

ENH: DMD: enable operations on patch fields
ENH: DMD: enable postProcess utility
---
 src/functionObjects/field/DMD/DMD.C           |  234 +++
 src/functionObjects/field/DMD/DMD.H           |  278 ++++
 .../field/DMD/DMDModels/DMDModel/DMDModel.C   |   78 +
 .../field/DMD/DMDModels/DMDModel/DMDModel.H   |  203 +++
 .../DMD/DMDModels/DMDModel/DMDModelNew.C      |   58 +
 .../DMDModels/DMDModel/DMDModelTemplates.C    |   54 +
 .../field/DMD/DMDModels/derived/STDMD/STDMD.C |  988 +++++++++++++
 .../field/DMD/DMDModels/derived/STDMD/STDMD.H |  530 +++++++
 .../DMDModels/derived/STDMD/STDMDTemplates.C  |  236 +++
 .../STDMDTemplates.C => DMD/DMDTemplates.C}   |  161 +--
 src/functionObjects/field/Make/files          |    5 +-
 src/functionObjects/field/Make/options        |    2 +
 src/functionObjects/field/STDMD/STDMD.C       | 1280 -----------------
 src/functionObjects/field/STDMD/STDMD.H       |  645 ---------
 14 files changed, 2731 insertions(+), 2021 deletions(-)
 create mode 100644 src/functionObjects/field/DMD/DMD.C
 create mode 100644 src/functionObjects/field/DMD/DMD.H
 create mode 100644 src/functionObjects/field/DMD/DMDModels/DMDModel/DMDModel.C
 create mode 100644 src/functionObjects/field/DMD/DMDModels/DMDModel/DMDModel.H
 create mode 100644 src/functionObjects/field/DMD/DMDModels/DMDModel/DMDModelNew.C
 create mode 100644 src/functionObjects/field/DMD/DMDModels/DMDModel/DMDModelTemplates.C
 create mode 100644 src/functionObjects/field/DMD/DMDModels/derived/STDMD/STDMD.C
 create mode 100644 src/functionObjects/field/DMD/DMDModels/derived/STDMD/STDMD.H
 create mode 100644 src/functionObjects/field/DMD/DMDModels/derived/STDMD/STDMDTemplates.C
 rename src/functionObjects/field/{STDMD/STDMDTemplates.C => DMD/DMDTemplates.C} (51%)
 delete mode 100644 src/functionObjects/field/STDMD/STDMD.C
 delete mode 100644 src/functionObjects/field/STDMD/STDMD.H

diff --git a/src/functionObjects/field/DMD/DMD.C b/src/functionObjects/field/DMD/DMD.C
new file mode 100644
index 0000000000..e225d31cb4
--- /dev/null
+++ b/src/functionObjects/field/DMD/DMD.C
@@ -0,0 +1,234 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020-2021 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 "DMD.H"
+#include "DMDModel.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+    defineTypeNameAndDebug(DMD, 0);
+    addToRunTimeSelectionTable(functionObject, DMD, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::functionObjects::DMD::snapshot()
+{
+    bool processed = false;
+    processed = processed || getSnapshot<scalar>();
+    processed = processed || getSnapshot<vector>();
+    processed = processed || getSnapshot<sphericalTensor>();
+    processed = processed || getSnapshot<symmTensor>();
+    processed = processed || getSnapshot<tensor>();
+
+    if (!processed)
+    {
+        FatalErrorInFunction
+            << "    functionObjects::" << type() << " " << name() << ":"
+            << " cannot find required input field during snapshot loading: "
+            << fieldName_ << nl
+            << "    Do you execute required functionObjects"
+            << " before executing DMD, e.g. mapFields?"
+            << exit(FatalError);
+    }
+}
+
+
+void Foam::functionObjects::DMD::initialise()
+{
+    const label nComps = DMDModelPtr_->nComponents(fieldName_);
+
+    if (patch_.empty())
+    {
+        nSnap_ = nComps*mesh_.nCells();
+    }
+    else
+    {
+        const label patchi = mesh_.boundaryMesh().findPatchID(patch_);
+
+        if (patchi < 0)
+        {
+            FatalErrorInFunction
+                << "Cannot find patch " << patch_
+                << exit(FatalError);
+        }
+
+        nSnap_ = nComps*(mesh_.C().boundaryField()[patchi]).size();
+    }
+
+    const label nSnapTotal = returnReduce(nSnap_, sumOp<label>());
+
+    if (nSnapTotal <= 0)
+    {
+        FatalErrorInFunction
+            << "  # Zero-size input field = " << fieldName_ << " #"
+            << exit(FatalError);
+    }
+
+    if (nSnap_ > 0)
+    {
+        z_ = RMatrix(2*nSnap_, 1, Zero);
+    }
+    else
+    {
+        z_ = RMatrix(1, 1, Zero);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::functionObjects::DMD::DMD
+(
+    const word& name,
+    const Time& runTime,
+    const dictionary& dict
+)
+:
+    fvMeshFunctionObject(name, runTime, dict),
+    DMDModelPtr_(DMDModel::New(mesh_, name, dict)),
+    z_(),
+    fieldName_(dict.get<word>("field")),
+    patch_(dict.getOrDefault<word>("patch", word::null)),
+    nSnap_(0),
+    step_(0)
+{
+    // Check if computation time-step is varying
+    if (runTime.isAdjustTimeStep())
+    {
+        WarningInFunction
+            << "  # DMD: Available only for fixed time-step computations. #"
+            << endl;
+    }
+
+    // Check if mesh topology is changing
+    if (mesh_.topoChanging())
+    {
+        FatalErrorInFunction
+            << "  # DMD: Available only for non-changing mesh topology. #"
+            << exit(FatalError);
+    }
+
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::functionObjects::DMD::read(const dictionary& dict)
+{
+    Info<< type() << " " << name() << ":" << endl;
+
+    if (fvMeshFunctionObject::read(dict) && DMDModelPtr_->read(dict))
+    {
+        return true;
+    }
+
+    return false;
+}
+
+
+bool Foam::functionObjects::DMD::execute()
+{
+    Log << type() << " " << name() << " execute:" << endl;
+
+    snapshot();
+
+    if (step_ == 1)
+    {
+        DMDModelPtr_->initialise(z_);
+    }
+
+    if (step_ > 1)
+    {
+        DMDModelPtr_->update(z_);
+    }
+
+    ++step_;
+
+    Log << tab << "Execution index = " << step_ << " for field: " << fieldName_
+        << endl;
+
+    return true;
+}
+
+
+bool Foam::functionObjects::DMD::write()
+{
+    if (postProcess)
+    {
+        return true;
+    }
+
+    return end();
+}
+
+
+bool Foam::functionObjects::DMD::end()
+{
+    if (step_ == 0)
+    {
+        // Avoid double execution of write() when writeControl=onEnd
+        return true;
+    }
+
+    Log << type() << " " << name() << " write:" << endl;
+
+    if (step_ < 2)
+    {
+        WarningInFunction
+            << "  # DMD needs at least three snapshots to produce output #"
+            << nl
+            << "  # Only " << step_ + 1 << " snapshots are available #"
+            << nl
+            << "  # Skipping DMD output calculation and write #"
+            << endl;
+
+        return false;
+    }
+
+    z_.clear();
+
+    DMDModelPtr_->fit();
+
+    mesh_.time().printExecutionTime(Info);
+
+    // Restart the incremental orthonormal basis update
+    step_ = 0;
+
+    return true;
+}
+
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/DMD/DMD.H b/src/functionObjects/field/DMD/DMD.H
new file mode 100644
index 0000000000..b468e85f86
--- /dev/null
+++ b/src/functionObjects/field/DMD/DMD.H
@@ -0,0 +1,278 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020-2021 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::DMD
+
+Group
+    grpFieldFunctionObjects
+
+Description
+    Computes a dynamic mode decomposition model on a specified field.
+
+    Dynamic mode decomposition (i.e. DMD) is a data-driven
+    dimensionality reduction method. DMD is being used as a mathematical
+    processing tool to reveal dominant modes out of a given field (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,
+    tractable, and even reproducible without computing simulations.
+    DMD only relies on input data, therefore it is an equation-free approach.
+
+    References:
+    \verbatim
+        DMD characteristics:
+            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(s) <!--
+                   -->                    | \<time\>/\<inputField\>(s)
+      output file  | dat           | postProcessing/\<FO\>/\<time\>/\<file\>(s)
+      output field | volVectorField(s)    | \<time\>/\<outputField\>(s)
+    \endtable
+
+    where \c \<Type\>=Scalar/Vector/SphericalTensor/SymmTensor/Tensor.
+
+    Output fields:
+    \verbatim
+      modeRe_<modeIndex>_<field>_<FO> | Real part of a mode field
+      modeIm_<modeIndex>_<field>_<FO> | Imaginary part of a mode field
+    \endverbatim
+
+    Output files:
+    \verbatim
+      dynamics_<field>.dat            | Dynamics data for each mode
+      filteredDynamics_<field>.dat    | Filtered dynamics data for each mode
+    \endverbatim
+
+    wherein for each mode, the following quantities are output into files:
+    \vartable
+      freq        | Frequency
+      mag         | Magnitude
+      ampRe       | Amplitude coefficients (real part)
+      ampIm       | Amplitude coefficients (imaginary part)
+      evalRe      | Eigenvalue (real part)
+      evalIm      | Eigenvalue (imaginary part)
+    \endvartable
+
+Usage
+    Minimal example by using \c system/controlDict.functions:
+    \verbatim
+    DMD1
+    {
+        // Mandatory entries (unmodifiable)
+        type                DMD;
+        libs                (fieldFunctionObjects);
+        DMDModel            <DMDModel>;
+        field               <inputField>;
+
+        // Optional entries (unmodifiable)
+        patch               <patchName>;
+
+        // Mandatory/Optional (inherited) entries
+        ...
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property     | Description                        | Type | Reqd | Deflt
+      type         | Type name: DMD                     | word | yes  | -
+      libs         | Library name: fieldFunctionObjects | word | yes  | -
+      DMDModel     | Name of specified DMD model        | word | yes  | -
+      field        | Name of operand field              | word | yes  | -
+      patch        | Name of operand patch              | word | no   | null
+    \endtable
+
+    Options for the \c DMDModel entry:
+    \verbatim
+      STDMD    | Streaming total dynamic mode decomposition
+    \endverbatim
+
+    The inherited entries are elaborated in:
+     - \link functionObject.H \endlink
+     - \link writeFile.H \endlink
+
+    Minimal example by using the \c postProcess utility:
+    \verbatim
+        <solver> -postProcess -fields '(U p)' -time '10:'
+    \endverbatim
+
+Note
+  - Warning: DMD is an active research area at the time of writing;
+    therefore, there could be cases whereat oddities can be seen.
+
+See also
+  - Foam::functionObjects::fvMeshFunctionObject
+  - Foam::functionObjects::writeFile
+  - Foam::DMDModel
+  - Foam::DMDModels::STDMD
+
+SourceFiles
+    DMD.C
+    DMDTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef functionObjects_DMD_H
+#define functionObjects_DMD_H
+
+#include "fvMeshFunctionObject.H"
+#include "RectangularMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward Declarations
+class DMDModel;
+
+namespace functionObjects
+{
+
+/*---------------------------------------------------------------------------*\
+                              Class DMD Declaration
+\*---------------------------------------------------------------------------*/
+
+class DMD
+:
+    public fvMeshFunctionObject
+{
+    typedef RectangularMatrix<scalar> RMatrix;
+
+    // Private Data
+
+        //- Dynamic mode decomposition model
+        autoPtr<DMDModel> DMDModelPtr_;
+
+        //- Snapshot matrix (effectively a column vector)
+        //  Upper half = current-time snapshot slot
+        //  Lower half = previous-time snapshot slot
+        //  A snapshot is an input dataset to be processed per execution step
+        RMatrix z_;
+
+        //- Name of operand field
+        const word fieldName_;
+
+        //- Name of operand patch
+        const word patch_;
+
+        //- Number of elements in a snapshot
+        label nSnap_;
+
+        //- Current execution-step index of DMD,
+        //- not necessarily that of the simulation
+        label step_;
+
+
+    // Private Member Functions
+
+        // Process
+
+            //- Initialise snapshot at the first-execution step
+            //  Initialisation at the ctor or read level is not possible
+            //  since the operand field is not available in the database
+            void initialise();
+
+            //- Create operand snapshot by using
+            //- current-time and previous-time operand fields
+            void snapshot();
+
+            //- Get operand field based on its base type
+            template<class Type>
+            bool getSnapshot();
+
+            //- Get operand field based on its geometric field type
+            //  Move previous-time field into previous-time slot in snapshot
+            //  copy new current-time field into current-time slot in snapshot
+            template<class GeoFieldType>
+            bool getSnapshotField();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("DMD");
+
+
+    // Constructors
+
+        //- Construct from Time and dictionary
+        DMD
+        (
+            const word& name,
+            const Time& runTime,
+            const dictionary& dict
+        );
+
+        //- No copy construct
+        DMD(const DMD&) = delete;
+
+        //- No copy assignment
+        void operator=(const DMD&) = delete;
+
+
+    //- Destructor
+    virtual ~DMD() = default;
+
+
+    // Member Functions
+
+        //- Read DMD settings
+        virtual bool read(const dictionary& dict);
+
+        //- Execute DMD
+        virtual bool execute();
+
+        //- Write DMD results
+        virtual bool write();
+
+        //- Write DMD results
+        virtual bool end();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace functionObjects
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "DMDTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/DMD/DMDModels/DMDModel/DMDModel.C b/src/functionObjects/field/DMD/DMDModels/DMDModel/DMDModel.C
new file mode 100644
index 0000000000..c43727548c
--- /dev/null
+++ b/src/functionObjects/field/DMD/DMDModels/DMDModel/DMDModel.C
@@ -0,0 +1,78 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020-2021 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 "DMDModel.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(DMDModel, 0);
+    defineRunTimeSelectionTable(DMDModel, dictionary);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::DMDModel::DMDModel
+(
+    const fvMesh& mesh,
+    const word& name,
+    const dictionary& dict
+)
+:
+    writeFile(mesh, name, typeName, dict, false),
+    mesh_(mesh),
+    name_(name)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::label Foam::DMDModel::nComponents(const word& fieldName) const
+{
+    label nComps = 0;
+    bool processed = false;
+    processed = processed || nComponents<scalar>(fieldName, nComps);
+    processed = processed || nComponents<vector>(fieldName, nComps);
+    processed = processed || nComponents<sphericalTensor>(fieldName, nComps);
+    processed = processed || nComponents<symmTensor>(fieldName, nComps);
+    processed = processed || nComponents<tensor>(fieldName, nComps);
+
+    if (!processed)
+    {
+        FatalErrorInFunction
+            << "  # Unknown type of input field during initialisation = "
+            << fieldName << " #" << nl
+            << exit(FatalError);
+    }
+
+    return nComps;
+}
+
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/DMD/DMDModels/DMDModel/DMDModel.H b/src/functionObjects/field/DMD/DMDModels/DMDModel/DMDModel.H
new file mode 100644
index 0000000000..8081db5a01
--- /dev/null
+++ b/src/functionObjects/field/DMD/DMDModels/DMDModel/DMDModel.H
@@ -0,0 +1,203 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020-2021 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/>.
+
+Namespace
+    Foam::DMDModels
+
+Description
+    A namespace for various dynamic mode
+    decomposition (DMD) model implementations.
+
+Class
+    Foam::DMDModel
+
+Description
+    Abstract base class for DMD models to handle DMD
+    characteristics for the \c DMD function object.
+
+SourceFiles
+    DMDModel.C
+    DMDModelNew.C
+    DMDModelTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef DMDModel_H
+#define DMDModel_H
+
+#include "fvMesh.H"
+#include "dictionary.H"
+#include "HashSet.H"
+#include "runTimeSelectionTables.H"
+#include "writeFile.H"
+#include "OFstream.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                            Class DMDModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class DMDModel
+:
+    public functionObjects::writeFile
+{
+    typedef SquareMatrix<scalar> SMatrix;
+    typedef RectangularMatrix<scalar> RMatrix;
+
+
+protected:
+
+    // Protected Data
+
+        //- Reference to the mesh
+        const fvMesh& mesh_;
+
+        //- Name of operand function object
+        const word name_;
+
+
+    // Protected Member Functions
+
+        // Evaluation
+
+            //- Compute and write mode dynamics
+            virtual bool dynamics() = 0;
+
+            //- Compute and write modes
+            virtual bool modes() = 0;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("DMDModel");
+
+
+    // Declare runtime constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            DMDModel,
+            dictionary,
+            (
+                const fvMesh& mesh,
+                const word& name,
+                const dictionary& dict
+            ),
+            (mesh, name, dict)
+        );
+
+
+    // Selectors
+
+        //- Return a reference to the selected DMD model
+        static autoPtr<DMDModel> New
+        (
+            const fvMesh& mesh,
+            const word& name,
+            const dictionary& dict
+        );
+
+
+    // Constructors
+
+        //- Construct from components
+        DMDModel
+        (
+            const fvMesh& mesh,
+            const word& name,
+            const dictionary& dict
+        );
+
+        //- No copy construct
+        DMDModel(const DMDModel&) = delete;
+
+        //- No copy assignment
+        void operator=(const DMDModel&) = delete;
+
+
+    //- Destructor
+    virtual ~DMDModel() = default;
+
+
+    // Member Functions
+
+        // Process
+
+            //- Initialise model data members with a given snapshot
+            virtual bool initialise(const RMatrix& snapshot) = 0;
+
+            //- Update model data members with a given snapshot
+            virtual bool update(const RMatrix& snapshot) = 0;
+
+            //- Compute and write modes and
+            //- mode dynamics of model data members
+            virtual bool fit() = 0;
+
+            //- Compute and write a reconstruction of flow field
+            //- based on given modes and mode dynamics (currently no-op)
+            virtual void reconstruct(const wordList modes)
+            {
+                NotImplemented;
+            }
+
+
+        // Access
+
+            //- Return number of components of the base type of a given field
+            label nComponents(const word& fieldName) const;
+
+            //- Get the number of components of the base type of a given field
+            template<class Type>
+            bool nComponents(const word& fieldName, label& nComps) const;
+
+
+        // IO
+
+            //- Read model settings
+            virtual bool read(const dictionary& dict) = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "DMDModelTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/DMD/DMDModels/DMDModel/DMDModelNew.C b/src/functionObjects/field/DMD/DMDModels/DMDModel/DMDModelNew.C
new file mode 100644
index 0000000000..e3a72a4a5c
--- /dev/null
+++ b/src/functionObjects/field/DMD/DMDModels/DMDModel/DMDModelNew.C
@@ -0,0 +1,58 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020-2021 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 "DMDModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::DMDModel> Foam::DMDModel::New
+(
+    const fvMesh& mesh,
+    const word& name,
+    const dictionary& dict
+)
+{
+    const word modelType(dict.get<word>("DMDModel"));
+
+    auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);
+
+    if (!cstrIter.found())
+    {
+        FatalIOErrorInLookup
+        (
+            dict,
+            "DMDModel",
+            modelType,
+            *dictionaryConstructorTablePtr_
+        ) << exit(FatalIOError);
+    }
+
+    return autoPtr<DMDModel>(cstrIter()(mesh, name, dict));
+}
+
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/DMD/DMDModels/DMDModel/DMDModelTemplates.C b/src/functionObjects/field/DMD/DMDModels/DMDModel/DMDModelTemplates.C
new file mode 100644
index 0000000000..9281dd2385
--- /dev/null
+++ b/src/functionObjects/field/DMD/DMDModels/DMDModel/DMDModelTemplates.C
@@ -0,0 +1,54 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 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 Type>
+bool Foam::DMDModel::nComponents(const word& fieldName, label& nComps) const
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
+    typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
+
+    if (mesh_.foundObject<VolFieldType>(fieldName))
+    {
+        nComps = pTraits<typename VolFieldType::value_type>::nComponents;
+        return true;
+    }
+    else if (mesh_.foundObject<SurfaceFieldType>(fieldName))
+    {
+        nComps = pTraits<typename SurfaceFieldType::value_type>::nComponents;
+        return true;
+    }
+
+    return false;
+}
+
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/DMD/DMDModels/derived/STDMD/STDMD.C b/src/functionObjects/field/DMD/DMDModels/derived/STDMD/STDMD.C
new file mode 100644
index 0000000000..ac180feb7a
--- /dev/null
+++ b/src/functionObjects/field/DMD/DMDModels/derived/STDMD/STDMD.C
@@ -0,0 +1,988 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020-2021 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 "EigenMatrix.H"
+#include "QRMatrix.H"
+#include "addToRunTimeSelectionTable.H"
+
+using namespace Foam::constant::mathematical;
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace DMDModels
+{
+    defineTypeNameAndDebug(STDMD, 0);
+    addToRunTimeSelectionTable(DMDModel, STDMD, dictionary);
+}
+}
+
+const Foam::Enum
+<
+    Foam::DMDModels::STDMD::modeSorterType
+>
+Foam::DMDModels::STDMD::modeSorterTypeNames
+({
+    { modeSorterType::FIRST_SNAPSHOT, "firstSnapshot" },
+    { modeSorterType::KIEWAT , "kiewat" },
+    { modeSorterType::KOU_ZHANG , "kouZhang" }
+});
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::scalar Foam::DMDModels::STDMD::L2norm(const RMatrix& z) const
+{
+    #ifdef FULLDEBUG
+    // Check if the given RectangularMatrix is effectively a column vector
+    if (z.n() != 1)
+    {
+        FatalErrorInFunction
+            << "  # Input matrix is not a column vector. #"
+            << exit(FatalError);
+    }
+    #endif
+    const bool noSqrt = true;
+    scalar result = z.columnNorm(0, noSqrt);
+    reduce(result, sumOp<scalar>());
+
+    // Heuristic addition to avoid very small or zero norm
+    return max(SMALL, Foam::sqrt(result));
+}
+
+
+Foam::RectangularMatrix<Foam::scalar> Foam::DMDModels::STDMD::orthonormalise
+(
+    RMatrix ez
+) const
+{
+    RMatrix dz(Q_.n(), 1, Zero);
+
+    for (label i = 0; i < nGramSchmidt_; ++i)
+    {
+        dz = Q_ & ez;
+        reduce(dz, sumOp<RMatrix>());
+        ez -= Q_*dz;
+    }
+
+    return ez;
+}
+
+
+void Foam::DMDModels::STDMD::expand(const RMatrix& ez, const scalar ezNorm)
+{
+    Info<< tab << "Expanding orthonormal basis for field: " << fieldName_
+        << endl;
+
+    // Stack a column "(ez/ezNorm)" to "Q"
+    Q_.resize(Q_.m(), Q_.n() + 1);
+    Q_.subColumn(Q_.n() - 1) = ez/ezNorm;
+
+    // Stack a row (Zero) and column (Zero) to "G"
+    G_.resize(G_.m() + 1);
+}
+
+
+void Foam::DMDModels::STDMD::compress()
+{
+    Info<< tab << "Compressing orthonormal basis for field: " << fieldName_
+        << endl;
+
+    RMatrix q(1, 1, Zero);
+
+    if (Pstream::master())
+    {
+        const bool symmetric = true;
+        const EigenMatrix<scalar> EM(G_, symmetric);
+        const SquareMatrix<scalar>& EVecs = EM.EVecs();
+        DiagonalMatrix<scalar> EVals(EM.EValsRe());
+
+        // Sort eigenvalues in descending order, and track indices
+        const auto descend = [&](scalar a, scalar b){ return a > b; };
+        const List<label> permutation(EVals.sortPermutation(descend));
+        EVals.applyPermutation(permutation);
+        EVals.resize(EVals.size() - 1);
+
+        // Update "G"
+        G_ = SMatrix(maxRank_, Zero);
+        G_.diag(EVals);
+
+        q.resize(Q_.n(), maxRank_);
+        for (label i = 0; i < maxRank_; ++i)
+        {
+            q.subColumn(i) = EVecs.subColumn(permutation[i]);
+        }
+    }
+    Pstream::scatter(G_);
+    Pstream::scatter(q);
+
+    // Update "Q"
+    Q_ = Q_*q;
+}
+
+
+Foam::SquareMatrix<Foam::scalar> Foam::DMDModels::STDMD::
+reducedKoopmanOperator()
+{
+    Info<< tab << "Computing Atilde" << endl;
+
+    Info<< tab << "Computing local Rx" << endl;
+
+    RMatrix Rx(1, 1, Zero);
+    {
+        QRMatrix<RMatrix> QRM
+        (
+            Qupper_,
+            QRMatrix<RMatrix>::outputTypes::REDUCED_R,
+            QRMatrix<RMatrix>::colPivoting::FALSE
+        );
+        Rx = QRM.R();
+    }
+    Rx.round();
+
+    // Convenience objects A1, A2, A3 to compute "Atilde"
+    RMatrix A1(1, 1, Zero);
+
+    if (Pstream::parRun())
+    {
+        // Parallel direct tall-skinny QR decomposition
+        // (BGD:Fig. 5) & (DGHL:Fig. 2)
+        // Tests revealed that the distribution of "Q" does not affect
+        // the final outcome of TSQR decomposition up to sign
+
+        PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
+
+        const label myProcNo = Pstream::myProcNo();
+        const label procNoInSubset = myProcNo % nAgglomerationProcs_;
+
+        // Send Rx from sender subset neighbours to receiver subset master
+        if (procNoInSubset != 0)
+        {
+            const label procNoSubsetMaster = myProcNo - procNoInSubset;
+
+            UOPstream toSubsetMaster(procNoSubsetMaster, pBufs);
+            toSubsetMaster << Rx;
+
+            Rx.clear();
+        }
+
+        pBufs.finishedSends();
+
+        // Receive Rx by receiver subset masters
+        if (procNoInSubset == 0)
+        {
+            // Accept only subset masters possessing sender subset neighbours
+            if (myProcNo + 1 < Pstream::nProcs())
+            {
+                for
+                (
+                    label nbr = myProcNo + 1;
+                    (
+                        nbr < myProcNo + nAgglomerationProcs_
+                     && nbr < Pstream::nProcs()
+                    );
+                    ++nbr
+                )
+                {
+                    RMatrix recvMtrx;
+
+                    UIPstream fromNbr(nbr, pBufs);
+                    fromNbr >> recvMtrx;
+
+                    // Append received Rx to Rx of receiver subset masters
+                    if (recvMtrx.size() > 0)
+                    {
+                        Rx.resize(Rx.m() + recvMtrx.m(), Rx.n());
+                        Rx.subMatrix
+                        (
+                            Rx.m() - recvMtrx.m(),
+                            Rx.n() - recvMtrx.n()
+                        ) = recvMtrx;
+                    }
+                }
+            }
+
+            // Apply interim QR decomposition on Rx of receiver subset masters
+            QRMatrix<RMatrix> QRM
+            (
+                Rx,
+                QRMatrix<RMatrix>::outputTypes::REDUCED_R,
+                QRMatrix<RMatrix>::storeMethods::IN_PLACE,
+                QRMatrix<RMatrix>::colPivoting::FALSE
+            );
+            Rx.round();
+        }
+
+        pBufs.clear();
+
+        // Send interim Rx from subset masters to the master
+        if (procNoInSubset == 0)
+        {
+            if (!Pstream::master())
+            {
+                UOPstream toMaster(Pstream::masterNo(), pBufs);
+                toMaster << Rx;
+            }
+        }
+
+        pBufs.finishedSends();
+
+        // Receive interim Rx by the master, and apply final operations
+        if (Pstream::master())
+        {
+            for
+            (
+                label nbr = Pstream::masterNo() + nAgglomerationProcs_;
+                nbr < Pstream::nProcs();
+                nbr += nAgglomerationProcs_
+            )
+            {
+                RMatrix recvMtrx;
+
+                UIPstream fromSubsetMaster(nbr, pBufs);
+                fromSubsetMaster >> recvMtrx;
+
+                // Append received Rx to Rx of the master
+                if (recvMtrx.size() > 0)
+                {
+                    Rx.resize(Rx.m() + recvMtrx.m(), Rx.n());
+                    Rx.subMatrix
+                    (
+                        Rx.m() - recvMtrx.m(),
+                        Rx.n() - recvMtrx.n()
+                    ) = recvMtrx;
+                }
+            }
+
+            Info<< tab << "Computing TSQR" << endl;
+            QRMatrix<RMatrix> QRM
+            (
+                Rx,
+                QRMatrix<RMatrix>::outputTypes::REDUCED_R,
+                QRMatrix<RMatrix>::storeMethods::IN_PLACE,
+                QRMatrix<RMatrix>::colPivoting::FALSE
+            );
+            Rx.round();
+
+            // Rx produced by TSQR is unique up to the sign, hence the revert
+            for (scalar& x : Rx)
+            {
+                x *= -1;
+            }
+
+            Info<< tab << "Computing RxInv" << endl;
+            RxInv_ = pinv(Rx);
+
+            Info<< tab << "Computing A1" << endl;
+            A1 = RxInv_*RMatrix(pinv(Rx*(G_^Rx)));
+            Rx.clear();
+        }
+        Pstream::scatter(RxInv_);
+        Pstream::scatter(A1);
+
+        Info<< tab << "Computing A2" << endl;
+        SMatrix A2(Qupper_ & Qlower_);
+        reduce(A2, sumOp<SMatrix>());
+        Qlower_.clear();
+
+        Info<< tab << "Computing A3" << endl;
+        SMatrix A3(Qupper_ & Qupper_);
+        reduce(A3, sumOp<SMatrix>());
+
+        Info<< tab << "Computing Atilde" << endl;
+        // by optimized matrix chain multiplication
+        // obtained by dynamic programming memoisation
+        return SMatrix(RxInv_ & ((A2*(G_*A3))*A1));
+    }
+    else
+    {
+        Info<< tab << "Computing RxInv" << endl;
+        RxInv_ = pinv(Rx);
+
+        Info<< tab << "Computing A1" << endl;
+        A1 = RxInv_*RMatrix(pinv(Rx*(G_^Rx)));
+        Rx.clear();
+
+        Info<< tab << "Computing A2" << endl;
+        SMatrix A2(Qupper_ & Qlower_);
+        Qlower_.clear();
+
+        Info<< tab << "Computing A3" << endl;
+        SMatrix A3(Qupper_ & Qupper_);
+
+        Info<< tab << "Computing Atilde" << endl;
+        return SMatrix(RxInv_ & ((A2*(G_*A3))*A1));
+    }
+}
+
+
+bool Foam::DMDModels::STDMD::eigendecomposition(SMatrix& Atilde)
+{
+    bool fail = false;
+
+    // Compute eigenvalues, and clip eigenvalues with values < "minEval"
+    if (Pstream::master())
+    {
+        Info<< tab << "Computing eigendecomposition" << endl;
+
+        // Replace "Atilde" by eigenvalues (in-place eigendecomposition)
+        const EigenMatrix<scalar> EM(Atilde);
+        const DiagonalMatrix<scalar>& evalsRe = EM.EValsRe();
+        const DiagonalMatrix<scalar>& evalsIm = EM.EValsIm();
+
+        evals_.resize(evalsRe.size());
+        evecs_ = RCMatrix(EM.complexEVecs());
+
+        // Convert scalar eigenvalue pairs to complex eigenvalues
+        label i = 0;
+        for (auto& eval : evals_)
+        {
+            eval = complex(evalsRe[i], evalsIm[i]);
+            ++i;
+        }
+
+        Info<< tab << "Filtering eigenvalues" << endl;
+
+        List<complex> cp(evals_.size());
+        auto it =
+            std::copy_if
+            (
+                evals_.cbegin(),
+                evals_.cend(),
+                cp.begin(),
+                [&](const complex& x){ return mag(x) > minEval_; }
+            );
+        cp.resize(std::distance(cp.begin(), it));
+
+        if (cp.size() == 0)
+        {
+            WarningInFunction
+                << "No eigenvalue with mag(eigenvalue) larger than "
+                << "minEval = " << minEval_ << " was found." << nl
+                << "    Input field may contain only zero-value elements," << nl
+                << "    e.g. no-slip velocity condition on a given patch." << nl
+                << "    Otherwise, please decrease the value of minEval." << nl
+                << "    Skipping further dynamics/mode computations." << nl
+                << endl;
+
+            fail = true;
+        }
+        else
+        {
+            evals_ = cp;
+        }
+    }
+    Pstream::scatter(fail);
+
+    if (fail)
+    {
+        return false;
+    }
+
+    Pstream::scatter(evals_);
+    Pstream::scatter(evecs_);
+
+    Atilde.clear();
+
+    return true;
+}
+
+
+void Foam::DMDModels::STDMD::frequencies()
+{
+    if (Pstream::master())
+    {
+        Info<< tab << "Computing frequencies" << endl;
+
+        freqs_.resize(evals_.size());
+
+        // Frequency equation (K:Eq. 81)
+        auto frequencyEquation =
+            [&](const complex& eval)
+            {
+                return Foam::log(max(eval, complex(SMALL))).imag()/(twoPi*dt_);
+            };
+
+        // Compute frequencies
+        std::transform
+        (
+            evals_.cbegin(),
+            evals_.cend(),
+            freqs_.begin(),
+            frequencyEquation
+        );
+
+        Info<< tab << "Computing frequency indices" << endl;
+
+        // Initialise iterator by the first search
+        auto margin = [&](const scalar& x){ return (fMin_ <= x && x < fMax_); };
+
+        auto it = std::find_if(freqs_.cbegin(), freqs_.cend(), margin);
+
+        while (it != freqs_.end())
+        {
+            freqsi_.append(std::distance(freqs_.cbegin(), it));
+            it = std::find_if(std::next(it), freqs_.cend(), margin);
+        }
+    }
+    Pstream::scatter(freqs_);
+    Pstream::scatter(freqsi_);
+}
+
+
+void Foam::DMDModels::STDMD::amplitudes()
+{
+    const IOField<scalar> snapshot0
+    (
+        IOobject
+        (
+            "snapshot0_" + name_ + "_" + fieldName_,
+            timeName0_,
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            false
+        )
+    );
+
+    RMatrix snapshot(1, 1, Zero);
+    if (!empty_)
+    {
+        snapshot.resize(Qupper_.m(), 1);
+        std::copy(snapshot0.cbegin(), snapshot0.cend(), snapshot.begin());
+    }
+
+    RMatrix R((RxInv_.T()^Qupper_)*snapshot);
+    reduce(R, sumOp<RMatrix>());
+
+    if (Pstream::master())
+    {
+        Info<< tab << "Computing amplitudes" << endl;
+
+        amps_.resize(R.m());
+        const RCMatrix pEvecs(pinv(evecs_));
+
+        // amps_ = pEvecs*R;
+        for (label i = 0; i < amps_.size(); ++i)
+        {
+            for (label j = 0; j < R.m(); ++j)
+            {
+                amps_[i] += pEvecs(i, j)*R(j, 0);
+            }
+        }
+    }
+    Pstream::scatter(amps_);
+}
+
+
+void Foam::DMDModels::STDMD::magnitudes()
+{
+    if (Pstream::master())
+    {
+        Info<< tab << "Computing magnitudes" << endl;
+
+        mags_.resize(amps_.size());
+
+        Info<< tab << "Sorting modes with ";
+
+        switch (modeSorter_)
+        {
+            case modeSorterType::FIRST_SNAPSHOT:
+            {
+                Info<< "method of first snapshot" << endl;
+
+                std::transform
+                (
+                    amps_.cbegin(),
+                    amps_.cend(),
+                    mags_.begin(),
+                    [&](const complex& val){ return mag(val); }
+                );
+                break;
+            }
+
+            case modeSorterType::KIEWAT:
+            {
+                Info<< "modified weighted amplitude scaling method" << endl;
+
+                //  Eigendecomposition returns evecs with
+                //  the unity norm, hence "modeNorm = 1"
+                const scalar modeNorm = 1;
+                const scalar pr = 1;
+                List<scalar> w(step_);
+                std::iota(w.begin(), w.end(), 1);
+                w = sin(twoPi/step_*(w - 1 - 0.25*step_))*pr + pr;
+
+                forAll(amps_, i)
+                {
+                    mags_[i] = sorter(w, amps_[i], evals_[i], modeNorm);
+                }
+
+                break;
+            }
+
+            case modeSorterType::KOU_ZHANG:
+            {
+                Info<< "weighted amplitude scaling method" << endl;
+
+                const scalar modeNorm = 1;
+                const List<scalar> w(step_, 1.0);
+
+                forAll(amps_, i)
+                {
+                    mags_[i] = sorter(w, amps_[i], evals_[i], modeNorm);
+                }
+
+                break;
+            }
+
+            default:
+                break;
+        }
+
+        Info<< tab << "Computing magnitude indices" << endl;
+
+        magsi_ = freqsi_;
+
+        auto descend =
+            [&](const label i1, const label i2)
+            {
+                return !(mags_[i1] < mags_[i2]);
+            };
+
+        std::sort(magsi_.begin(), magsi_.end(), descend);
+    }
+    Pstream::scatter(mags_);
+    Pstream::scatter(magsi_);
+}
+
+
+Foam::scalar Foam::DMDModels::STDMD::sorter
+(
+    const List<scalar>& weight,
+    const complex& amplitude,
+    const complex& eigenvalue,
+    const scalar modeNorm
+) const
+{
+    // Omit eigenvalues with very large or very small mags
+    if (!(mag(eigenvalue) < GREAT && mag(eigenvalue) > VSMALL))
+    {
+        Info<< "  Returning zero magnitude for mag(eigenvalue) = "
+            << mag(eigenvalue) << endl;
+
+        return 0;
+    }
+
+    // Omit eigenvalue-STDMD step combinations that pose a risk of overflow
+    if (mag(eigenvalue)*step_ > sortLimiter_)
+    {
+        Info<< "  Returning zero magnitude for"
+            << " mag(eigenvalue) = " << mag(eigenvalue)
+            << " current index = " << step_
+            << " sortLimiter = " << sortLimiter_
+            << endl;
+
+        return 0;
+    }
+
+    scalar magnitude = 0;
+
+    for (label j = 0; j < step_; ++j)
+    {
+        magnitude += weight[j]*modeNorm*mag(amplitude*pow(eigenvalue, j + 1));
+    }
+
+    return magnitude;
+}
+
+
+bool Foam::DMDModels::STDMD::dynamics()
+{
+    SMatrix Atilde(reducedKoopmanOperator());
+    G_.clear();
+
+    if(!eigendecomposition(Atilde))
+    {
+        return false;
+    }
+
+    frequencies();
+
+    amplitudes();
+
+    magnitudes();
+
+    return true;
+}
+
+
+bool Foam::DMDModels::STDMD::modes()
+{
+    Info<< tab << "Computing modes" << endl;
+
+    bool processed = false;
+    processed = processed || modes<scalar>();
+    processed = processed || modes<vector>();
+    processed = processed || modes<sphericalTensor>();
+    processed = processed || modes<symmTensor>();
+    processed = processed || modes<tensor>();
+
+    if (!processed)
+    {
+        return false;
+    }
+
+    return true;
+}
+
+
+void Foam::DMDModels::STDMD::writeToFile(const word& fileName) const
+{
+    Info<< tab << "Writing objects of dynamics" << endl;
+
+    // Write objects of dynamics
+    {
+        autoPtr<OFstream> osPtr =
+            createFile
+            (
+                fileName + "_" + fieldName_,
+                mesh_.time().timeOutputValue()
+            );
+        OFstream& os = osPtr.ref();
+
+        writeFileHeader(os);
+
+        for (const auto& i : labelRange(0, freqs_.size()))
+        {
+            os  << freqs_[i] << tab
+                << mags_[i] << tab
+                << amps_[i].real() << tab
+                << amps_[i].imag() << tab
+                << evals_[i].real() << tab
+                << evals_[i].imag() << endl;
+        }
+    }
+
+    Info<< tab << "Ends output processing for field: " << fieldName_
+        << endl;
+}
+
+
+void Foam::DMDModels::STDMD::writeFileHeader(Ostream& os) const
+{
+    writeHeader(os, "DMD 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::DMDModels::STDMD::filter()
+{
+    Info<< tab << "Filtering objects of dynamics" << endl;
+
+    // Filter objects according to iMags
+    filterIndexed(evals_, magsi_);
+    filterIndexed(evecs_, magsi_);
+    filterIndexed(freqs_, magsi_);
+    filterIndexed(amps_, magsi_);
+    filterIndexed(mags_, magsi_);
+
+    // Clip objects if need be (assuming objects have the same len)
+    if (freqs_.size() > nModes_)
+    {
+        evals_.resize(nModes_);
+        evecs_.resize(evecs_.m(), nModes_);
+        freqs_.resize(nModes_);
+        amps_.resize(nModes_);
+        mags_.resize(nModes_);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::DMDModels::STDMD::STDMD
+(
+    const fvMesh& mesh,
+    const word& name,
+    const dictionary& dict
+)
+:
+    DMDModel(mesh, name, dict),
+    modeSorter_
+    (
+        modeSorterTypeNames.getOrDefault
+        (
+            "modeSorter",
+            dict,
+            modeSorterType::FIRST_SNAPSHOT
+        )
+    ),
+    Q_(),
+    G_(),
+    Qupper_(1, 1, Zero),
+    Qlower_(1, 1, Zero),
+    RxInv_(1, 1, Zero),
+    evals_(Zero),
+    evecs_(1, 1, Zero),
+    freqs_(Zero),
+    freqsi_(),
+    amps_(Zero),
+    mags_(Zero),
+    magsi_(Zero),
+    fieldName_(dict.get<word>("field")),
+    patch_(dict.getOrDefault<word>("patch", word::null)),
+    timeName0_(),
+    minBasis_(0),
+    minEval_(0),
+    dt_(0),
+    sortLimiter_(500),
+    fMin_(0),
+    fMax_(pTraits<label>::max),
+    nGramSchmidt_(5),
+    maxRank_(pTraits<label>::max),
+    step_(0),
+    nModes_(pTraits<label>::max),
+    nAgglomerationProcs_(20),
+    empty_(false)
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+bool Foam::DMDModels::STDMD::read(const dictionary& dict)
+{
+    writeToFile_ = dict.getOrDefault<bool>("writeToFile", true);
+
+    modeSorter_ =
+        modeSorterTypeNames.getOrDefault
+        (
+            "modeSorter",
+            dict,
+            modeSorterType::FIRST_SNAPSHOT
+        );
+
+    minBasis_ =
+        dict.getCheckOrDefault<scalar>("minBasis", 1e-8, scalarMinMax::ge(0));
+
+    minEval_ =
+        dict.getCheckOrDefault<scalar>("minEVal", 1e-8, scalarMinMax::ge(0));
+
+    dt_ =
+        dict.getCheckOrDefault<scalar>
+        (
+            "interval",
+            (
+                dict.getCheck<label>("executeInterval", labelMinMax::ge(1))
+               *mesh_.time().deltaT().value()
+            ),
+            scalarMinMax::ge(0)
+        );
+
+    sortLimiter_ =
+        dict.getCheckOrDefault<scalar>("sortLimiter", 500, scalarMinMax::ge(1));
+
+    nGramSchmidt_ =
+        dict.getCheckOrDefault<label>("nGramSchmidt", 5, labelMinMax::ge(1));
+
+    maxRank_ =
+        dict.getCheckOrDefault<label>
+        (
+            "maxRank",
+            pTraits<label>::max,
+            labelMinMax::ge(1)
+        );
+
+    nModes_ =
+        dict.getCheckOrDefault<label>
+        (
+            "nModes",
+            pTraits<label>::max,
+            labelMinMax::ge(1)
+        );
+
+    fMin_ = dict.getCheckOrDefault<label>("fMin", 0, labelMinMax::ge(0));
+
+    fMax_ =
+        dict.getCheckOrDefault<label>
+        (
+            "fMax",
+            pTraits<label>::max,
+            labelMinMax::ge(fMin_ + 1)
+        );
+
+    nAgglomerationProcs_ =
+        dict.getCheckOrDefault<label>
+        (
+            "nAgglomerationProcs",
+            20,
+            labelMinMax::ge(1)
+        );
+
+    Info<< tab << "Settings are read for:" << nl
+        << "    field: " << fieldName_ << nl
+        << "    modeSorter: " << modeSorterTypeNames[modeSorter_] << nl
+        << "    nModes: " << nModes_ << nl
+        << "    maxRank: " << maxRank_ << nl
+        << "    nGramSchmidt: " << nGramSchmidt_ << nl
+        << "    fMin: " << fMin_ << nl
+        << "    fMax: " << fMax_ << nl
+        << "    minBasis: " << minBasis_ << nl
+        << "    minEVal: " << minEval_ << nl
+        << "    sortLimiter: " << sortLimiter_ << nl
+        << "    nAgglomerationProcs: " << nAgglomerationProcs_ << nl
+        << endl;
+
+    return true;
+}
+
+
+bool Foam::DMDModels::STDMD::initialise(const RMatrix& z)
+{
+    const scalar norm = L2norm(z);
+
+    if (mag(norm) > 0)
+    {
+        // First-processed snapshot required by the mode-sorting
+        // algorithms at the final output computations (K:p. 43)
+        {
+            const label nSnap = z.m()/2;
+
+            if (nSnap == 0)
+            {
+                empty_ = true;
+            }
+
+            scalarField snapshot0(nSnap);
+            std::copy(z.cbegin(), z.cbegin() + nSnap, snapshot0.begin());
+            timeName0_ = mesh_.time().timeName();
+
+            IOField<scalar>
+            (
+                IOobject
+                (
+                    "snapshot0_" + name_ + "_" + fieldName_,
+                    timeName0_,
+                    mesh_,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                snapshot0
+            ).write();
+        }
+
+        Q_ = z/norm;
+        G_ = SMatrix(1);
+        G_(0,0) = sqr(norm);
+
+        ++step_;
+
+        return true;
+    }
+
+    return false;
+}
+
+
+bool Foam::DMDModels::STDMD::update(const RMatrix& z)
+{
+    {
+        //- Working copy of the augmented snapshot matrix "z"
+        //- being used in the classical Gram-Schmidt process
+        const RMatrix ez(orthonormalise(z));
+
+        // Check basis for "z" and, if necessary, expand "Q" and "G"
+        const scalar ezNorm = L2norm(ez);
+        const scalar zNorm = L2norm(z);
+
+        if (ezNorm/zNorm > minBasis_)
+        {
+            expand(ez, ezNorm);
+        }
+    }
+
+    // Update "G" before the potential orthonormal basis compression
+    {
+        RMatrix zTilde(Q_ & z);
+        reduce(zTilde, sumOp<RMatrix>());
+
+        G_ += SMatrix(zTilde^zTilde);
+    }
+
+    // Compress the orthonormal basis if required
+    if (Q_.n() >= maxRank_)
+    {
+        compress();
+    }
+
+    ++step_;
+
+    return true;
+}
+
+
+bool Foam::DMDModels::STDMD::fit()
+{
+    // DMD statistics and mode evaluation (K:Fig. 16)
+    const label nSnap = Q_.m()/2;
+
+    // Move upper and lower halves of "Q" to new containers
+    Qupper_ = RMatrix(Q_.subMatrix(0, 0, max(nSnap, 1)));
+    Qlower_ = RMatrix(Q_.subMatrix(nSnap, 0, max(nSnap, 1)));
+    Q_.clear();
+
+    if (!dynamics())
+    {
+        return true;
+    }
+
+    modes();
+
+    if (Pstream::master() && writeToFile_)
+    {
+        writeToFile(word("dynamics"));
+
+        filter();
+
+        writeToFile(word("filteredDynamics"));
+    }
+
+    step_ = 0;
+
+    return true;
+}
+
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/DMD/DMDModels/derived/STDMD/STDMD.H b/src/functionObjects/field/DMD/DMDModels/derived/STDMD/STDMD.H
new file mode 100644
index 0000000000..0b86e137fc
--- /dev/null
+++ b/src/functionObjects/field/DMD/DMDModels/derived/STDMD/STDMD.H
@@ -0,0 +1,530 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020-2021 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::DMDModels::STDMD
+
+Description
+    Streaming Total Dynamic Mode Decomposition (i.e. \c STDMD)
+    is a variant of dynamic mode decomposition.
+
+    Among other Dynamic Mode Decomposition (DMD) variants, \c 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
+        DMD 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
+    \endverbatim
+
+Usage
+    Minimal example by using \c system/controlDict.functions:
+    \verbatim
+    DMD1
+    {
+        // Mandatory/Optional entries
+        ...
+
+        // Mandatory entries (unmodifiable)
+        DMDModel            STDMD;
+
+        // Conditional mandatory entries (runtime modifiable)
+
+            // Option-1
+            interval            5.5;
+
+            // Option-2
+            executeInterval     10;
+
+        // Optional entries (runtime modifiable)
+        modeSorter          kiewat;
+        nGramSchmidt        5;
+        maxRank             50;
+        nModes              50;
+        fMin                0;
+        fMax                1000000000;
+        nAgglomerationProcs 20;
+
+        // Optional entries (runtime modifiable, yet not recommended)
+        minBasis            0.00000001;
+        minEVal             0.00000001;
+        sortLimiter         500.0;
+
+        // Mandatory/Optional (inherited) entries
+        ...
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property     | Description                        | Type | Reqd  | Deflt
+      DMDModel     | Type: STDMD                        | word | yes   | -
+      interval     | STDMD time-step size [s]           | scalar | cndtnl <!--
+               --> | executeInterval*(current time-step of the simulation)
+      modeSorter   | Mode-sorting algorithm model       | word | no    <!--
+               --> | firstSnapshot
+      nModes       | Number of output modes in input frequency 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
+      nAgglomerationProcs | Number of processors at each agglomeration <!--
+               --> unit during the computation of reduced Koopman      <!--
+               --> operator                               | label | no | 20
+      minBasis | Orthogonal basis expansion threshold     | scalar| no | 1e-8
+      minEVal  | Min eigenvalue for below eigenvalues are omitted      <!--
+               -->                                        | scalar| no | 1e-8
+      sortLimiter | Max allowable magnitude for <!--
+               -->  mag(eigenvalue)*(number of DMD 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
+      firstSnapshot | First-snapshot amplitude magnitude method
+    \endverbatim
+
+Note
+  - To specify the STDMD time-step size (not necessarily equal to the
+    time step of the simulation), entries of either \c interval or
+    \c executeInterval must be available (highest to lowest precedence). While
+    \c interval allows users to directly specify the STDMD time-step size
+    in seconds, in absence of \c interval, 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) can upset STDMD.
+  - Warning: In \c modeSorter algorithms of \c kiewat and \c kouZhang, a
+    function of \f$power(x,y)\f$ exists 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::DMD
+  - Foam::DMDModel
+
+SourceFiles
+    STDMD.C
+    STDMDTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef DMDModels_STDMD_H
+#define DMDModels_STDMD_H
+
+#include "DMDModel.H"
+#include "SquareMatrix.H"
+#include "RectangularMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace DMDModels
+{
+
+/*---------------------------------------------------------------------------*\
+                              Class STDMD Declaration
+\*---------------------------------------------------------------------------*/
+
+class STDMD
+:
+    public DMDModel
+{
+    typedef SquareMatrix<scalar> SMatrix;
+    typedef RectangularMatrix<scalar> RMatrix;
+    typedef RectangularMatrix<complex> RCMatrix;
+
+
+    // Private Enumerations
+
+        //- Options for the mode-sorting algorithm
+        enum modeSorterType : char
+        {
+            FIRST_SNAPSHOT = 0, //!< "First-snapshot amplitude magnitude method"
+            KIEWAT,            //!< "Modified weighted-amplitude scaling method"
+            KOU_ZHANG          //!< "Weighted-amplitude scaling method"
+        };
+
+        //- Names for modeSorterType
+        static const Enum<modeSorterType> modeSorterTypeNames;
+
+
+    // Private Data
+
+        //- Mode-sorting algorithm
+        enum modeSorterType modeSorter_;
+
+        //- Accumulated-in-time unitary similarity matrix produced by the
+        //- orthonormal decomposition of the augmented snapshot matrix 'z'
+        //  (K:Eq. 60)
+        RMatrix Q_;
+
+        //- Accumulated-in-time (squared) upper triangular matrix produced by
+        //- the orthonormal decomposition of the augmented snapshot matrix 'z'
+        //  (K:Eq. 64)
+        SMatrix G_;
+
+        //- Upper half of 'Q'
+        RMatrix Qupper_;
+
+        //- Lower half of 'Q'
+        RMatrix Qlower_;
+
+        //- Moore-Penrose pseudo-inverse of 'R' produced by
+        //- the QR decomposition of the last time-step 'Q'
+        RMatrix RxInv_;
+
+        //- Eigenvalues of modes
+        List<complex> evals_;
+
+        //- Eigenvectors of modes
+        RCMatrix evecs_;
+
+        //- (Non-negative) frequencies of modes
+        List<scalar> freqs_;
+
+        //- Indices of 'freqs' where frequencies are
+        //- non-negative and within ['fMin', 'fMax']
+        DynamicList<label> freqsi_;
+
+        //- Amplitudes of modes
+        List<complex> amps_;
+
+        //- (Descending) magnitudes of (complex) amplitudes of modes
+        List<scalar> mags_;
+
+        //- Indices of 'mags'
+        List<label> magsi_;
+
+        //- Name of operand field
+        const word fieldName_;
+
+        //- Name of operand patch
+        const word patch_;
+
+        //- First-processed snapshot required by the mode-sorting
+        //- algorithms at the final output computations (K:p. 43)
+        word timeName0_;
+
+        //- Min value to execute orthogonal basis expansion of 'Q' and 'G'
+        scalar minBasis_;
+
+        //- Min eigenvalue magnitude below where 'evals' are omitted
+        scalar minEval_;
+
+        //- STDMD time-step size that equals to
+        //- (executeInterval of DMD)*(deltaT of simulation) [s]
+        scalar dt_;
+
+        //- Maximum allowable magnitude for mag(eigenvalue)*(number of
+        //- STDMD steps) to be used in modeSorter={kiewat,kouZhang} to
+        //- avoid overflow errors
+        scalar sortLimiter_;
+
+        //- Min frequency: Output only entries corresponding to freqs >= 'fMin'
+        label fMin_;
+
+        //- Max frequency: Output only entries corresponding to freqs <= 'fMax'
+        label fMax_;
+
+        //- Number of maximum iterations of the classical
+        //- Gram-Schmidt procedure for the orthonormalisation
+        label nGramSchmidt_;
+
+        //- Maximum allowable matrix column for 'Q' and 'G'
+        //  'Q' is assumed to always have full-rank, thus 'Q.n() = rank'
+        label maxRank_;
+
+        //- Current execution-step index of STDMD,
+        //- not necessarily that of the simulation
+        label step_;
+
+        //- Number of output modes within input frequency range
+        //- starting from the most energetic mode
+        label nModes_;
+
+        //- Number of processors at each agglomeration unit
+        //- during the computation of reduced Koopman operator
+        label nAgglomerationProcs_;
+
+        //- (Internal) Flag to tag snapshots which are effectively empty
+        bool empty_;
+
+
+    // Private Member Functions
+
+        // Process
+
+            //- Return (parallel) L2-norm of a given column vector
+            scalar L2norm(const RMatrix& z) const;
+
+            //- Execute (parallel) classical Gram-Schmidt
+            //- process to orthonormalise 'ez' (Ka:Fig. 5)
+            RMatrix orthonormalise(RMatrix ez) const;
+
+            //- Expand orthonormal bases 'Q' and 'G' by stacking a column
+            //- '(ez/ezNorm)' to 'Q', and a row (Zero) and column (Zero)
+            //- to 'G' if '(ezNorm/zNorm > minBasis)'
+            void expand(const RMatrix& ez, const scalar ezNorm);
+
+            //- Compress orthonormal basis for 'Q' and 'G' if '(Q.n()>maxRank)'
+            void compress();
+
+
+        // Evaluation
+
+            //- Return reduced Koopman operator 'Atilde' (K:Eq. 78)
+            //  Also fills 'RxInv'.
+            //  The function was not divided into subsections to ensure
+            //  minimal usage of memory, hence the complexity of the function
+            SMatrix reducedKoopmanOperator();
+
+            //- Compute eigenvalues and eigenvectors of 'Atilde'
+            //  Remove any eigenvalues whose magnitude is smaller than
+            //  'minEVal' while keeping the order of elements the same
+            bool eigendecomposition(SMatrix& Atilde);
+
+            //- Compute and filter frequencies and its indices (K:Eq. 81)
+            //  Indices of 'freqs' where frequencies are
+            //  non-negative and within ['fMin', 'fMax']
+            void frequencies();
+
+            //- Compute amplitudes
+            void amplitudes();
+
+            //- Compute magnitudes and ordered magnitude indices
+            //  Includes advanced sorting methods using
+            //  the chosen weighted amplitude scaling method
+            void magnitudes();
+
+            //- Eigenvalue weighted amplitude scaling (KZ:Eq. 33)
+            //- Modified eigenvalue weighted amplitude scaling (K)
+            scalar sorter
+            (
+                const List<scalar>& weight,
+                const complex& amplitude,
+                const complex& eigenvalue,
+                const scalar modeNorm
+            ) const;
+
+            //- Compute and write mode dynamics
+            virtual bool dynamics();
+
+            //- Compute and write modes
+            virtual bool modes();
+
+            //- Select field type for modes
+            template<class Type>
+            bool modes();
+
+            //- Compute modes based on input field type
+            template<class GeoFieldType>
+            bool calcModes();
+
+            //- Compute a mode for a given scalar-based input field
+            template<class GeoFieldType>
+            typename std::enable_if
+            <
+                std::is_same<scalar, typename GeoFieldType::value_type>::value,
+                void
+            >::type calcMode
+            (
+                GeoFieldType& modeRe,
+                GeoFieldType& modeIm,
+                const RMatrix& primitiveMode,
+                const label i
+            );
+
+            //- Compute a mode for a given non-scalar-based input field
+            template<class GeoFieldType>
+            typename std::enable_if
+            <
+                !std::is_same<scalar, typename GeoFieldType::value_type>::value,
+                void
+            >::type calcMode
+            (
+                GeoFieldType& modeRe,
+                GeoFieldType& modeIm,
+                const RMatrix& primitiveMode,
+                const label i
+            );
+
+
+        // IO
+
+            //- Write objects of dynamics
+            void writeToFile(const word& fileName) const;
+
+            // Notifying the compiler that we want both writeToFile functions
+            using Foam::functionObjects::writeFile::writeToFile;
+
+            //- Write file header information
+            void writeFileHeader(Ostream& os) const;
+
+            //- Filter objects of dynamics according to 'freqsi' and 'magsi'
+            void filter();
+
+            //- 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
+            );
+
+
+public:
+
+    //- Runtime type information
+    TypeName("STDMD");
+
+
+    // Constructors
+
+        //- Construct from components
+        STDMD
+        (
+            const fvMesh& mesh,
+            const word& name,
+            const dictionary& dict
+        );
+
+        //- No copy construct
+        STDMD(const STDMD&) = delete;
+
+        //- No copy assignment
+        void operator=(const STDMD&) = delete;
+
+
+    //- Destructor
+    virtual ~STDMD() = default;
+
+
+    // Member Functions
+
+        // Evaluation
+
+            //- Initialise 'Q' and 'G' (both require the first two snapshots)
+            virtual bool initialise(const RMatrix& z);
+
+            //- Incremental orthonormal basis update (K:Fig. 15)
+            virtual bool update(const RMatrix& z);
+
+            //- Compute and write modes and
+            //- mode dynamics of model data members
+            virtual bool fit();
+
+
+        // IO
+
+            //- Read STDMD settings
+            virtual bool read(const dictionary& dict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace DMDModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "STDMDTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/DMD/DMDModels/derived/STDMD/STDMDTemplates.C b/src/functionObjects/field/DMD/DMDModels/derived/STDMD/STDMDTemplates.C
new file mode 100644
index 0000000000..408f2075fe
--- /dev/null
+++ b/src/functionObjects/field/DMD/DMDModels/derived/STDMD/STDMDTemplates.C
@@ -0,0 +1,236 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020-2021 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"
+#include "fixedValueFvPatchField.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+void Foam::DMDModels::STDMD::filterIndexed
+(
+    List<Type>& lst,
+    const UList<label>& indices
+)
+{
+    // Elements within [a, b]
+    List<Type> lstWithin(indices.size());
+
+    // Copy if frequency of element 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::DMDModels::STDMD::filterIndexed
+(
+    MatrixType& mat,
+    const UList<label>& indices
+)
+{
+    // Elements within [a, b]
+    MatrixType matWithin(labelPair(mat.m(), indices.size()));
+
+    // Copy if frequency of element is within [a, b]
+    label j = 0;
+    for (const auto& i : indices)
+    {
+        matWithin.subColumn(j) = mat.subColumn(i);
+        ++j;
+    }
+    mat.transfer(matWithin);
+}
+
+
+template<class Type>
+bool Foam::DMDModels::STDMD::modes()
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
+    typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
+
+    if (mesh_.foundObject<VolFieldType>(fieldName_))
+    {
+        return calcModes<VolFieldType>();
+    }
+    else if (mesh_.foundObject<SurfaceFieldType>(fieldName_))
+    {
+        return calcModes<SurfaceFieldType>();
+    }
+
+    return false;
+}
+
+
+template<class GeoFieldType>
+bool Foam::DMDModels::STDMD::calcModes()
+{
+    typedef typename GeoFieldType::value_type Type;
+
+    // Resize the number of output variables to "nModes" if requested
+    if (magsi_.size() > nModes_)
+    {
+        magsi_.resize(nModes_);
+    }
+
+    // Compute and write modes one by one
+    const RMatrix primitiveMode(Qupper_*RxInv_);
+    Qupper_.clear();
+    RxInv_.clear();
+
+    label modei = 0;
+    for (const label i : magsi_)
+    {
+        GeoFieldType modeRe
+        (
+            IOobject
+            (
+                "modeRe_" + std::to_string(modei)
+                + "_" + fieldName_ + "_" + name_,
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh_,
+            dimensioned<Type>(dimless, Zero),
+            fixedValueFvPatchField<Type>::typeName
+        );
+
+        GeoFieldType modeIm
+        (
+            IOobject
+            (
+                "modeIm_" + std::to_string(modei)
+                + "_" + fieldName_ + "_" + name_,
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh_,
+            dimensioned<Type>(dimless, Zero),
+            fixedValueFvPatchField<Type>::typeName
+        );
+
+        if (modeRe.size() != 0 && !empty_)
+        {
+            if (patch_.empty())
+            {
+                auto& inModeRe = modeRe.primitiveFieldRef();
+                auto& inModeIm = modeIm.primitiveFieldRef();
+
+                calcMode(inModeRe, inModeIm, primitiveMode, i);
+            }
+            else
+            {
+                const label patchi = mesh_.boundaryMesh().findPatchID(patch_);
+
+                auto& bfModeRe = modeRe.boundaryFieldRef()[patchi];
+                auto& bfModeIm = modeIm.boundaryFieldRef()[patchi];
+
+                calcMode(bfModeRe, bfModeIm, primitiveMode, i);
+            }
+        }
+
+        modeRe.write();
+        modeIm.write();
+        ++modei;
+    }
+
+    return true;
+}
+
+
+template<class GeoFieldType>
+typename std::enable_if
+<
+    std::is_same<Foam::scalar, typename GeoFieldType::value_type>::value,
+    void
+>::type Foam::DMDModels::STDMD::calcMode
+(
+    GeoFieldType& modeRe,
+    GeoFieldType& modeIm,
+    const RMatrix& primitiveMode,
+    const label i
+)
+{
+    const label fieldSize = modeRe.size();
+
+    for (label p = 0; p < primitiveMode.m(); ++p)
+    {
+        complex mode(Zero);
+        for (label q = 0; q < evecs_.m(); ++q)
+        {
+            mode += primitiveMode(p, q)*evecs_(q, i);
+        }
+        label p1 = p%fieldSize;
+        modeRe[p1] = mode.real();
+        modeIm[p1] = mode.imag();
+    }
+}
+
+
+template<class GeoFieldType>
+typename std::enable_if
+<
+    !std::is_same<Foam::scalar, typename GeoFieldType::value_type>::value,
+    void
+>::type Foam::DMDModels::STDMD::calcMode
+(
+    GeoFieldType& modeRe,
+    GeoFieldType& modeIm,
+    const RMatrix& primitiveMode,
+    const label i
+)
+{
+    const label fieldSize = modeRe.size();
+
+    for (label p = 0; p < primitiveMode.m(); ++p)
+    {
+        complex mode(Zero);
+        for (label q = 0; q < evecs_.m(); ++q)
+        {
+            mode += primitiveMode(p, q)*evecs_(q, i);
+        }
+        label p1 = p%fieldSize;
+        label p2 = p/fieldSize;
+        modeRe[p1][p2] = mode.real();
+        modeIm[p1][p2] = mode.imag();
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/STDMD/STDMDTemplates.C b/src/functionObjects/field/DMD/DMDTemplates.C
similarity index 51%
rename from src/functionObjects/field/STDMD/STDMDTemplates.C
rename to src/functionObjects/field/DMD/DMDTemplates.C
index e32bb2ddc3..a08fc5c02e 100644
--- a/src/functionObjects/field/STDMD/STDMDTemplates.C
+++ b/src/functionObjects/field/DMD/DMDTemplates.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -30,112 +30,83 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-template<class GeoFieldType>
-bool Foam::functionObjects::STDMD::getSnapshot()
+template<class Type>
+bool Foam::functionObjects::DMD::getSnapshot()
 {
-    if (!initialised_)
+    typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
+    typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
+
+    if (foundObject<VolFieldType>(fieldName_))
     {
-        init();
+        return getSnapshotField<VolFieldType>();
+    }
+    else if (foundObject<SurfaceFieldType>(fieldName_))
+    {
+        return getSnapshotField<SurfaceFieldType>();
     }
 
-    // Move previous-time snapshot into previous-time slot in z_
-    // Effectively moves the lower half of z_ to its upper half
+    return false;
+}
+
+
+template<class GeoFieldType>
+bool Foam::functionObjects::DMD::getSnapshotField()
+{
+    if (step_ == 0)
+    {
+        initialise();
+    }
+
+    if (z_.size() == 1)
+    {
+        return true;
+    }
+
+    // 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();
+    // Copy new current-time snapshot into current-time slot in "z"
+    // Effectively copies the new field elements into the lower half of "z"
+    const label nComps =
+        pTraits<typename GeoFieldType::value_type>::nComponents;
 
-    for (direction dir = 0; dir < nComps_; ++dir)
+    const GeoFieldType& field = lookupObject<GeoFieldType>(fieldName_);
+
+    if (patch_.empty())
     {
-        z_.subColumn(0, nSnap_ + dir*nField, nField) = Field.component(dir);
+        const label nField = field.size();
+
+        for (direction dir = 0; dir < nComps; ++dir)
+        {
+            z_.subColumn(0, nSnap_ + dir*nField, nField) = field.component(dir);
+        }
+    }
+    else
+    {
+        const label patchi = mesh_.boundaryMesh().findPatchID(patch_);
+
+        if (patchi < 0)
+        {
+            FatalErrorInFunction
+                << "Cannot find patch " << patch_
+                << exit(FatalError);
+        }
+
+        const typename GeoFieldType::Boundary& bf = field.boundaryField();
+
+        const Field<typename GeoFieldType::value_type>& pbf = bf[patchi];
+
+        const label nField = pbf.size();
+
+        for (direction dir = 0; dir < nComps; ++dir)
+        {
+            z_.subColumn(0, nSnap_ + dir*nField, nField) = pbf.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/src/functionObjects/field/Make/files b/src/functionObjects/field/Make/files
index 2b84b1f263..1ca2168217 100644
--- a/src/functionObjects/field/Make/files
+++ b/src/functionObjects/field/Make/files
@@ -121,6 +121,9 @@ stabilityBlendingFactor/stabilityBlendingFactor.C
 
 interfaceHeight/interfaceHeight.C
 
-STDMD/STDMD.C
+DMD/DMD.C
+DMD/DMDModels/DMDModel/DMDModel.C
+DMD/DMDModels/DMDModel/DMDModelNew.C
+DMD/DMDModels/derived/STDMD/STDMD.C
 
 LIB = $(FOAM_LIBBIN)/libfieldFunctionObjects
diff --git a/src/functionObjects/field/Make/options b/src/functionObjects/field/Make/options
index 57ea31c1cf..93f082e386 100644
--- a/src/functionObjects/field/Make/options
+++ b/src/functionObjects/field/Make/options
@@ -1,5 +1,6 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
     -I$(LIB_SRC)/fileFormats/lnInclude \
     -I$(LIB_SRC)/surfMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
@@ -22,6 +23,7 @@ EXE_INC = \
 
 LIB_LIBS = \
     -lfiniteVolume \
+    -lfiniteArea \
     -lfileFormats \
     -lsurfMesh \
     -lmeshTools \
diff --git a/src/functionObjects/field/STDMD/STDMD.C b/src/functionObjects/field/STDMD/STDMD.C
deleted file mode 100644
index 4fc192d670..0000000000
--- a/src/functionObjects/field/STDMD/STDMD.C
+++ /dev/null
@@ -1,1280 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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 (const int i : Pstream::allProcs())
-            {
-                const 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 label i : labelRange(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
deleted file mode 100644
index 42d137b970..0000000000
--- a/src/functionObjects/field/STDMD/STDMD.H
+++ /dev/null
@@ -1,645 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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 \c \<Type\>=Scalar/Vector/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
-
-// ************************************************************************* //