diff --git a/src/functionObjects/field/DMD/DMD.C b/src/functionObjects/field/DMD/DMD.C
new file mode 100644
index 0000000000000000000000000000000000000000..e225d31cb484fe14bcfb7139007a36aba99b833a
--- /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 0000000000000000000000000000000000000000..b468e85f8648b73667729ade808b387c6b32f0cf
--- /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 0000000000000000000000000000000000000000..c43727548c83d8f42e0f75c289ba04cd807c04f4
--- /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 0000000000000000000000000000000000000000..8081db5a01a88d06e68acf89467e1ed4383ed7eb
--- /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 0000000000000000000000000000000000000000..e3a72a4a5c2967a68c101ec4fa504d7914a7d9b4
--- /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 0000000000000000000000000000000000000000..9281dd2385a779afffcf74b165463eaaa82b70e5
--- /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 0000000000000000000000000000000000000000..ac180feb7ac0f29ed629b541825ac96ba2de3854
--- /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 0000000000000000000000000000000000000000..0b86e137fc861b0a6ba7bf3f362eb1aa768ad930
--- /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 0000000000000000000000000000000000000000..408f2075fe2fd2e51294f43b01957b6a6cea5a92
--- /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 e32bb2ddc33ff8e8680071ae986b51e50a45826e..a08fc5c02ed69cb05c2906762784cb3e665f1051 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,111 +30,82 @@ License
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
-template<class GeoFieldType>
-bool Foam::functionObjects::STDMD::getSnapshot()
-{
- if (!initialised_)
- {
- init();
- }
-
- // Move previous-time snapshot into previous-time slot in z_
- // Effectively moves the lower half of z_ to its upper half
- std::rotate(z_.begin(), z_.begin() + nSnap_, z_.end());
-
- // Copy new current-time snapshot into current-time slot in z_
- // Effectively copies the new field elements into the lower half of z_
- const GeoFieldType& Field = lookupObject<GeoFieldType>(fieldName_);
- const label nField = Field.size();
-
- for (direction dir = 0; dir < nComps_; ++dir)
- {
- z_.subColumn(0, nSnap_ + dir*nField, nField) = Field.component(dir);
- }
-
- return true;
-}
-
-
template<class Type>
-bool Foam::functionObjects::STDMD::getSnapshotType()
+bool Foam::functionObjects::DMD::getSnapshot()
{
typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
if (foundObject<VolFieldType>(fieldName_))
{
- return getSnapshot<VolFieldType>();
+ return getSnapshotField<VolFieldType>();
}
else if (foundObject<SurfaceFieldType>(fieldName_))
{
- return getSnapshot<SurfaceFieldType>();
+ return getSnapshotField<SurfaceFieldType>();
}
return false;
}
-template<class Type>
-bool Foam::functionObjects::STDMD::getComps()
+template<class GeoFieldType>
+bool Foam::functionObjects::DMD::getSnapshotField()
{
- typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
- typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
-
- if (foundObject<VolFieldType>(fieldName_))
+ if (step_ == 0)
{
- nComps_ = pTraits<typename VolFieldType::value_type>::nComponents;
- return true;
+ initialise();
}
- else if (foundObject<SurfaceFieldType>(fieldName_))
+
+ if (z_.size() == 1)
{
- nComps_ = pTraits<typename SurfaceFieldType::value_type>::nComponents;
return true;
}
- return false;
-}
+ // 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 label nComps =
+ pTraits<typename GeoFieldType::value_type>::nComponents;
-template<class Type>
-void Foam::functionObjects::STDMD::filterIndexed
-(
- List<Type>& lst,
- const UList<label>& indices
-)
-{
- // Elems within [a, b]
- List<Type> lstWithin(indices.size());
+ const GeoFieldType& field = lookupObject<GeoFieldType>(fieldName_);
- // Copy if frequency of elem is within [a, b]
- label j = 0;
- for (const auto& i : indices)
+ if (patch_.empty())
{
- lstWithin[j] = lst[i];
- ++j;
+ const label nField = field.size();
+
+ for (direction dir = 0; dir < nComps; ++dir)
+ {
+ z_.subColumn(0, nSnap_ + dir*nField, nField) = field.component(dir);
+ }
}
- lst.transfer(lstWithin);
-}
+ else
+ {
+ const label patchi = mesh_.boundaryMesh().findPatchID(patch_);
+ if (patchi < 0)
+ {
+ FatalErrorInFunction
+ << "Cannot find patch " << patch_
+ << exit(FatalError);
+ }
-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()));
+ const typename GeoFieldType::Boundary& bf = field.boundaryField();
- // Copy if frequency of elem is within [a, b]
- label j = 0;
- for (const auto& i : indices)
- {
- matWithin.subColumn(j) = mat.subColumn(i);
- ++j;
+ 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);
+ }
}
- mat.transfer(matWithin);
+
+ return true;
}
diff --git a/src/functionObjects/field/Make/files b/src/functionObjects/field/Make/files
index 2b84b1f263356a10c2baf91aaf1c9c6973f009fc..1ca2168217c384887478b799111d13973fa2f704 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 57ea31c1cfb026a75c927836f43aed29ef24bf1d..93f082e386b2940c9517351ed9fb8f9bd156ded9 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 4fc192d670a6d399ae95c10633c283db3b0f8353..0000000000000000000000000000000000000000
--- 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 42d137b97062e718903944ccc7c529191048adb7..0000000000000000000000000000000000000000
--- 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
-
-// ************************************************************************* //