STDMD.H
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2020-2021 OpenCFD Ltd.
9-------------------------------------------------------------------------------
10License
11 This file is part of OpenFOAM.
12
13 OpenFOAM is free software: you can redistribute it and/or modify it
14 under the terms of the GNU General Public License as published by
15 the Free Software Foundation, either version 3 of the License, or
16 (at your option) any later version.
17
18 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 for more details.
22
23 You should have received a copy of the GNU General Public License
24 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25
26Class
27 Foam::DMDModels::STDMD
28
29Description
30 Streaming Total Dynamic Mode Decomposition (i.e. \c STDMD)
31 is a variant of dynamic mode decomposition.
32
33 Among other Dynamic Mode Decomposition (DMD) variants, \c STDMD is
34 presumed to provide the general DMD method capabilities alongside
35 economised and feasible memory and CPU usage.
36
37 The code implementation corresponds to Figs. 15-16 of the first citation
38 below, more broadly to Section 2.4.
39
40 References:
41 \verbatim
42 DMD and mode-sorting algorithms (tags:K, HRDC, KZ, HWR):
43 Kiewat, M. (2019).
44 Streaming modal decomposition approaches for vehicle aerodynamics.
45 PhD thesis. Munich: Technical University of Munich.
46 URL:mediatum.ub.tum.de/doc/1482652/1482652.pdf
47
48 Hemati, M. S., Rowley, C. W.,
49 Deem, E. A., & Cattafesta, L. N. (2017).
50 De-biasing the dynamic mode decomposition
51 for applied Koopman spectral analysis of noisy datasets.
52 Theoretical and Computational Fluid Dynamics, 31(4), 349-368.
53 DOI:10.1007/s00162-017-0432-2
54
55 Kou, J., & Zhang, W. (2017).
56 An improved criterion to select
57 dominant modes from dynamic mode decomposition.
58 European Journal of Mechanics-B/Fluids, 62, 109-129.
59 DOI:10.1016/j.euromechflu.2016.11.015
60
61 Hemati, M. S., Williams, M. O., & Rowley, C. W. (2014).
62 Dynamic mode decomposition for large and streaming datasets.
63 Physics of Fluids, 26(11), 111701.
64 DOI:10.1063/1.4901016
65
66 Parallel classical Gram-Schmidt process (tag:Ka):
67 Katagiri, T. (2003).
68 Performance evaluation of parallel
69 Gram-Schmidt re-orthogonalization methods.
70 In: Palma J. M. L. M., Sousa A. A., Dongarra J., Hernández V. (eds)
71 High Performance Computing for Computational Science — VECPAR 2002.
72 Lecture Notes in Computer Science, vol 2565, p. 302-314.
73 Berlin, Heidelberg: Springer.
74 DOI:10.1007/3-540-36569-9_19
75
76 Parallel direct tall-skinny QR decomposition (tags:BGD, DGHL):
77 Benson, A. R., Gleich, D. F., & Demmel, J. (2013).
78 Direct QR factorizations for
79 tall-and-skinny matrices in MapReduce architectures.
80 2013 IEEE International Conference on Big Data.
81 DOI:10.1109/bigdata.2013.6691583
82
83 Demmel, J., Grigori, L., Hoemmen, M., & Langou, J. (2012).
84 Communication-optimal parallel
85 and sequential QR and LU factorizations.
86 SIAM Journal on Scientific Computing, 34(1), A206-A239.
87 DOI:10.1137/080731992
88 \endverbatim
89
90Usage
91 Minimal example by using \c system/controlDict.functions:
92 \verbatim
93 DMD1
94 {
95 // Mandatory entries
96 DMDModel STDMD;
97
98 // Conditional mandatory entries
99
100 // Option-1
101 interval 5.5;
102
103 // Option-2
104 executeInterval 10;
105
106 // Optional entries
107 modeSorter kiewat;
108 nGramSchmidt 5;
109 maxRank 50;
110 nModes 50;
111 fMin 0;
112 fMax 1000000000;
113 nAgglomerationProcs 20;
114
115 // Optional entries (not recommended to change)
116 minBasis 0.00000001;
117 minEVal 0.00000001;
118 sortLimiter 500.0;
119
120 // Inherited entries
121 ...
122 }
123 \endverbatim
124
125 where the entries mean:
126 \table
127 Property | Description | Type | Reqd | Deflt
128 DMDModel | Type: STDMD | word | yes | -
129 interval | STDMD time-step size [s] | scalar | cndtnl <!--
130 --> | executeInterval*(current time-step of the simulation)
131 modeSorter | Mode-sorting algorithm model | word | no <!--
132 --> | firstSnapshot
133 nModes | Number of output modes in input frequency range <!--
134 --> | label | no | GREAT
135 maxRank | Max columns in accumulated matrices | label | no | GREAT
136 nGramSchmidt | Number of Gram-Schmidt iterations | label | no | 5
137 fMin | Min (non-negative) output frequency | label | no | 0
138 fMax | Max output frequency | label | no | GREAT
139 nAgglomerationProcs | Number of processors at each agglomeration <!--
140 --> unit during the computation of reduced Koopman <!--
141 --> operator | label | no | 20
142 minBasis | Orthogonal basis expansion threshold | scalar| no | 1e-8
143 minEVal | Min eigenvalue for below eigenvalues are omitted <!--
144 --> | scalar| no | 1e-8
145 sortLimiter | Max allowable magnitude for <!--
146 --> mag(eigenvalue)*(number of DMD steps) to be used in <!--
147 --> modeSorter={kiewat,kouZhang} to avoid overflow errors <!--
148 --> | scalar | no | 500.0
149 \endtable
150
151 Options for the \c modeSorter entry:
152 \verbatim
153 kiewat | Modified weighted-amplitude scaling method
154 kouZhang | Weighted-amplitude scaling method
155 firstSnapshot | First-snapshot amplitude magnitude method
156 \endverbatim
157
158Note
159 - To specify the STDMD time-step size (not necessarily equal to the
160 time step of the simulation), entries of either \c interval or
161 \c executeInterval must be available (highest to lowest precedence). While
162 \c interval allows users to directly specify the STDMD time-step size
163 in seconds, in absence of \c interval, for convenience,
164 \c executeInterval allows users to compute the STDMD time-step internally
165 by multiplying itself with the current time-step size of the simulation.
166 - Limitation: Restart is currently not available since intermediate writing
167 of STDMD matrices are not supported.
168 - Limitation: Non-physical input (e.g. full-zero fields) can upset STDMD.
169 - Warning: In \c modeSorter algorithms of \c kiewat and \c kouZhang, a
170 function of \f$power(x,y)\f$ exists where \c x is the magnitude of an
171 eigenvalue, and \c y is the number of STDMD snapshots. This function poses
172 a risk of overflow errors since, for example, if \c x ends up above 1.5 with
173 500 snapshots or more, this function automatically throws the floating point
174 error meaning overflow. Therefore, the domain-range of this function is
175 heuristically constrained by the optional entry \c sortLimiter which skips
176 the evaluation of this function for a given
177 mag(eigenvalue)*(no. of STDMD steps), i.e. x*y, whose magnitude is larger
178 than \c sortLimiter.
179
180See also
181 - Foam::functionObjects::DMD
182 - Foam::DMDModel
183
184SourceFiles
185 STDMD.C
186 STDMDTemplates.C
187
188\*---------------------------------------------------------------------------*/
189
190#ifndef DMDModels_STDMD_H
191#define DMDModels_STDMD_H
192
193#include "DMDModel.H"
194#include "SquareMatrix.H"
195#include "RectangularMatrix.H"
196
197// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
198
199namespace Foam
200{
201namespace DMDModels
202{
203
204/*---------------------------------------------------------------------------*\
205 Class STDMD Declaration
206\*---------------------------------------------------------------------------*/
207
208class STDMD
209:
210 public DMDModel
211{
212 typedef SquareMatrix<scalar> SMatrix;
213 typedef RectangularMatrix<scalar> RMatrix;
214 typedef RectangularMatrix<complex> RCMatrix;
215
216
217 // Private Enumerations
218
219 //- Options for the mode-sorting algorithm
220 enum modeSorterType : char
221 {
222 FIRST_SNAPSHOT = 0,
223 KIEWAT,
224 KOU_ZHANG
225 };
226
227 //- Names for modeSorterType
228 static const Enum<modeSorterType> modeSorterTypeNames;
229
230
231 // Private Data
232
233 //- Mode-sorting algorithm
234 enum modeSorterType modeSorter_;
235
236 //- Accumulated-in-time unitary similarity matrix produced by the
237 //- orthonormal decomposition of the augmented snapshot matrix 'z'
238 // (K:Eq. 60)
239 RMatrix Q_;
240
241 //- Accumulated-in-time (squared) upper triangular matrix produced by
242 //- the orthonormal decomposition of the augmented snapshot matrix 'z'
243 // (K:Eq. 64)
244 SMatrix G_;
245
246 //- Upper half of 'Q'
247 RMatrix Qupper_;
248
249 //- Lower half of 'Q'
250 RMatrix Qlower_;
251
252 //- Moore-Penrose pseudo-inverse of 'R' produced by
253 //- the QR decomposition of the last time-step 'Q'
254 RMatrix RxInv_;
255
256 //- Eigenvalues of modes
257 List<complex> evals_;
258
259 //- Eigenvectors of modes
260 RCMatrix evecs_;
261
262 //- (Non-negative) frequencies of modes
263 List<scalar> freqs_;
264
265 //- Indices of 'freqs' where frequencies are
266 //- non-negative and within ['fMin', 'fMax']
267 DynamicList<label> freqsi_;
268
269 //- Amplitudes of modes
270 List<complex> amps_;
271
272 //- (Descending) magnitudes of (complex) amplitudes of modes
273 List<scalar> mags_;
274
275 //- Indices of 'mags'
276 List<label> magsi_;
277
278 //- Names of operand patches
279 const wordRes patches_;
280
281 //- Name of operand field
282 const word fieldName_;
283
284 //- First-processed snapshot required by the mode-sorting
285 //- algorithms at the final output computations (K:p. 43)
286 word timeName0_;
287
288 //- Min value to execute orthogonal basis expansion of 'Q' and 'G'
289 scalar minBasis_;
290
291 //- Min eigenvalue magnitude below where 'evals' are omitted
292 scalar minEval_;
293
294 //- STDMD time-step size that equals to
295 //- (executeInterval of DMD)*(deltaT of simulation) [s]
296 scalar dt_;
297
298 //- Maximum allowable magnitude for mag(eigenvalue)*(number of
299 //- STDMD steps) to be used in modeSorter={kiewat,kouZhang} to
300 //- avoid overflow errors
301 scalar sortLimiter_;
302
303 //- Min frequency: Output only entries corresponding to freqs >= 'fMin'
304 label fMin_;
305
306 //- Max frequency: Output only entries corresponding to freqs <= 'fMax'
307 label fMax_;
308
309 //- Number of maximum iterations of the classical
310 //- Gram-Schmidt procedure for the orthonormalisation
311 label nGramSchmidt_;
312
313 //- Maximum allowable matrix column for 'Q' and 'G'
314 // 'Q' is assumed to always have full-rank, thus 'Q.n() = rank'
315 label maxRank_;
316
317 //- Current execution-step index of STDMD,
318 //- not necessarily that of the simulation
319 label step_;
320
321 //- Number of output modes within input frequency range
322 //- starting from the most energetic mode
323 label nModes_;
324
325 //- Number of processors at each agglomeration unit
326 //- during the computation of reduced Koopman operator
327 label nAgglomerationProcs_;
328
329 //- (Internal) Flag to tag snapshots which are effectively empty
330 bool empty_;
331
332
333 // Private Member Functions
334
335 // Evaluation
336
337 //- Return (parallel) L2-norm of a given column vector
338 scalar L2norm(const RMatrix& z) const;
339
340 //- Execute (parallel) classical Gram-Schmidt
341 //- process to orthonormalise 'ez' (Ka:Fig. 5)
342 RMatrix orthonormalise(RMatrix ez) const;
343
344 //- Expand orthonormal bases 'Q' and 'G' by stacking a column
345 //- '(ez/ezNorm)' to 'Q', and a row (Zero) and column (Zero)
346 //- to 'G' if '(ezNorm/zNorm > minBasis)'
347 void expand(const RMatrix& ez, const scalar ezNorm);
348
349 //- Update 'G' before the potential orthonormal basis compression
350 void updateG(const RMatrix& z);
351
352 //- Compress orthonormal basis for 'Q' and 'G' if '(Q.n()>maxRank)'
353 void compress();
354
355 //- Return reduced Koopman operator 'Atilde' (K:Eq. 78)
356 // Also fills 'RxInv'.
357 // The function was not divided into subsections to ensure
358 // minimal usage of memory, hence the complexity of the function
359 SMatrix reducedKoopmanOperator();
360
361 //- Compute eigenvalues and eigenvectors of 'Atilde'
362 // Remove any eigenvalues whose magnitude is smaller than
363 // 'minEVal' while keeping the order of elements the same
364 bool eigendecomposition(SMatrix& Atilde);
365
366 //- Compute and filter frequencies and its indices (K:Eq. 81)
367 // Indices of 'freqs' where frequencies are
368 // non-negative and within ['fMin', 'fMax']
369 void frequencies();
370
371 //- Compute amplitudes
372 void amplitudes();
373
374 //- Compute magnitudes and ordered magnitude indices
375 // Includes advanced sorting methods using
376 // the chosen weighted amplitude scaling method
377 void magnitudes();
378
379 //- Eigenvalue weighted amplitude scaling (KZ:Eq. 33)
380 //- Modified eigenvalue weighted amplitude scaling (K)
381 scalar sorter
382 (
383 const List<scalar>& weight,
384 const complex& amplitude,
385 const complex& eigenvalue,
386 const scalar modeNorm
387 ) const;
388
389 //- Compute and write mode dynamics
390 virtual bool dynamics();
391
392 //- Compute and write modes
393 virtual bool modes();
394
395 //- Select field type for modes
396 template<class Type>
397 bool modes();
398
399 //- Compute modes based on input field type
400 template<class GeoFieldType>
401 bool calcModes();
402
403 //- Compute a mode for a given scalar-based input field
404 template<class GeoFieldType>
405 typename std::enable_if
406 <
407 std::is_same<scalar, typename GeoFieldType::value_type>::value,
408 void
409 >::type calcMode
410 (
411 GeoFieldType& modeRe,
412 GeoFieldType& modeIm,
413 const RMatrix& primitiveMode,
414 const label magi,
415 const label rowi = 0
416 );
417
418 //- Compute a mode for a given non-scalar-based input field
419 template<class GeoFieldType>
420 typename std::enable_if
421 <
422 !std::is_same<scalar, typename GeoFieldType::value_type>::value,
423 void
424 >::type calcMode
425 (
426 GeoFieldType& modeRe,
427 GeoFieldType& modeIm,
428 const RMatrix& primitiveMode,
429 const label magi,
430 const label rowi = 0
431 );
432
433
434 // I-O
435
436 //- Write objects of dynamics
437 void writeToFile(const word& fileName) const;
438
439 // Notifying the compiler that we want both writeToFile functions
441
442 //- Write file header information
443 void writeFileHeader(Ostream& os) const;
444
445 //- Filter objects of dynamics according to 'freqsi' and 'magsi'
446 void filter();
447
448 //- Return a new List containing elems of List at 'indices'
449 template<class Type>
450 void filterIndexed
451 (
452 List<Type>& lst,
453 const UList<label>& indices
454 );
455
456 //- Return a new Matrix containing columns of Matrix at 'indices'
457 template<class MatrixType>
458 void filterIndexed
459 (
460 MatrixType& lst,
461 const UList<label>& indices
462 );
463
464
465public:
466
467 //- Runtime type information
468 TypeName("STDMD");
469
470
471 // Constructors
472
473 //- Construct from components
474 STDMD
475 (
476 const fvMesh& mesh,
477 const word& name,
478 const dictionary& dict
479 );
480
481 //- No copy construct
482 STDMD(const STDMD&) = delete;
483
484 //- No copy assignment
485 void operator=(const STDMD&) = delete;
486
487
488 //- Destructor
489 virtual ~STDMD() = default;
490
491
492 // Member Functions
493
494 // Evaluation
495
496 //- Initialise 'Q' and 'G' (both require the first two snapshots)
497 virtual bool initialise(const RMatrix& z);
498
499 //- Incremental orthonormal basis update (K:Fig. 15)
500 virtual bool update(const RMatrix& z);
501
502 //- Compute and write modes and
503 //- mode dynamics of model data members
504 virtual bool fit();
505
506
507 // IO
508
509 //- Read STDMD settings
510 virtual bool read(const dictionary& dict);
511};
512
513
514// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
515
516} // End namespace DMDModels
517} // End namespace Foam
518
519// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
520
521#ifdef NoRepository
522 #include "STDMDTemplates.C"
523#endif
524
525// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
526
527#endif
528
529// ************************************************************************* //
Abstract base class for DMD models to handle DMD characteristics for the DMD function object.
Definition: DMDModel.H:68
Streaming Total Dynamic Mode Decomposition (i.e. STDMD) is a variant of dynamic mode decomposition.
Definition: STDMD.H:306
virtual bool read(const dictionary &dict)
Read STDMD settings.
Definition: STDMD.C:864
virtual bool initialise(const RMatrix &z)
Initialise 'Q' and 'G' (both require the first two snapshots)
Definition: STDMD.C:951
TypeName("STDMD")
Runtime type information.
void operator=(const STDMD &)=delete
No copy assignment.
virtual bool fit()
Definition: STDMD.C:1044
STDMD(const STDMD &)=delete
No copy construct.
virtual ~STDMD()=default
Destructor.
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
Enum is a wrapper around a list of names/values that represent particular enumeration (or int) values...
Definition: Enum.H:61
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: List.H:77
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
Definition: Ostream.H:62
A templated (M x N) rectangular matrix of objects of <Type>, containing M*N elements,...
A templated (N x N) square matrix of objects of <Type>, containing N*N elements, derived from Matrix.
Definition: SquareMatrix.H:66
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: UList.H:94
A complex number, similar to the C++ complex type.
Definition: complex.H:83
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
A class for handling file names.
Definition: fileName.H:76
virtual bool writeToFile() const
Flag to allow writing to file.
Definition: writeFile.C:250
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:91
A List of wordRe with additional matching capabilities.
Definition: wordRes.H:54
A class for handling words, derived from Foam::string.
Definition: word.H:68
mesh update()
dynamicFvMesh & mesh
OBJstream os(runTime.globalPath()/outputName)
Namespace for OpenFOAM.
fileName::Type type(const fileName &name, const bool followLink=true)
Return the file type: DIRECTORY or FILE, normally following symbolic links.
Definition: MSwindows.C:598
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
dictionary dict
#define TypeName(TypeNameString)
Declare a ClassName() with extra virtual type info.
Definition: typeInfo.H:73