Commit 730233cd authored by Kutalmis Bercin's avatar Kutalmis Bercin Committed by Andrew Heather
Browse files

ENH: add new FO Streaming-Total Dynamic Mode Decomposition (STDMD)

    STDMD (i.e. Streaming Total Dynamic Mode Decomposition) is a variant of
    a data-driven dimensionality reduction method.

    STDMD is being used as a mathematical post-processing tool to compute
    a set of dominant modes out of a given flow (or dataset) each of which is
    associated with a constant frequency and decay rate, so that dynamic
    features of a given flow may become interpretable, and tractable.
    Among other Dynamic Mode Decomposition (DMD) variants, STDMD is presumed
    to provide the general DMD method capabilities alongside economised and
    feasible memory and CPU usage.

    Please refer to the header file documentation for further details.

  ENH: add new STDMD tutorial, pimpleFoam/laminar/cylinder2D
parent ef9ee7a8
......@@ -118,4 +118,6 @@ stabilityBlendingFactor/stabilityBlendingFactor.C
interfaceHeight/interfaceHeight.C
STDMD/STDMD.C
LIB = $(FOAM_LIBBIN)/libfieldFunctionObjects
EXE_INC = \
-I$(LIB_SRC)/OpenFOAM/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/fileFormats/lnInclude \
-I$(LIB_SRC)/surfMesh/lnInclude \
......@@ -22,6 +23,7 @@ EXE_INC = \
-I$(LIB_SRC)/phaseSystemModels/reactingEulerFoam/phaseSystems/lnInclude
LIB_LIBS = \
-lOpenFOAM \
-lfiniteVolume \
-lfileFormats \
-lsurfMesh \
......
This diff is collapsed.
This diff is collapsed.
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "volFields.H"
#include "surfaceFields.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class GeoFieldType>
bool Foam::functionObjects::STDMD::getSnapshot()
{
if (!initialised_)
{
init();
}
// Move previous-time snapshot into previous-time slot in z_
// Effectively moves the lower half of z_ to its upper half
std::rotate(z_.begin(), z_.begin() + nSnap_, z_.end());
// Copy new current-time snapshot into current-time slot in z_
// Effectively copies the new field elements into the lower half of z_
const GeoFieldType& Field = lookupObject<GeoFieldType>(fieldName_);
const label nField = Field.size();
for (direction dir = 0; dir < nComps_; ++dir)
{
z_.subColumn(0, nSnap_ + dir*nField, nField) = Field.component(dir);
}
return true;
}
template<class Type>
bool Foam::functionObjects::STDMD::getSnapshotType()
{
typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
if (foundObject<VolFieldType>(fieldName_))
{
return getSnapshot<VolFieldType>();
}
else if (foundObject<SurfaceFieldType>(fieldName_))
{
return getSnapshot<SurfaceFieldType>();
}
return false;
}
template<class Type>
bool Foam::functionObjects::STDMD::getComps()
{
typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
if (foundObject<VolFieldType>(fieldName_))
{
nComps_ = pTraits<typename VolFieldType::value_type>::nComponents;
return true;
}
else if (foundObject<SurfaceFieldType>(fieldName_))
{
nComps_ = pTraits<typename SurfaceFieldType::value_type>::nComponents;
return true;
}
return false;
}
template<class Type>
void Foam::functionObjects::STDMD::filterIndexed
(
List<Type>& lst,
const UList<label>& indices
)
{
// Elems within [a, b]
List<Type> lstWithin(indices.size());
// Copy if frequency of elem is within [a, b]
label j = 0;
for (const auto& i : indices)
{
lstWithin[j] = lst[i];
++j;
}
lst.transfer(lstWithin);
}
template<class MatrixType>
void Foam::functionObjects::STDMD::filterIndexed
(
MatrixType& mat,
const UList<label>& indices
)
{
// Elems within [a, b]
MatrixType matWithin(labelPair(mat.m(), indices.size()));
// Copy if frequency of elem is within [a, b]
label j = 0;
for (const auto& i : indices)
{
matWithin.subColumn(j) = mat.subColumn(i);
++j;
}
mat.transfer(matWithin);
}
// ************************************************************************* //
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v1912 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0.012984 0 0);
boundaryField
{
inlet
{
type fixedValue;
value $internalField;
}
outlet
{
type zeroGradient;
}
topAndBottom
{
type symmetry;
}
"left|right"
{
type empty;
}
cylinder
{
type noSlip;
}
}
// ************************************************************************* //
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v1912 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
{
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value $internalField;
}
topAndBottom
{
type symmetry;
}
"left|right"
{
type empty;
}
cylinder
{
type zeroGradient;
}
}
// ************************************************************************* //
#!/bin/sh
cd "${0%/*}" || exit # Run from this directory
. ${WM_PROJECT_DIR:?}/bin/tools/CleanFunctions # Tutorial clean functions
#------------------------------------------------------------------------------
cleanCase0
rm -rf dynamicCode
rm -rf constant/coarseMesh
rm -f system/coarseMesh/decomposeParDict
#------------------------------------------------------------------------------
#!/bin/sh
cd "${0%/*}" || exit # Run from this directory
. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions # Tutorial run functions
#------------------------------------------------------------------------------
./Allrun.pre
runApplication decomposePar
runParallel renumberMesh -overwrite
cp -f system/decomposeParDict system/coarseMesh
runApplication -s coarseMesh decomposePar -region coarseMesh
runParallel -s coarseMesh renumberMesh -overwrite -region coarseMesh
runParallel $(getApplication)
runParallel -s main redistributePar -reconstruct -latestTime
runParallel -s coarseMesh redistributePar -reconstruct \
-region coarseMesh -time '50,100,200'
#------------------------------------------------------------------------------
#!/bin/sh
cd "${0%/*}" || exit # Run from this directory
. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions # Tutorial run functions
#------------------------------------------------------------------------------
runApplication -s coarseMesh blockMesh -dict system/blockMeshDict.coarse
runApplication snappyHexMesh -overwrite
mkdir -p constant/coarseMesh
mv -f constant/polyMesh constant/coarseMesh
runApplication -s main blockMesh -dict system/blockMeshDict.main
runApplication mirrorMesh -overwrite
restore0Dir
#------------------------------------------------------------------------------
\ No newline at end of file
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v1912 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
transportModel Newtonian;
nu 1.558e-5;
// ************************************************************************* //
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v1912 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object turbulenceProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
simulationType laminar;
// ************************************************************************* //
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v1912 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
scale 1;
x0 -0.16;
x1 0.8;
y0 -0.24;
y1 0.24;
z0 -0.00375;
z1 0.00375;
nx 90;
ny 60;
nz 1;
vertices
(
($x0 $y0 $z0)
($x1 $y0 $z0)
($x1 $y1 $z0)
($x0 $y1 $z0)
($x0 $y0 $z1)
($x1 $y0 $z1)
($x1 $y1 $z1)
($x0 $y1 $z1)
);
blocks
(
hex (0 1 2 3 4 5 6 7) ($nx $ny $nz) simpleGrading (1 1 1)
);
edges
(
);
boundary
(
allPatches
{
type empty;
inGroups (allPatches);
faces
(
(3 7 6 2)
(1 5 4 0)
(0 4 7 3)
(2 6 5 1)
(0 3 2 1)
(4 5 6 7)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v1912 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
scale 1;
R 0.06; // Cylinder radius
x0 #eval{ 0.5*$R };
x2 #eval{ 2.38732*$R };
x3 #eval{ 10.0*$R };
xOutlet #eval{ 18.6667*$R };
xInlet #eval{ -10.125*$R };
RsinPi8 #eval{ $R*sin(0.125*pi()) };
RsinPi8n #eval{ -$RsinPi8 };
RcosPi8 #eval{ $R*cos(0.125*pi()) };
RcosPi8n #eval{ -$RcosPi8 };
RsinPi4 #eval{ $R*sin(0.25*pi()) };
x2sinPi8 #eval{ $x2*sin(0.125*pi()) };
x2sinPi8n #eval{ -$x2sinPi8 };
x2cosPi8 #eval{ $x2*cos(0.125*pi()) };
x2cosPi8n #eval{ -$x2cosPi8 };
x2sinPi4 #eval{ $x2*sin(0.25*pi()) };
z0 -0.0075;
z1 0.0075;
nz 1;
vertices #codeStream
{
codeInclude
#{
#include "pointField.H"
#};
code
#{
pointField points(19);
points[0] = point($R, 0, $z0);
points[1] = point($x2, 0, $z0);
points[2] = point($x3, 0, $z0);
points[3] = point($x3, $x2sinPi4, $z0);
points[4] = point($x2sinPi4, $x2sinPi4, $z0);
points[5] = point($RsinPi4, $RsinPi4, $z0);
points[6] = point($x3, $x3, $z0);
points[7] = point($x2sinPi4, $x3, $z0);
// Mirror +x points to -x side
points[11] = point(-points[0].x(), points[0].y(), points[0].z());
points[12] = point(-points[1].x(), points[1].y(), points[1].z());
points[13] = point(-points[2].x(), points[2].y(), points[2].z());
points[14] = point(-points[3].x(), points[3].y(), points[3].z());
points[15] = point(-points[4].x(), points[4].y(), points[4].z());
points[16] = point(-points[5].x(), points[5].y(), points[5].z());
points[17] = point(-points[6].x(), points[6].y(), points[6].z());
points[18] = point(-points[7].x(), points[7].y(), points[7].z());
// Points on the y-z-plane
points[8] = point(0, $x3, $z0);
points[9] = point(0, $x2, $z0);
points[10] = point(0, $R, $z0);
// Mirror -z points to +z side
label sz = points.size();
points.setSize(2*sz);
for (label i = 0; i < sz; ++i)
{
const point& pt = points[i];
points[i + sz] = point(pt.x(), pt.y(), $z1);
}
// Add an inner cylinder
sz = points.size();
label nAdd = 6;
points.setSize(sz + nAdd);
// Points within the inner cylinder
points[sz] = point(0, 0, $z0);
points[sz + 1] = point(0, $x0, $z0);
points[sz + 2] = point($x0, $x0, $z0);
points[sz + 3] = point($x0, 0, $z0);
// Mirror points from +x side to -x side
points[sz + 4] =
point(-points[sz + 2].x(), points[sz + 2].y(), points[sz + 2].z());
points[sz + 5] =
point(-points[sz + 3].x(), points[sz + 3].y(), points[sz + 3].z());
// Mirror -z points to +z side
sz = points.size();
points.setSize(sz + nAdd);
for (label i = 0; i < nAdd; ++i)
{
const point& pt = points[i+sz-nAdd];
points[i+sz] = point(pt.x(), pt.y(), $z1);
}
// Add downstream and upstream blocks
sz = points.size();
nAdd = 6;
points.setSize(sz + nAdd);
// Points on outlet
points[sz] = point($xOutlet, 0, $z0);
points[sz + 1] = point($xOutlet, $x3, $z0);
points[sz + 4] = point($xOutlet, $x2sinPi4, $z0);
// Points on inlet
points[sz + 2] = point($xInlet, 0, $z0);
points[sz + 3] = point($xInlet, $x3, $z0);
points[sz + 5] = point($xInlet, $x2sinPi4, $z0);
// Mirror -z points to +z side
sz = points.size();
points.setSize(sz + nAdd);
for (label i = 0; i < nAdd; ++i)
{
const point& pt = points[i + sz - nAdd];
points[i + sz] = point(pt.x(), pt.y(), $z1);
}
os << points;
#};
};
blocks
(
hex ( 5 4 9 10 24 23 28 29) (16 15 $nz) simpleGrading (1.6 1 1)
hex ( 0 1 4 5 19 20 23 24) (16 15 $nz) simpleGrading (1.6 1 1)
hex ( 1 2 3 4 20 21 22 23) (61 15 $nz) simpleGrading (1 1 1)
hex ( 4 3 6 7 23 22 25 26) (61 61 $nz) simpleGrading (1 1 1)
hex ( 9 4 7 8 28 23 26 27) (15 61 $nz) simpleGrading (1 1 1)
hex (15 16 10 9 34 35 29 28) (16 15 $nz) simpleGrading (0.625 1 1)
hex (12 11 16 15 31 30 35 34) (16 15 $nz) simpleGrading (0.625 1 1)
hex (13 12 15 14 32 31 34 33) (61 15 $nz) simpleGrading (1 1 1)
hex (14 15 18 17 33 34 37 36) (61 61 $nz) simpleGrading (1 1 1)
hex (15 9 8 18 34 28 27 37) (15 61 $nz) simpleGrading (1 1 1)
hex (2 50 54 3 21 56 60 22) (69 15 $nz) simpleGrading (1 1 1) // downstream
hex (3 54 51 6 22 60 57 25) (69 61 $nz) simpleGrading (1 1 1)