mirror of
https://develop.openfoam.com/Development/openfoam.git
synced 2025-11-28 03:28:01 +00:00
ENH: DMD: refactor dynamic mode decomposition FO
ENH: DMD: enable operations on patch fields ENH: DMD: enable postProcess utility
This commit is contained in:
Kutalmis Bercin
committed by
Andrew Heather
Andrew Heather
parent
8bfab7aa80
commit
2a8e1c0865
234
src/functionObjects/field/DMD/DMD.C
Normal file
234
src/functionObjects/field/DMD/DMD.C
Normal file
@ -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;
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
278
src/functionObjects/field/DMD/DMD.H
Normal file
278
src/functionObjects/field/DMD/DMD.H
Normal file
@ -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
|
||||
|
||||
// ************************************************************************* //
|
||||
78
src/functionObjects/field/DMD/DMDModels/DMDModel/DMDModel.C
Normal file
78
src/functionObjects/field/DMD/DMDModels/DMDModel/DMDModel.C
Normal file
@ -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;
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
203
src/functionObjects/field/DMD/DMDModels/DMDModel/DMDModel.H
Normal file
203
src/functionObjects/field/DMD/DMDModels/DMDModel/DMDModel.H
Normal file
@ -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
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -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));
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -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;
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
988
src/functionObjects/field/DMD/DMDModels/derived/STDMD/STDMD.C
Normal file
988
src/functionObjects/field/DMD/DMDModels/derived/STDMD/STDMD.C
Normal file
@ -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;
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
530
src/functionObjects/field/DMD/DMDModels/derived/STDMD/STDMD.H
Normal file
530
src/functionObjects/field/DMD/DMDModels/derived/STDMD/STDMD.H
Normal file
@ -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
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -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();
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -5,7 +5,7 @@
|
||||
\\ / A nd | www.openfoam.com
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
Copyright (C) 2020 OpenCFD Ltd.
|
||||
Copyright (C) 2020-2021 OpenCFD Ltd.
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
@ -30,112 +30,83 @@ License
|
||||
|
||||
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
|
||||
|
||||
template<class GeoFieldType>
|
||||
bool Foam::functionObjects::STDMD::getSnapshot()
|
||||
template<class Type>
|
||||
bool Foam::functionObjects::DMD::getSnapshot()
|
||||
{
|
||||
if (!initialised_)
|
||||
typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
|
||||
typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
|
||||
|
||||
if (foundObject<VolFieldType>(fieldName_))
|
||||
{
|
||||
init();
|
||||
return getSnapshotField<VolFieldType>();
|
||||
}
|
||||
else if (foundObject<SurfaceFieldType>(fieldName_))
|
||||
{
|
||||
return getSnapshotField<SurfaceFieldType>();
|
||||
}
|
||||
|
||||
// Move previous-time snapshot into previous-time slot in z_
|
||||
// Effectively moves the lower half of z_ to its upper half
|
||||
return false;
|
||||
}
|
||||
|
||||
|
||||
template<class GeoFieldType>
|
||||
bool Foam::functionObjects::DMD::getSnapshotField()
|
||||
{
|
||||
if (step_ == 0)
|
||||
{
|
||||
initialise();
|
||||
}
|
||||
|
||||
if (z_.size() == 1)
|
||||
{
|
||||
return true;
|
||||
}
|
||||
|
||||
// Move previous-time snapshot into previous-time slot in "z"
|
||||
// Effectively moves the lower half of "z" to its upper half
|
||||
std::rotate(z_.begin(), z_.begin() + nSnap_, z_.end());
|
||||
|
||||
// Copy new current-time snapshot into current-time slot in z_
|
||||
// Effectively copies the new field elements into the lower half of z_
|
||||
const GeoFieldType& Field = lookupObject<GeoFieldType>(fieldName_);
|
||||
const label nField = Field.size();
|
||||
// Copy new current-time snapshot into current-time slot in "z"
|
||||
// Effectively copies the new field elements into the lower half of "z"
|
||||
const label nComps =
|
||||
pTraits<typename GeoFieldType::value_type>::nComponents;
|
||||
|
||||
for (direction dir = 0; dir < nComps_; ++dir)
|
||||
const GeoFieldType& field = lookupObject<GeoFieldType>(fieldName_);
|
||||
|
||||
if (patch_.empty())
|
||||
{
|
||||
z_.subColumn(0, nSnap_ + dir*nField, nField) = Field.component(dir);
|
||||
const label nField = field.size();
|
||||
|
||||
for (direction dir = 0; dir < nComps; ++dir)
|
||||
{
|
||||
z_.subColumn(0, nSnap_ + dir*nField, nField) = field.component(dir);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
const label patchi = mesh_.boundaryMesh().findPatchID(patch_);
|
||||
|
||||
if (patchi < 0)
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "Cannot find patch " << patch_
|
||||
<< exit(FatalError);
|
||||
}
|
||||
|
||||
const typename GeoFieldType::Boundary& bf = field.boundaryField();
|
||||
|
||||
const Field<typename GeoFieldType::value_type>& pbf = bf[patchi];
|
||||
|
||||
const label nField = pbf.size();
|
||||
|
||||
for (direction dir = 0; dir < nComps; ++dir)
|
||||
{
|
||||
z_.subColumn(0, nSnap_ + dir*nField, nField) = pbf.component(dir);
|
||||
}
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
template<class Type>
|
||||
bool Foam::functionObjects::STDMD::getSnapshotType()
|
||||
{
|
||||
typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
|
||||
typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
|
||||
|
||||
if (foundObject<VolFieldType>(fieldName_))
|
||||
{
|
||||
return getSnapshot<VolFieldType>();
|
||||
}
|
||||
else if (foundObject<SurfaceFieldType>(fieldName_))
|
||||
{
|
||||
return getSnapshot<SurfaceFieldType>();
|
||||
}
|
||||
|
||||
return false;
|
||||
}
|
||||
|
||||
|
||||
template<class Type>
|
||||
bool Foam::functionObjects::STDMD::getComps()
|
||||
{
|
||||
typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
|
||||
typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
|
||||
|
||||
if (foundObject<VolFieldType>(fieldName_))
|
||||
{
|
||||
nComps_ = pTraits<typename VolFieldType::value_type>::nComponents;
|
||||
return true;
|
||||
}
|
||||
else if (foundObject<SurfaceFieldType>(fieldName_))
|
||||
{
|
||||
nComps_ = pTraits<typename SurfaceFieldType::value_type>::nComponents;
|
||||
return true;
|
||||
}
|
||||
|
||||
return false;
|
||||
}
|
||||
|
||||
|
||||
template<class Type>
|
||||
void Foam::functionObjects::STDMD::filterIndexed
|
||||
(
|
||||
List<Type>& lst,
|
||||
const UList<label>& indices
|
||||
)
|
||||
{
|
||||
// Elems within [a, b]
|
||||
List<Type> lstWithin(indices.size());
|
||||
|
||||
// Copy if frequency of elem is within [a, b]
|
||||
label j = 0;
|
||||
for (const auto& i : indices)
|
||||
{
|
||||
lstWithin[j] = lst[i];
|
||||
++j;
|
||||
}
|
||||
lst.transfer(lstWithin);
|
||||
}
|
||||
|
||||
|
||||
template<class MatrixType>
|
||||
void Foam::functionObjects::STDMD::filterIndexed
|
||||
(
|
||||
MatrixType& mat,
|
||||
const UList<label>& indices
|
||||
)
|
||||
{
|
||||
// Elems within [a, b]
|
||||
MatrixType matWithin(labelPair(mat.m(), indices.size()));
|
||||
|
||||
// Copy if frequency of elem is within [a, b]
|
||||
label j = 0;
|
||||
for (const auto& i : indices)
|
||||
{
|
||||
matWithin.subColumn(j) = mat.subColumn(i);
|
||||
++j;
|
||||
}
|
||||
mat.transfer(matWithin);
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -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
|
||||
|
||||
@ -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 \
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
@ -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
|
||||
|
||||
// ************************************************************************* //
|
||||
Reference in New Issue
Block a user