# Analyzing Nonlinear Dynamics via Data-Driven Dynamic Mode Decomposition-Like Methods

Authors: Soledad Le Clainche; José M. Vega
Journal: Complexity
 (2018)

---

## Abstract

This article presents a review on two methods based on dynamic mode decomposition and its multiple applications, focusing on higher order dynamic mode decomposition (which provides a purely temporal Fourier-like decomposition) and spatiotemporal Koopman decomposition (which gives a spatiotemporal Fourier-like decomposition). These methods are purely data-driven, using either numerical or experimental data, and permit reconstructing the given data and identifying the temporal growth rates and frequencies involved in the dynamics and the spatial growth rates and wavenumbers in the case of the spatiotemporal Koopman decomposition. Thus, they may be used to either identify and extrapolate the dynamics from transient behavior to permanent dynamics or construct efficient, purely data-driven reduced order models.


---

## Body

## 1. Introduction

Uncovering the quantitative essence of complex signals, either numerical or experimental, coming from nonlinear systems is an interesting topic indata science andfluid dynamics. Among its several applications in multidisciplinary sciences, the purely data-driven analysis of nonlinear dynamics has raised the interest of applied mathematicians (to understand the dynamics of ordinary or partial differential equations and systems, such as the complex Ginzburg-Landau equation), doctors (to find relevant patterns in medical images), and aeronautical, mechanical, and civil engineers (looking for flow patterns or global instabilities or computing aerodynamic forces on slender or bluff bodies). Even in simple, low-dimensional dynamical systems, such as the Lorenz system (3 ordinary differential equations), it is possible to find complex dynamics, including periodic, quasi-periodic, and chaotic behaviors, as well as complex transients.Fourier analysis has long been used as a data-processing tool for analyzing complex signals. In addition to improving the standard Fourier transform (e.g.,fast Fourier transform,power spectral density [1], and the so-calledLaskar method [2, 3]), several techniques have been developed in recent years that extend the scope of Fourier analysis by seeking not only the frequencies of the signal but also the growth rates, which are relevant to transient behavior. One of the most popular methods to obtain such an expansion isdynamic mode decomposition (DMD) [4] (see also [5–7]), which needs to be extended to give optimal results [8], as will be discussed in this paper. These DMD-like methods are based on the idea of decomposing a signal as an expansion in Fourier-like modes, as (1)vt≃∑n=1Nanuneδn+iωntfort1≤t≤t1+T,where v is a time-dependent vector signal, which can be either real or complex, un are conveniently normalized spatial modes, an>0 are the mode amplitudes, and δn and ωn are the associated growth rates and frequencies, respectively. Dynamics of the type (1) are paramount in applications, as it will be illustrated in the paper. A good example is the joint wake (accounting for interaction between the individual wakes of the various buildings) of urban topographies; see Figure 1. This example is interesting for two reasons. First, direct numerical simulation is not affordable in realistic conditions due to the huge value of the Reynolds number, meaning that numerical data must be obtained using turbulence modeling, which only captures the large scales. These are temporally quasi-periodic and may be described in the form (1); small scales, instead, exhibit spatiotemporal chaos. Experimental data, on the other hand, exhibit nonnegligible errors that must be filtered out and also describe the large scales only.Figure 1
An urban topography composed by three skyscrapers and four elongated buildings (left) and horizontal snapshots of the associated joint wake (right).Expansion (1) is appropriate to represent temporally growing (if δn>0) or decaying (if δn<0), periodic (if the frequencies ωn are commensurable) or quasi-periodic (if the frequencies are incommensurable) phenomena. Commensurability is a subtle matter in finite precision computations, where exact incommensurability is not possible. Instead, we shall guess that a set of frequencies are incommensurable when the computed values are not commensurable to some small tolerance. Namely, e.g., three frequencies are incommensurable if no relation of the form mω1±nω2±pω3=ε is possible for moderate values (say, not larger than 10) of the natural numbers m, n, and p and some small (say, ~10-5) ε. This can be discerned provided that the frequencies are computed with great accuracy. Moreover, the robustness and accuracy of the approximation (1) are a good means to ensure that the right periodic or quasi-periodic dynamics have been captured. In chaotic dynamics, instead, an infinite number of modes should be considered even for moderate accuracy and the analysis using these methods is subtle.In the data-driven methods described and used along the paper, space and time will both be discretized. However, it is enlightening to consider at the moment continuous time, with the vector snapshotsv of size J. The simplest problem giving the dynamics in the right hand side of (1) is the ODE linear system (2)dvtdt=Avt,for an appropriate J×J matrix A that may be computed via the pseudoinverse (applied to a discretized version of this equation, considering a collection of nearby snapshots and approximating the time derivative by finite differences).Obviously, (2) makes sense if the dynamics associated with the snapshots are linear. However, expansion (1) and thus the linear equation (2) also apply to nonlinear dynamics. This is because the linear system (2) only intends to approximate the particular trajectory in the left hand side of (1), not the complete dynamics of the nonlinear system. In other words, if a different trajectory of the nonlinear system is given, for different initial conditions, matrix A is also generally different. In any event, for generic A, all solutions of (2) are of the form displayed in the right hand side of (1), where the modes un are the eigenvectors of A and δn+iωn are the associated eigenvalues. However, not all expansions of form (1) are solutions of (2) because (2) provides at most J modes, while the number of terms involved in (1), N, which is known as thespectral complexity of expansion (1), can be larger than J. More precisely, defining thespatial complexity of expansion (1), M≤J, as the dimension of the span of the modes u1,…,uN, the linear problem (2) is appropriate to compute (1) only if M=N (note that M cannot be larger than N). We insist here that the spectral and spatial complexities are defined for expansion (1), not for the dynamical system from which the snapshots may come. However, N>M in many dynamical systems of scientific/technological interest, such as the following:(i)
TheLorenz system [9] (and other ODEs) in which M=3, while N may be quite large(ii)
Manynonlinear, infinite-dimensional dynamical systems (PDEs) exhibiting aninertial manifold [10] (limited M) but complex (e.g., quasi-periodic) temporal behavior(iii)
Extrapolating from transients to obtain attractors [11] requiring identifying additional (temporally decaying modes), which increases N(iv)
Experimental noise/numerical errors adding (unphysical) temporal complexity that increases N and must be identified toclean the postprocessed dataWhenN>M, the number of involved modes is increased by using a higher order ODE system. Specifically the linear system (1) can be replaced by (3)ddvdtd=A1dd-1vtdtd-1+A2dd-2vtdtd-2+⋯+Ad-1dvtdt+Advt,for appropriate J×J matrices A1,…,Ad. Obviously, this linear problem can also be written in the form (3), as (4)dv~dt=Av~.Here, v~ is a vector of size Jd that collects v, dv/dt,…, and dd-1v/dtd-1, and the (Jd)×(Jd)-matrix A is a sparse matrix whose blocks are the unit J×J-unit matrix and the matrices A1,…,Ad. The matrix A can be computed via the pseudoinverse applied to v~(t). Now, for an appropriate value of the index d, any dynamics of the form displayed in the right hand side of (1) can be captured. Note that the size of the matrix A is huge when J and d are both large, which would require very large computational resources. However, these resources are greatly decreased by applying this idea after appropriate dimension reduction; see below.In the discrete counterpart of (1), the time variable is discretized at K equispaced values of t, as tk=t1+(k-1)Δt, with Δt=T/K, and a set of snapshots is defined as vk=v(tk). The counterparts of (1), (2), and (3) are (5)vk≃∑n=1Nanuneδn+iωntkfork=1,…,K,(6)vk+1=Rtvkfork=1,…,K-1,(7)vk+d=R1tvk+d-1+R2tvk+d-2+⋯+Rdtvkfork=1,…,K-d,respectively. Obviously, (1) follows from (5) by just setting tk continuous (interpolation in t). Interpolation can also be made in the spatial variable, which yields the continuoustemporal DMD expansion (cf. (1)): (8)vx,t≃∑n=1Nanunxeδn+iωntfort1≤t≤t1+T,where the modes un depend on the spatial variable and the vector v may account for a vector state variable, such as the velocity vector in fluid flows.Equation (6) is the essence ofstandard DMD [6] that, as explained above, only applies when the spatial and spectral complexities coincide. Equation (7), instead, is the basis ofhigher order DMD (HODMD) [8]. This method permits computing the expansion (1) in the general case M<N choosing an appropriate value of the tunable index d. Note that the main difference between (6) and (7) is that the latter equation relates each snapshot with not only the last snapshot but also the former d-1 snapshots. In fact, HODMD reduces to standard DMD if d=1. For d>1, instead, HODMD can be seen as a result of applying standard DMD to a set of enlarged snapshots containing also the delayed snapshots. Thus, HODMD synergically combines the advantages of standard DMD and some consequences [12] of the delayed embedding theorem by Takens [13], who followed and formalized former seminal ideas by Packard et al. [14]. In fact, combination of DMD and delayed snapshots had previously been suggested [7] and performed [15], in a spirit different from that in [8]. Other applications of Takens’ theorem in different contexts (not related to DMD) include, e.g., modal parameter identification and model reduction [16], identification of invariant sets and bifurcations from time-series [17], and analysis of high-dimensional time-series combining the Laplacian eigenmaps with time-lagged embedding [18]. On the other hand, some of the many improvements of standard DMD based on different ideas (not related to time-lagged snapshots) are as follows:(i)
Optimized DMD [19], in which the expansion is computed via an optimization problem. Compared to this method, standard DMD and HODMD are simpler, purely linear algebra driven approaches. However, some linear algebra oriented (and thus less computationally expensive) approaches have been developed to treat separable nonlinear least squares problems [20], in particular, that involved in optimized DMD [21].(ii)
Sparsity promoting DMD [22], which uses a penalty to identify a smaller set of important modes via convex optimization techniques.(iii)
Extended DMD [23], which extends the standard DMD approximation to include more basis functions, allowing the method to capture more complex behavior.(iv)
Total least squares DMD [24], which combines total least squares and standard DMD to treat noisy data. Let us mention here that other improvements of standard DMD have been developed to correct the effect of noise. Among these, Dawson et al. [25] also tried to characterize the noise properties, Takeishi et al. [26] combined a Bayesian formulation with DMD, and Dicle et al. [27] computed low-rank factors of the DMD operator while satisfying the total least squares constraints.On the other hand, for spatially developing flows,t may be replaced in (8) by a distinguished spatial variable x, obtaining thespatial DMD expansion(9)vx,y,t=∑m=1Mamumy,teνm+iκmx,where y collects the remaining transversal spatial variables and νm and κm are the spatial growth rates and wavenumbers, respectively. Still, a simultaneous exponential behavior in space and time may be relevant, considering thespatiotemporal Koopman decomposition (STKD) [28](10)vx,y,t=∑m=1M∑n=1Namnumnyeνm+iκmx+δn+iωnt,which will be obtained by an appropriate STKD method described below. This method roughly consists in applying HODMD along the x and t coordinates simultaneously. Expansion (10) has been attempted in [29] in a data-driven fashion, through an approach based on Koopman operator theory plus the use of appropriate basis functions. See also the analysis in [30], where some aspects of the theory of operator-valued kernels for machine learning with delay-coordinate maps of dynamical systems are combined. Similar spatiotemporal expansions have also been already derived (via a not purely data-driven method) for the Navier-Stokes (NS) equations in [31] by first projecting the NS equations into the temporal modes and then applying spatial DMD to these modes using appropriate basis functions.The modesumn(y) are appropriately normalized and expansion (10) will be further truncated by retaining only those modes whose amplitudes amn are larger than a given threshold. In other words, expansion (10) will not generally contain the M×N modes indicated in (10), but a smaller number of modes, which is emphasized rewriting (10) as (11)vx,y,t=∑m,namnumnyeνm+iκmx+δn+iωnt.It is interesting to note that the methods considered above are related to former seminal ideas by Koopman [32]; see also [5, 31, 33].Koopman operator theory is concerned withobservability in dynamical systems. Specifically, for a dynamical system v(t), we may consider the temporal shift t→t+Δt, with fixed Δt. For an arbitrary vector observable function g, the Koopman operator Rt is defined as (12)Rt·gvt=gvt+Δt.This operator is linear, because obviously Rt·(α1g1+α2g2)=α1Rt·g1+α2Rt·g2, and infinite-dimensional, since it acts on the infinite-dimensional space of observables. Thus, its eigenvalues and eigenvectors allow for giving an exact representation of any observable. Setting g=the identity, the Koopman operator (12) reduces to (6) and is thus related to the standard DMD method. Similarly, considering an observable that also takes into account delayed snapshots, the Koopman operator reduces to (7) and is thus related to HODMD. Still, particular realizations of the spatiotemporal Koopman operator Rx,t, defined as (13)Rx,t·gvx,y,t)=vx+Δx,y,t+Δt,are related to the STKD method mentioned above. This connection is why we call this method spatiotemporal Koopman decomposition.Finally, we note that in the weakly nonlinear analysis of oscillatory pattern forming systems [34, 35], some equations similar to (9) and (10) are also obtained. However, these methods require both knowing the governing equations and assuming that the dynamics are small perturbations of a known steady or periodic solution. In contrast, DMD, HODMD, and STKD provide approximations of type (9) and (10) in a purely data-driven fashion for arbitrary fully nonlinear dynamics, which are great advantages, specially in the case of the analysis of experimental data. Note in this context that expansions (9) and (10) can be seen assemianalytical expressions that allow for the fast online computation of the associated spatiotemporal data, which leads to apurely data-driven reduced order model (ROM). Among the most important applications of such a ROM, (i) it avoids numerical computation, which may be quite computationally expensive in large scale systems, or even not possible when the model is obtained from experimental data; (ii) it predicts different flow states through available data reducing the number of data collected in experiments; and (iii) it provides efficient tools for various applications, including flow control [36] and optimization design [37].With the above in mind, the main object of this paper is to present the HODMD and STKD methods, illustrating them in simple toy models and applying them to various problems of scientific and technological interest. To this end, the HODMD method is described and illustrated in a simple toy model in Section2, while several applications of this method are considered in Section 3. The STKD method is developed in Section 4 and applied to various problems in Section 5. The paper ends with some concluding remarks in Section 6.


## 2. Higher Order Dynamic Mode Decomposition

As above, for a time-dependent dynamical system given by the vector state variablev(t), of size J, we consider a set of K snapshots at equispaced values of t, (14)vk≡vtk,withtk=t1+k-1Δtfork=1,…,K.If the dynamical system is infinite-dimensional, the finite-dimensional state vector v generally results from spatial discretization. The snapshots are organized in a snapshot matrix (whose columns are the snapshots), as (15)V1K=v1,v2,⋯,vK.In the sequel, the accuracy of the various approximations of the snapshots will be measured in terms of therelative root mean square (RMS) error, defined as (16)RRMSerror=∑k=1Kvkapprox-vk22∑k=1Kvk22,where ·2 denotes the usual Euclidean norm. The right hand side of (16) measures the collective RMS error in space and time normalized with the spatiotemporal RMS norm of the snapshots. Obviously, other norms (omitted here for the sake of brevity) would be more appropriate when the dynamics exhibit disparate localized behavior in space and/or time.For simplicity, only the temporal expansion (5) will be considered in this section. The spatial expansion (9) is obtained by just interchanging the roles of x and t.


### 2.1. The HODMD Method

The HODMD algorithm proceeds in two steps (see [8] for more details).Step 1 (dimension reduction).
Applying truncatedsingular value decomposition (SVD) [38] to the snapshot matrix V1K yields (17)V1K≃PΣQ⊤≡PT^1Korvk=Pt^k,withT^1K=ΣQ⊤,where the matrix T^1K is known as the dimension-reduced snapshot matrix and their columns, thedimension-reduced snapshots, t^k. The number of retained modes, which is the spatial complexity defined above, M, is determined by a (tunable) tolerance ε1, requiring that the RRMS error of the approximation, as defined in (16), be smaller than ε1; this error is easily computed in terms of the singular values using well-known SVD formulae [38]. For noisy data we may choose ε1 comparable to the noise level, which helps to clean the data [39]. Now we note that if we had a DMD expansion for the dimension-reduced snapshots, as (18)t^k≃∑n=1Nanu^neδn+iωntkfork=1,…,K,which in fact will be computed in the next step, substituting this expansion into (17) (namely, premultiplying (18) by the matrix P), then we obtain the DMD expansion (5) we are looking for.Step 2 (the DMD-d algorithm for the dimension-reduced snapshots).
Now, we must treat the standard DMD and the HODMD algorithms separately. For ordinary DMD, we consider the Koopman assumption (6) for the dimension-reduced snapshots, which reads (19)t^k+1≃R^tt^k,
where theM×M matrix R^t (reducedKoopman matrix) is computed from the dimension-reduced snapshots via the pseudoinverse. Now, the reduced expansion (18) is computed as follows. The reduced modes u^n are normalized eigenvectors of R^t, the exponents are related to the eigenvalues, μn, as δn+iωn=log⁡μn/Δt, and the amplitudes an are computed via least squares fitting (as in optimized DMD [19]) of the dimension-reduced snapshots in expansion (18). This completes the algorithm for standard DMD.
Concerning HODMD, the counterpart of (6) for the reduced snapshots, (20)t^k+d≃R^1t^k+d-1+R^2t^k+d-2+…+R^dt^k,is treated using theenlarged snapshots(21)t~k≡t^kt^k+1…t^k+d-2t^k+d-1.Applying Step 1 above (a second dimension reduction that eliminates redundancies between the snapshots and delayed snapshots) gives a new set of enlarged dimension-reduced snapshots, t~^k, which are assumed to satisfy the counterpart of (19), namely, t~^k+1≃R~tt~^k, where the matrix R~t is computed via the pseudo inverse, as above. The eigenvalues of eigenvectors of R~t permit obtaining the counterpart of (18) for the enlarged dimension-reduced snapshots. Considering only the first M components in this expansion, we obtain the dimension-reduced expansion (18), where the mode amplitudes an are computed via least square fitting with the dimension-reduced snapshots, which completes the derivation of (18) in the HODMD method.
Once expansion (18) has been computed either for standard DMD or for HODMD, the expansion is sorted in decreasing order of the mode amplitudes and is further truncated by eliminating those modes such that (22)ana1<ε2,
for some tunable parameterε2.Some remarks about the method are now in order:(i)
A solver for the DMD-d algorithm described above can be found in [40].(ii)
It is interesting to note that the enlarged snapshots, containing also the delayed snapshots, are defined for the reduced snapshots (whose size isM), not for the (usually much larger) original snapshots, which, as anticipated, would greatly increase the computational cost of the method. In fact, the ratio of the computational cost of HODMD to that of standard DMD scales with a first or second power (depending on the comparative values of J, Md, and K) of J/(Md), which is of order one in most applications (but obviously not always).(iii)
Considering periodic or quasi-periodic dynamics, and measuring the accuracy of the DMD reconstruction with the RRMS error defined in (16), the method is quite robust in connection with the tunable parameters d, ε1, and ε2. In particular, the plot of the RRMS error versus d is fairly flat near the minimum, which means that the selection of d is not critical. In principle, if the data were exact, the smaller ε1 and ε2, the better the approximation. However, in the presence of noise or discretization/truncation errors, ε1 and ε2 must be not too small to avoid that the method captures unphysical modes. In other words, these two parameters must be comparable to the noise/errors level, but see the next remark. Finally, the timespan tK-t1 and the time shift between snapshots, Δt, must be selected as somewhat larger and smaller than the largest and smallest involved periods, respectively.(iv)
An obvious trade-off of the HODMD method is that we have an additional tunable parameterd, which is no present in standard DMD and must be chosen via trial-and-error. However, as explained above, the method is robust in connection with the parameter d.(v)
The parametersε1 and ε2 can be seen as governing spatial and temporal truncation, respectively. For finite accuracy and given values of the tunable parameters, a new application of the method to the HODMD reconstructed solution does not necessarily give the same reconstruction. This is specially interesting for noisy data, in which the iterative application of the method (iterative HODMD, see [41]) improves noise filtering (beyond the well-known [39] noise filtering effect associated with the application of truncated SVD in Step 1). It is to be noted that iterating HODMD with fixed tolerances ε1 and ε2 eliminates both spurious SVD modes and spurious DMD frequencies, decreasing along the iterations both the spatial and spectral complexities, M and N, respectively. The iteration is terminated when M and N do not decrease further. By choosing ε1 and ε2 comparable to the noise level, the method permits uncovering physical phenomena that were masked by noise in the original data. A striking example of the ability of the iterative HODMD method to filter noise and uncovering the physically relevant dynamics will be given in Section 3.2.(vi)
When more than one spatial coordinate is present, and spatial discretization is made in a structured mesh in some of the spatial coordinates, the given data can be organized in asnapshot tensor instead of a snapshot matrix. In this case, standard DMD can be substituted in Step 1 (see [41]) by HOSVD [42], a well-known extension of standard DMD by Tucker [43]. HOSVD decomposes the tensor as a linear combination of tensor products of modes in the various directions. The DMD-d algorithm is applied separately (perhaps with different values of d) to the various sets of modes.


### 2.2. Illustration of the HODMD Method in a Toy Model

Let us first illustrate the method in the following quasi-periodic simple toy model(23)vt=1+12sin⁡tsin⁡3t.This model exhibits the incommensurable fundamental frequencies 3±1, which are somewhat disparate from each other and their joint harmonics. The spatial complexity is M=1, while the spectral complexity N (which increases as the accuracy of the reconstruction is required to be more and more strict) is larger than one. The convenience of using delayed snapshots when the spatial complexity is 1 was already mentioned in [7]. The HODMD algorithm DMD-800 is applied to a set of K=2000 snapshots in the timespan 0≤t≤32, with tolerances ε1=10-10 and ε2=10-4. The method identifies N=27 frequencies and reconstructs the snapshots with RRMS error, as defined in (16), ~2·10-4 in this interval. Thus, the original quasi-periodic function and its reconstruction are plot-indistinguishable (see Figure 2-left). Moreover, the extrapolation of this reconstruction to the interval 180≤t≤200 (Figure 2-right) is still quite good, since the RRMS error is ~2.7·10-3. On the contrary, the standard DMD algorithm DMD-1 (Figure 2-left, thin dashed line) identifies just one spurious frequency and gives a O(1) RRMS error.Figure 2
Toy model function defined in (16) (thin solid lines) and the reconstruction using DMD-800 (thick dashed lines) and DMD-1 (thin dashed line).More realistic applications of the HODMD method and its iterative extension are considered in the next section.


## 2.1. The HODMD Method

The HODMD algorithm proceeds in two steps (see [8] for more details).Step 1 (dimension reduction).
Applying truncatedsingular value decomposition (SVD) [38] to the snapshot matrix V1K yields (17)V1K≃PΣQ⊤≡PT^1Korvk=Pt^k,withT^1K=ΣQ⊤,where the matrix T^1K is known as the dimension-reduced snapshot matrix and their columns, thedimension-reduced snapshots, t^k. The number of retained modes, which is the spatial complexity defined above, M, is determined by a (tunable) tolerance ε1, requiring that the RRMS error of the approximation, as defined in (16), be smaller than ε1; this error is easily computed in terms of the singular values using well-known SVD formulae [38]. For noisy data we may choose ε1 comparable to the noise level, which helps to clean the data [39]. Now we note that if we had a DMD expansion for the dimension-reduced snapshots, as (18)t^k≃∑n=1Nanu^neδn+iωntkfork=1,…,K,which in fact will be computed in the next step, substituting this expansion into (17) (namely, premultiplying (18) by the matrix P), then we obtain the DMD expansion (5) we are looking for.Step 2 (the DMD-d algorithm for the dimension-reduced snapshots).
Now, we must treat the standard DMD and the HODMD algorithms separately. For ordinary DMD, we consider the Koopman assumption (6) for the dimension-reduced snapshots, which reads (19)t^k+1≃R^tt^k,
where theM×M matrix R^t (reducedKoopman matrix) is computed from the dimension-reduced snapshots via the pseudoinverse. Now, the reduced expansion (18) is computed as follows. The reduced modes u^n are normalized eigenvectors of R^t, the exponents are related to the eigenvalues, μn, as δn+iωn=log⁡μn/Δt, and the amplitudes an are computed via least squares fitting (as in optimized DMD [19]) of the dimension-reduced snapshots in expansion (18). This completes the algorithm for standard DMD.
Concerning HODMD, the counterpart of (6) for the reduced snapshots, (20)t^k+d≃R^1t^k+d-1+R^2t^k+d-2+…+R^dt^k,is treated using theenlarged snapshots(21)t~k≡t^kt^k+1…t^k+d-2t^k+d-1.Applying Step 1 above (a second dimension reduction that eliminates redundancies between the snapshots and delayed snapshots) gives a new set of enlarged dimension-reduced snapshots, t~^k, which are assumed to satisfy the counterpart of (19), namely, t~^k+1≃R~tt~^k, where the matrix R~t is computed via the pseudo inverse, as above. The eigenvalues of eigenvectors of R~t permit obtaining the counterpart of (18) for the enlarged dimension-reduced snapshots. Considering only the first M components in this expansion, we obtain the dimension-reduced expansion (18), where the mode amplitudes an are computed via least square fitting with the dimension-reduced snapshots, which completes the derivation of (18) in the HODMD method.
Once expansion (18) has been computed either for standard DMD or for HODMD, the expansion is sorted in decreasing order of the mode amplitudes and is further truncated by eliminating those modes such that (22)ana1<ε2,
for some tunable parameterε2.Some remarks about the method are now in order:(i)
A solver for the DMD-d algorithm described above can be found in [40].(ii)
It is interesting to note that the enlarged snapshots, containing also the delayed snapshots, are defined for the reduced snapshots (whose size isM), not for the (usually much larger) original snapshots, which, as anticipated, would greatly increase the computational cost of the method. In fact, the ratio of the computational cost of HODMD to that of standard DMD scales with a first or second power (depending on the comparative values of J, Md, and K) of J/(Md), which is of order one in most applications (but obviously not always).(iii)
Considering periodic or quasi-periodic dynamics, and measuring the accuracy of the DMD reconstruction with the RRMS error defined in (16), the method is quite robust in connection with the tunable parameters d, ε1, and ε2. In particular, the plot of the RRMS error versus d is fairly flat near the minimum, which means that the selection of d is not critical. In principle, if the data were exact, the smaller ε1 and ε2, the better the approximation. However, in the presence of noise or discretization/truncation errors, ε1 and ε2 must be not too small to avoid that the method captures unphysical modes. In other words, these two parameters must be comparable to the noise/errors level, but see the next remark. Finally, the timespan tK-t1 and the time shift between snapshots, Δt, must be selected as somewhat larger and smaller than the largest and smallest involved periods, respectively.(iv)
An obvious trade-off of the HODMD method is that we have an additional tunable parameterd, which is no present in standard DMD and must be chosen via trial-and-error. However, as explained above, the method is robust in connection with the parameter d.(v)
The parametersε1 and ε2 can be seen as governing spatial and temporal truncation, respectively. For finite accuracy and given values of the tunable parameters, a new application of the method to the HODMD reconstructed solution does not necessarily give the same reconstruction. This is specially interesting for noisy data, in which the iterative application of the method (iterative HODMD, see [41]) improves noise filtering (beyond the well-known [39] noise filtering effect associated with the application of truncated SVD in Step 1). It is to be noted that iterating HODMD with fixed tolerances ε1 and ε2 eliminates both spurious SVD modes and spurious DMD frequencies, decreasing along the iterations both the spatial and spectral complexities, M and N, respectively. The iteration is terminated when M and N do not decrease further. By choosing ε1 and ε2 comparable to the noise level, the method permits uncovering physical phenomena that were masked by noise in the original data. A striking example of the ability of the iterative HODMD method to filter noise and uncovering the physically relevant dynamics will be given in Section 3.2.(vi)
When more than one spatial coordinate is present, and spatial discretization is made in a structured mesh in some of the spatial coordinates, the given data can be organized in asnapshot tensor instead of a snapshot matrix. In this case, standard DMD can be substituted in Step 1 (see [41]) by HOSVD [42], a well-known extension of standard DMD by Tucker [43]. HOSVD decomposes the tensor as a linear combination of tensor products of modes in the various directions. The DMD-d algorithm is applied separately (perhaps with different values of d) to the various sets of modes.


## 2.2. Illustration of the HODMD Method in a Toy Model

Let us first illustrate the method in the following quasi-periodic simple toy model(23)vt=1+12sin⁡tsin⁡3t.This model exhibits the incommensurable fundamental frequencies 3±1, which are somewhat disparate from each other and their joint harmonics. The spatial complexity is M=1, while the spectral complexity N (which increases as the accuracy of the reconstruction is required to be more and more strict) is larger than one. The convenience of using delayed snapshots when the spatial complexity is 1 was already mentioned in [7]. The HODMD algorithm DMD-800 is applied to a set of K=2000 snapshots in the timespan 0≤t≤32, with tolerances ε1=10-10 and ε2=10-4. The method identifies N=27 frequencies and reconstructs the snapshots with RRMS error, as defined in (16), ~2·10-4 in this interval. Thus, the original quasi-periodic function and its reconstruction are plot-indistinguishable (see Figure 2-left). Moreover, the extrapolation of this reconstruction to the interval 180≤t≤200 (Figure 2-right) is still quite good, since the RRMS error is ~2.7·10-3. On the contrary, the standard DMD algorithm DMD-1 (Figure 2-left, thin dashed line) identifies just one spurious frequency and gives a O(1) RRMS error.Figure 2
Toy model function defined in (16) (thin solid lines) and the reconstruction using DMD-800 (thick dashed lines) and DMD-1 (thin dashed line).More realistic applications of the HODMD method and its iterative extension are considered in the next section.


## 3. Applications of HODMD

The discussion about spectral and spatial complexities is a key point to understand the wide range of applications of HODMD. As anticipated in Section1 and further explained in [11], there are several natural cases in which the spectral complexity is larger than the spatial complexity. Here, we shall concentrate on two such cases, one driven by numerical data from the complex Ginzburg-Landau equation (CGLE) and another by experimental data from the zero-net-mass-flux (ZNMF) jet, which will be considered in Sections 3.1 and 3.2, respectively. Additional applications will be briefly addressed in Section 3.3.


### 3.1. Application to the CGLE

Let us consider the following CGLE with periodic boundary conditions:(24)∂tA=1+iα∂xx2A+μu-1+iβA2A,(25)A-1,t=A1,t,∂xA-1,t=∂xA1,t.The dependent variable A is complex and can be seen as thecomplex amplitude when using (24) as the weakly nonlinear “normal form”, which applies at the onset of oscillatory instabilities in many relevant physical systems [34, 35]. The CGLE is a well-known paradigm of pattern forming systems that exhibits complex dynamics [44] due to the modulational instability if αβ<-1 (Newell’s condition [44]) and μ exceeds a threshold value. For typical values of α, β, and μ, the CGLE possesses an inertial manifold contained in a linear subspace of dimension 20-30 (upper limit of the spatial complexity), while the dynamics may be complex, meaning that the spectral complexity can be large. Thus, the standard DMD-1 algorithm fails except for the simplest monochromatic attractors [8].Problem (24)-(25) is invariant under the O(2)×SO(2) symmetry group generated by spatial translations and reflections, and phase shifts in A, (26)x→x+c1,x→-x,A→Aeic2.Because of the latter symmetry, it is convenient to write down expansion (1) for the CGLE as (27)Ax,t=eiω0t∑n=1Nanunxeiωnt,where the frequencies ω1,ω2,… are linear combinations with integer coefficients of some fundamental frequencies, which are generally incommensurable with ω0. Thus, the dynamics of A are more complex than those of A. In other words, A is periodic if A is stationary, and A represents a 2-torus if A is periodic. In the sequel, we consider just two attractors, one such that |A| is periodic and another such that A is quasi-periodic. Both attractors are obtained for (28)α=-10,β=10,and two values of μ. Note that since α and β are both somewhat large, dispersion and nonlinear detuning are large compared to diffusion and nonlinear damping, respectively, which is demanding from the computational point of view and promotes complex dynamics. For the construction of the databases, the equation is solved numerically using a Crank-Nicolson plus Adams-Bashforth scheme [45] for temporal discretization, with time step Δt=10-5, while the second-order spatial derivative is discretized using centered finite differences in a uniform grid of 1000 points. HODMD is applied to a set of 2000 snapshots, collected in an appropriate time interval of length 1, disregarding a transient time-stage. In the two cases considered below, a good value of the index d when applying the HODMD algorithm is d=50, while the remaining tunable parameters are ε1=10-8 and two values of ε2, namely, (29)ε2=5·10-3highaccuracy,ε2=10-4lowaccuracy.Selecting these values of the tunable parameters as common for the two cases considered below illustrates the robustness of the method, which is further checked noting that low and high accuracy computations are consistent among each other, and the results coincide when different timespans are used in the computations.The first case to be considered in that forμ=20, in which using the DMD-50 algorithm, the RRMS reconstruction error (see (16)) is ~6·10-3 and ~2·10-5 for low and high accuracy, respectively, retaining 6 and 12 DMD modes, respectively. The semilogarithmic plot of the DMD amplitude versus frequency is plotted in Figure 3-left, where it can be seen that the relevant points are approximately contained in two straight lines centered at ω0=11.6, which means that expansion (27) converges spectrally; using the standard DMD method (DMD-1 algorithm), instead, gives spurious results, plotted with red symbols in the figure. The remaining points obtained via the DMD-50 algorithm lie equispaced in frequency, with a distance between them equal to ω=118.7, which means that the frequencies ωn appearing in expansion (27) are all (positive or negative) harmonics of ω. Thus, A is time-periodic (see Figure 3, middle), with a period 2π/ω, while A is quasi-periodic (2-torus), with fundamental frequencies ω0 and ω. Since ω0 is small compared to ω the pattern is slowly modulated (Figure 3, right). The spatial modes appearing in the expansion (27), un, which are not plotted for the sake of brevity, exhibit a reflection symmetry around x=-0.17 (indicated with a vertical white line in Figure 3), meaning that the whole pattern also exhibits this symmetry, consistently with the invariance of the problem under the actions (26). This implies that the pattern can be seen as a standing wave (SW), as indicated with the vertical while solid line in Figure 3, but it can also be considered as either a right or left traveling wave (TW), as indicated with the oblique dashed lines in this figure. This dual interpretation cannot be ascertained from the shape of the modes un, but this issue will be revisited and further checked in Section 5.1.Figure 3
The CGLE forμ=20, plotting the amplitudes versus frequencies in the DMD expansion (left plot) for low (black circles) and high (black crosses) accuracy computations using DMD-50, and using DMD-1 (red symbols), and the color map of R(A) (middle plot) and |A| (right plot).As a second case, we consider the valueμ=22 (Figure 4), in which the dynamics are more complex. Again, the DMD-1 algorithm yields spurious results. Instead, for low and high accuracy (see (29)), the DMD-50 algorithm retains 17 and 49 modes, giving a RRMS reconstruction error (see (16)) ~8·10-3 and ~2·10-4, respectively. The main difference with the case considered in Figure 3 is that now the plot of the amplitude versus the frequency forms a triangular pattern that also decays spectrally but shows that, in addition to ω0=-21.2, the frequencies in the expansion in the right hand side of (27) are linear combinations with integer coefficients of the fundamental frequencies ω1=156.6 and ω2=22.2, which are not easily identified in Figure 4, left. Thus, guessing that the three fundamental frequencies are incommensurable, the attractor for |A| is contained in a 2-torus while that for A lies in a 3-torus. Finally, the middle and right plots in Figure 4 show that the pattern is a TW but, as in the former case, elucidating this issue from the shape of the modes un(x) appearing in (27) turns out to be a difficult (and unnatural) issue. The STKD method, which will be applied to the CGLE in Section 5.1, instead, will solve this issue in a natural way.Figure 4
Counterpart of Figure3 for =22.In summary, the HODMD method is able to precisely compute the temporal DMD expansion for the CGLE, giving robust results and distinguishing between periodic and quasi-periodic phenomena. However, identifying SWs and TWs requires additional ingredients, which will be included in the STKD method, developed in Section4 and applied to the CGLE in Section 5.1.


### 3.2. Application to the Zero-Net-Mass-Flux Jet

A zero-net-mass-flux (ZNMF) jet [46, 47] is a pulsatile jet promoted by a piston or membrane that interchanges fluid flow through an orifice with an outer ambient fluid that is at rest far from the orifice. The net mass flux is zero, but at large Reynolds number there is a nonzero (and large) net flux of momentum that, in addition to its own scientific interest, is useful in various fields. For instance, synthetic jet actuators promote pulsatile jets [48] that are useful in, e.g., mixing enhancement [49], jet vectoring [50], heat transfer [51], and active flow control of boundary layer separation [52]. These jets are also relevant in nature, in connection with squid, octopus, salp, and jellyfish swimming [53].It must be noted that the nonzero momentum flux at large Reynolds number is due to the fact that, even if the periodic forcing is temporally reflection symmetric, the flow is highly nonreflection symmetric, as illustrated in Figure5, where an approximately rotationally symmetric jet produced by a circular orifice is considered. As can be seen, a vortex ring is formed in the injection phase that travels downstream, while the suction phase is characterized by a saddle point that separates the fluid reentering the cavity from that still moving downstream.Figure 5
Pseudo-streamlines in a meridional section for two representative snapshots of the injection (left) and suction (right) phases in the ZNMF jet.The ZNMF jet promotes a complex flow structure that involves various, disparate spatiotemporal scales, whose understanding remains as an open topic. Thus, this is a convenient problem to test the performance of the HODMD method. Also, in order to test the noise filtering ability of the iterative HODMD method, the analysis will be based on experimental data resulting from particle image velocimetry (PIV), as obtained [46, 47] in the experimental facility for ZNMF jets of the Laboratory for Turbulence Research and Combustion at Monash University. Such measurements exhibit an error ~2.4% at the 95% confidence level. Two nondimensional parameters characterize the flow conditions, namely theReynolds and Strouhal numbers, defined as Re=UD/ν and St=fD/U, respectively, where D is the jet orifice diameter, ν is the kinematic viscosity, U is related to the peak oscillation velocity of the piston V^p as U=5V^p/2 (see [41] for details), and f is the piston oscillation frequency. In the present case, we have Re=13329 and St=0.03, which give a flow lying in the laminar-transitional regime. The flow is fully turbulent in the far field, but laminar in the near field (near the orifice), which will be the region of the jet considered here; see [41] for further details on the experimental facility and the PIV measurement.The orifice is circular and the resulting near field flow is approximately axi-symmetric, at least in the large scale. The experimental data consists inK=1872 snapshots representing 12 piston oscillation cycles, giving the streamwise and radial velocity components, u(xi,rj,tk) and v(xi,rj,tk), respectively, at I×J=52×68 grid points in a meridional rectangle (interrogation window), where x and r are the streamwise and radial coordinates, respectively, with origin at the center of the orifice and t is the time variable, with origin at the beginning of the suction phase. The time variable has been rescaled such that the (nondimensional) oscillating frequencies are measured in terms of the associated Strouhal number defined above, meaning that the piston oscillating period is T=2π/0.03. Thus, we have a snapshot tensor (instead of a matrix) and HOSVD must be used instead of standard SVD in the first step of the HODMD method (see [41] for details). Assuming permanent dynamics (zero or very small growth rates) the aim of the iterative HODMD method is to obtain an expansion of the form (30)vx,r,t≃∑n=1Nunx,reiωnt,The iterative DMD-d algorithm has been applied to obtain this expansion using two values of the index d, namely, d=1 (standard DMD) and d=700 (strictly higher order DMD), with two tolerance sets, namely, (31)ε1=ε2=2.4·10-2lowaccuracy,ε1=ε2=10-2highaccuracy,which correspond to the estimated experimental error and more precise computations, respectively. Using the latter tolerances, below the experimental error, will permit elucidating the ability of the iterative HODMD method to uncover physical phenomena that are masked by the experimental errors in the given data.Figure6 shows the amplitude versus frequency (nondimensionalized as the Strouhal number) obtained in the analyses. DMD-1 captures 11 modes, all of them spurious.Figure 6
Amplitude versus St obtained for low accuracy via DMD-1 (red crosses) and DMD-700 (blue circles) and for high accuracy using DMD-700 (black circles).On the contrary, DMD-700 capturesN=19 and N=49 DMD modes using low and high accuracy, respectively. It is remarkable that the spatial complexity M is much smaller than the spectral complexity, N, since M=5 and M=11 for low and high accuracy, respectively, which justifies the failure of DMD-1. The method shows that the flow is periodic in the near field, with dominant frequency St=0.03 (piston oscillation frequency). Since the data are real, the DMD modes appear in complex-conjugate pairs, meaning that the retained modes contain the mean flow (St=0) and 9 and 24 nonzero harmonics for low and high accuracy, respectively.The suction phase (see Figure5-left) is not quite demanding, since it just shows a saddle point at the axis of the jet. The reconstructions using DMD-700 with low and high accuracies are qualitatively similar to their counterparts for the original PIV data. The main difference is that the iterative HODMD method gives smoother results with high accuracy than with low accuracy, and both are smoother than the original PIV data. This is an indication of the ability of the method to filter experimental errors, but this issue cannot be confirmed because the actual physical data are not available. Concentrating on the reconstructions with low accuracy, it is seen that the saddle point is created near the orifice at the beginning of this phase and travels downstream to a distance that is approximately equal to 2D, where it remains steady for most part of the suction phase (Figure 7, left and middle plots). This saddle point is blown downstream at the beginning of the injection phase, where a new saddle point is created near the orifice that travels dowstream very fast (Figure 7, right plot) and leaves the near field, in a short timespan of the order of T/100. This saddle point must be present for topological reasons, but it had never been seen in studies previous to [41], just because it is quite short living. In fact, it is not present in any of the snapshots used to apply the HODMD method. It must be noted that this saddle point would be present (with some spatial noise) in the experimental snapshots if the sampling frequency was sufficiently large. Thus, the identification of this saddle point illustrates well the advantages of the temporal interpolation that is implicit in the HODMD method.Figure 7
Two snapshots in the suction phase, fort=T/18 (left plot) and t=2T/9 (middle plot), and a snapshot at the very beginning of the injection phase, for t=T/2+T/72 (right plot).The injection phase (see Figure5, right) is more interesting in the present context. Figure 8 shows a comparison of the original data with the reconstructions provided by DMD-700 for low and high accuracy. As can be seen, the reconstructions are much smoother than the original PIV data, as it happened in the suction phase. However, the reconstruction with high accuracy (recall that this is below the experimental error) is qualitatively different, since it shows some small vortices leaving the edge of the cavity that cannot be even guessed in the original PIV data and do not appear in the low accuracy reconstruction because they involve a quite small energy. These small vortices are smooth, which may indicate that they are not artifacts associated with higher order modes resulting from experimental errors. However, this cannot be ascertained since the actual physical dynamics is not available. Therefore, we have conducted some numerical simulations to elucidate this point. The value of the Reynolds number in the experimental data considered above, Re~13000, is a quite large value that would need huge numerical resources, even using the computationally efficient solver Nek5000 [54] (an open source spectral element code). Therefore, we have performed computations using the lower Reynolds number Re=1000 and the same Strouhal number as in the experimental data considered above, namely, St=0.03. The numerical counterpart of the snapshot considered in Figure 8 is given in Figure 9, which confirms the presence of the small vortices.Figure 8
Instantaneous pseudo-streamlines and streamwise velocity contours for a representative snapshot (t=T/9) in the injection phase, showing the original PIV data (left plot) and the reconstructions using DMD-700 for low (middle plot) and high (right plot) accuracy.Figure 9
Counterpart of the snapshot considered in Figure5 but as numerically computed for Re=1000.In summary, the results above illustrate well the ability of the iterative HODMD method to both (i) take advantage of temporal interpolation to uncover fast events (such as the saddle point at the beginning of the injection phase) that are not present in the snapshots and (ii) filter experimental errors so well that some small patterns (such as the small vortices in the injection phase) can be identified that were completely masked by the experimental noise. In addition, the analysis above permits constructing a purely data-driven (from experimental data) ROM that reproduces well the periodic behavior of the ZNMF jet and shows a fast online operation.


### 3.3. Additional Applications

Some additional applications are now briefly quoted.The HODMD method permits identifying flow instabilities by monitoring those DMD modes whose growth rate changes sign. On the other hand, as anticipated, retaining only those terms appearing in the DMD expansion (1) that exhibit zero (or very small) growth rate, the final attractor may be obtained from transient data. This latter idea can be used to strongly decrease the CPU time required to compute the final attractor, avoiding continuing the numerical integration until final attractor is reached, which may require an extremely large timespan, specially near bifurcation points. Transient dynamics are usually found in the initial stage of a numerical simulation, where a large number of modes (either physical or spurious) appear simultaneously.Let us illustrate these two applications by considering the cylinder wake, namely, the wake behind a circular cylinder immersed in a free stream with constant, spatially uniform velocity, which is a basic fluid dynamics problem [55, 56]. This problem depends on a unique nondimensional parameter, the Reynolds number Re=UD/ν, where U is the free stream velocity, D is the cylinder diameter, and ν is the dynamic viscosity of the fluid, assumed to be incompressible. Thus, the nondimensional velocity and pressure obey the continuity and Navier-Stokes equations: (32)∇·v=0,(33)∂v∂t+v·∇v=-∇p+1ReΔv,which are usually integrated imposing v is equal to the free stream velocity at the inlet boundary, no-slip boundary conditions at the boundary of the body, appropriate (nonreflecting) boundary conditions at the sides and exit of the computational domain, and periodic boundary conditions in the spanwise direction, with a period Lz. The wake is two-dimensional and steady for small Re, but at Re≃46 it exhibits a first unsteady instability (Hopf bifurcation [57, 58]) that produces a still two-dimensional but unsteady, periodic von Karman vortex street flow. The flow remains orbitally stable up to Re≃190, where it suffers a secondary bifurcation (Floquet multiplier=1) and becomes three-dimensional [59] for spanwise period Lz>3D (with D the cylinder diameter). The resulting mode is a long-wave, periodic, flow known asmode A. For larger values of Re, at Re=259, a second periodic short-wave mode appears if Lz>1.5D that is known asmode B. Modes A and B are illustrated in Figure 10. These modes are usually obtained via Floquet stability analysis of the basic two-dimensional periodic solution [59–61], but they can also be obtained via HODMD applied to the whole dynamics obtained via numerical simulations [62]. It is relevant to note that for higher values of the Reynolds number, at Re ≃380, a new instability takes place that is associated with a pair of complex-conjugate Floquet multipliers that give a family ofquasi-periodic solutions [61]. These modes were also identified in [62].Figure 10
Streamwise velocity contours for modes A (left) and B (right) obtained at Re=280, with and Lz=6.99.As for temporal extrapolation via the HODMD method, this idea has been used [11] precisely in the three-dimensional cylinder wake, with Re=220 and Lz=4. The numerical data was obtained using the above-mentioned numerical solver Nek5000 [54]. Temporal HODMD extrapolation divided by five the CPU time required computing the final attractor with a reasonable accuracy; see [11] for further details.Additional applications involving numerically generated data include the analysis of the turbulent wake of a cross flow wind turbine [63] and an off-shore wind turbine [64], as well as the analysis of wakes interaction [65].Concerning analyses based on experimental data, a good example, linked to the field of renewable energies, is the postprocess of light detection and ranging (LiDAR) [66, 67] measurements. LiDAR is an experimental method for the remote measurement of the line-of-sight component of wind speed, usually employed in the wind energy industry with different goals, including the maximization of the power generated by the wind turbine. Measurements are made at a few (say, M) points upstream the wind turbine, based on detection of the Doppler shift for light backscattered from natural aerosols transported by the wind in the atmosphere. LiDAR measurements offer M time-dependent signals that present high noise levels (15-20%). Thus, the spatial complexity M is small, making HODMD a suitable tool to analyze the data and detect the dominant frequencies. A great difficulty is that LiDAR experiments only measure at distances between 20 m and 200 m from the measurement device, which is a restriction in certain situations. Figure 11 shows six planes of LiDAR measurements located at distances between 33 m and 201 m upstream the wind turbine. HODMD is used to analyze such data and to predict the wind velocity at a seventh distance (228 m), which is well beyond the LiDAR device measurement range. The idea used in [68] is to extrapolate the spatial modes in the DMD expansion, guessing the data at the seventh distance with relative errors ~ 2%. The computational cost for predicting these measurements is negligible, making it possible to easily update the model in real-time (leading to a fast, purely data-driven ROM).Figure 11
Measurements upstream the wind turbine in the LiDAR experimental campaign. Data are available in the six planes that are closer to the wind turbine, while the data on the seventh plane are to be predicted.Additional applications using experimental data where the HODMD method have proved to give good results including nonisothermal flows around square objects in a channel [69] and the analysis of aeroelastic data in aircraft flight tests [70]. The aim in the latter application is to obtain the natural aeroelastic frequencies, and the associated damping rates, of the various aircraft parts (e.g., wings or fuselage) from data obtained at some sensors along the aircraft when this is subject to various types of impulsive excitation.


## 3.1. Application to the CGLE

Let us consider the following CGLE with periodic boundary conditions:(24)∂tA=1+iα∂xx2A+μu-1+iβA2A,(25)A-1,t=A1,t,∂xA-1,t=∂xA1,t.The dependent variable A is complex and can be seen as thecomplex amplitude when using (24) as the weakly nonlinear “normal form”, which applies at the onset of oscillatory instabilities in many relevant physical systems [34, 35]. The CGLE is a well-known paradigm of pattern forming systems that exhibits complex dynamics [44] due to the modulational instability if αβ<-1 (Newell’s condition [44]) and μ exceeds a threshold value. For typical values of α, β, and μ, the CGLE possesses an inertial manifold contained in a linear subspace of dimension 20-30 (upper limit of the spatial complexity), while the dynamics may be complex, meaning that the spectral complexity can be large. Thus, the standard DMD-1 algorithm fails except for the simplest monochromatic attractors [8].Problem (24)-(25) is invariant under the O(2)×SO(2) symmetry group generated by spatial translations and reflections, and phase shifts in A, (26)x→x+c1,x→-x,A→Aeic2.Because of the latter symmetry, it is convenient to write down expansion (1) for the CGLE as (27)Ax,t=eiω0t∑n=1Nanunxeiωnt,where the frequencies ω1,ω2,… are linear combinations with integer coefficients of some fundamental frequencies, which are generally incommensurable with ω0. Thus, the dynamics of A are more complex than those of A. In other words, A is periodic if A is stationary, and A represents a 2-torus if A is periodic. In the sequel, we consider just two attractors, one such that |A| is periodic and another such that A is quasi-periodic. Both attractors are obtained for (28)α=-10,β=10,and two values of μ. Note that since α and β are both somewhat large, dispersion and nonlinear detuning are large compared to diffusion and nonlinear damping, respectively, which is demanding from the computational point of view and promotes complex dynamics. For the construction of the databases, the equation is solved numerically using a Crank-Nicolson plus Adams-Bashforth scheme [45] for temporal discretization, with time step Δt=10-5, while the second-order spatial derivative is discretized using centered finite differences in a uniform grid of 1000 points. HODMD is applied to a set of 2000 snapshots, collected in an appropriate time interval of length 1, disregarding a transient time-stage. In the two cases considered below, a good value of the index d when applying the HODMD algorithm is d=50, while the remaining tunable parameters are ε1=10-8 and two values of ε2, namely, (29)ε2=5·10-3highaccuracy,ε2=10-4lowaccuracy.Selecting these values of the tunable parameters as common for the two cases considered below illustrates the robustness of the method, which is further checked noting that low and high accuracy computations are consistent among each other, and the results coincide when different timespans are used in the computations.The first case to be considered in that forμ=20, in which using the DMD-50 algorithm, the RRMS reconstruction error (see (16)) is ~6·10-3 and ~2·10-5 for low and high accuracy, respectively, retaining 6 and 12 DMD modes, respectively. The semilogarithmic plot of the DMD amplitude versus frequency is plotted in Figure 3-left, where it can be seen that the relevant points are approximately contained in two straight lines centered at ω0=11.6, which means that expansion (27) converges spectrally; using the standard DMD method (DMD-1 algorithm), instead, gives spurious results, plotted with red symbols in the figure. The remaining points obtained via the DMD-50 algorithm lie equispaced in frequency, with a distance between them equal to ω=118.7, which means that the frequencies ωn appearing in expansion (27) are all (positive or negative) harmonics of ω. Thus, A is time-periodic (see Figure 3, middle), with a period 2π/ω, while A is quasi-periodic (2-torus), with fundamental frequencies ω0 and ω. Since ω0 is small compared to ω the pattern is slowly modulated (Figure 3, right). The spatial modes appearing in the expansion (27), un, which are not plotted for the sake of brevity, exhibit a reflection symmetry around x=-0.17 (indicated with a vertical white line in Figure 3), meaning that the whole pattern also exhibits this symmetry, consistently with the invariance of the problem under the actions (26). This implies that the pattern can be seen as a standing wave (SW), as indicated with the vertical while solid line in Figure 3, but it can also be considered as either a right or left traveling wave (TW), as indicated with the oblique dashed lines in this figure. This dual interpretation cannot be ascertained from the shape of the modes un, but this issue will be revisited and further checked in Section 5.1.Figure 3
The CGLE forμ=20, plotting the amplitudes versus frequencies in the DMD expansion (left plot) for low (black circles) and high (black crosses) accuracy computations using DMD-50, and using DMD-1 (red symbols), and the color map of R(A) (middle plot) and |A| (right plot).As a second case, we consider the valueμ=22 (Figure 4), in which the dynamics are more complex. Again, the DMD-1 algorithm yields spurious results. Instead, for low and high accuracy (see (29)), the DMD-50 algorithm retains 17 and 49 modes, giving a RRMS reconstruction error (see (16)) ~8·10-3 and ~2·10-4, respectively. The main difference with the case considered in Figure 3 is that now the plot of the amplitude versus the frequency forms a triangular pattern that also decays spectrally but shows that, in addition to ω0=-21.2, the frequencies in the expansion in the right hand side of (27) are linear combinations with integer coefficients of the fundamental frequencies ω1=156.6 and ω2=22.2, which are not easily identified in Figure 4, left. Thus, guessing that the three fundamental frequencies are incommensurable, the attractor for |A| is contained in a 2-torus while that for A lies in a 3-torus. Finally, the middle and right plots in Figure 4 show that the pattern is a TW but, as in the former case, elucidating this issue from the shape of the modes un(x) appearing in (27) turns out to be a difficult (and unnatural) issue. The STKD method, which will be applied to the CGLE in Section 5.1, instead, will solve this issue in a natural way.Figure 4
Counterpart of Figure3 for =22.In summary, the HODMD method is able to precisely compute the temporal DMD expansion for the CGLE, giving robust results and distinguishing between periodic and quasi-periodic phenomena. However, identifying SWs and TWs requires additional ingredients, which will be included in the STKD method, developed in Section4 and applied to the CGLE in Section 5.1.


## 3.2. Application to the Zero-Net-Mass-Flux Jet

A zero-net-mass-flux (ZNMF) jet [46, 47] is a pulsatile jet promoted by a piston or membrane that interchanges fluid flow through an orifice with an outer ambient fluid that is at rest far from the orifice. The net mass flux is zero, but at large Reynolds number there is a nonzero (and large) net flux of momentum that, in addition to its own scientific interest, is useful in various fields. For instance, synthetic jet actuators promote pulsatile jets [48] that are useful in, e.g., mixing enhancement [49], jet vectoring [50], heat transfer [51], and active flow control of boundary layer separation [52]. These jets are also relevant in nature, in connection with squid, octopus, salp, and jellyfish swimming [53].It must be noted that the nonzero momentum flux at large Reynolds number is due to the fact that, even if the periodic forcing is temporally reflection symmetric, the flow is highly nonreflection symmetric, as illustrated in Figure5, where an approximately rotationally symmetric jet produced by a circular orifice is considered. As can be seen, a vortex ring is formed in the injection phase that travels downstream, while the suction phase is characterized by a saddle point that separates the fluid reentering the cavity from that still moving downstream.Figure 5
Pseudo-streamlines in a meridional section for two representative snapshots of the injection (left) and suction (right) phases in the ZNMF jet.The ZNMF jet promotes a complex flow structure that involves various, disparate spatiotemporal scales, whose understanding remains as an open topic. Thus, this is a convenient problem to test the performance of the HODMD method. Also, in order to test the noise filtering ability of the iterative HODMD method, the analysis will be based on experimental data resulting from particle image velocimetry (PIV), as obtained [46, 47] in the experimental facility for ZNMF jets of the Laboratory for Turbulence Research and Combustion at Monash University. Such measurements exhibit an error ~2.4% at the 95% confidence level. Two nondimensional parameters characterize the flow conditions, namely theReynolds and Strouhal numbers, defined as Re=UD/ν and St=fD/U, respectively, where D is the jet orifice diameter, ν is the kinematic viscosity, U is related to the peak oscillation velocity of the piston V^p as U=5V^p/2 (see [41] for details), and f is the piston oscillation frequency. In the present case, we have Re=13329 and St=0.03, which give a flow lying in the laminar-transitional regime. The flow is fully turbulent in the far field, but laminar in the near field (near the orifice), which will be the region of the jet considered here; see [41] for further details on the experimental facility and the PIV measurement.The orifice is circular and the resulting near field flow is approximately axi-symmetric, at least in the large scale. The experimental data consists inK=1872 snapshots representing 12 piston oscillation cycles, giving the streamwise and radial velocity components, u(xi,rj,tk) and v(xi,rj,tk), respectively, at I×J=52×68 grid points in a meridional rectangle (interrogation window), where x and r are the streamwise and radial coordinates, respectively, with origin at the center of the orifice and t is the time variable, with origin at the beginning of the suction phase. The time variable has been rescaled such that the (nondimensional) oscillating frequencies are measured in terms of the associated Strouhal number defined above, meaning that the piston oscillating period is T=2π/0.03. Thus, we have a snapshot tensor (instead of a matrix) and HOSVD must be used instead of standard SVD in the first step of the HODMD method (see [41] for details). Assuming permanent dynamics (zero or very small growth rates) the aim of the iterative HODMD method is to obtain an expansion of the form (30)vx,r,t≃∑n=1Nunx,reiωnt,The iterative DMD-d algorithm has been applied to obtain this expansion using two values of the index d, namely, d=1 (standard DMD) and d=700 (strictly higher order DMD), with two tolerance sets, namely, (31)ε1=ε2=2.4·10-2lowaccuracy,ε1=ε2=10-2highaccuracy,which correspond to the estimated experimental error and more precise computations, respectively. Using the latter tolerances, below the experimental error, will permit elucidating the ability of the iterative HODMD method to uncover physical phenomena that are masked by the experimental errors in the given data.Figure6 shows the amplitude versus frequency (nondimensionalized as the Strouhal number) obtained in the analyses. DMD-1 captures 11 modes, all of them spurious.Figure 6
Amplitude versus St obtained for low accuracy via DMD-1 (red crosses) and DMD-700 (blue circles) and for high accuracy using DMD-700 (black circles).On the contrary, DMD-700 capturesN=19 and N=49 DMD modes using low and high accuracy, respectively. It is remarkable that the spatial complexity M is much smaller than the spectral complexity, N, since M=5 and M=11 for low and high accuracy, respectively, which justifies the failure of DMD-1. The method shows that the flow is periodic in the near field, with dominant frequency St=0.03 (piston oscillation frequency). Since the data are real, the DMD modes appear in complex-conjugate pairs, meaning that the retained modes contain the mean flow (St=0) and 9 and 24 nonzero harmonics for low and high accuracy, respectively.The suction phase (see Figure5-left) is not quite demanding, since it just shows a saddle point at the axis of the jet. The reconstructions using DMD-700 with low and high accuracies are qualitatively similar to their counterparts for the original PIV data. The main difference is that the iterative HODMD method gives smoother results with high accuracy than with low accuracy, and both are smoother than the original PIV data. This is an indication of the ability of the method to filter experimental errors, but this issue cannot be confirmed because the actual physical data are not available. Concentrating on the reconstructions with low accuracy, it is seen that the saddle point is created near the orifice at the beginning of this phase and travels downstream to a distance that is approximately equal to 2D, where it remains steady for most part of the suction phase (Figure 7, left and middle plots). This saddle point is blown downstream at the beginning of the injection phase, where a new saddle point is created near the orifice that travels dowstream very fast (Figure 7, right plot) and leaves the near field, in a short timespan of the order of T/100. This saddle point must be present for topological reasons, but it had never been seen in studies previous to [41], just because it is quite short living. In fact, it is not present in any of the snapshots used to apply the HODMD method. It must be noted that this saddle point would be present (with some spatial noise) in the experimental snapshots if the sampling frequency was sufficiently large. Thus, the identification of this saddle point illustrates well the advantages of the temporal interpolation that is implicit in the HODMD method.Figure 7
Two snapshots in the suction phase, fort=T/18 (left plot) and t=2T/9 (middle plot), and a snapshot at the very beginning of the injection phase, for t=T/2+T/72 (right plot).The injection phase (see Figure5, right) is more interesting in the present context. Figure 8 shows a comparison of the original data with the reconstructions provided by DMD-700 for low and high accuracy. As can be seen, the reconstructions are much smoother than the original PIV data, as it happened in the suction phase. However, the reconstruction with high accuracy (recall that this is below the experimental error) is qualitatively different, since it shows some small vortices leaving the edge of the cavity that cannot be even guessed in the original PIV data and do not appear in the low accuracy reconstruction because they involve a quite small energy. These small vortices are smooth, which may indicate that they are not artifacts associated with higher order modes resulting from experimental errors. However, this cannot be ascertained since the actual physical dynamics is not available. Therefore, we have conducted some numerical simulations to elucidate this point. The value of the Reynolds number in the experimental data considered above, Re~13000, is a quite large value that would need huge numerical resources, even using the computationally efficient solver Nek5000 [54] (an open source spectral element code). Therefore, we have performed computations using the lower Reynolds number Re=1000 and the same Strouhal number as in the experimental data considered above, namely, St=0.03. The numerical counterpart of the snapshot considered in Figure 8 is given in Figure 9, which confirms the presence of the small vortices.Figure 8
Instantaneous pseudo-streamlines and streamwise velocity contours for a representative snapshot (t=T/9) in the injection phase, showing the original PIV data (left plot) and the reconstructions using DMD-700 for low (middle plot) and high (right plot) accuracy.Figure 9
Counterpart of the snapshot considered in Figure5 but as numerically computed for Re=1000.In summary, the results above illustrate well the ability of the iterative HODMD method to both (i) take advantage of temporal interpolation to uncover fast events (such as the saddle point at the beginning of the injection phase) that are not present in the snapshots and (ii) filter experimental errors so well that some small patterns (such as the small vortices in the injection phase) can be identified that were completely masked by the experimental noise. In addition, the analysis above permits constructing a purely data-driven (from experimental data) ROM that reproduces well the periodic behavior of the ZNMF jet and shows a fast online operation.


## 3.3. Additional Applications

Some additional applications are now briefly quoted.The HODMD method permits identifying flow instabilities by monitoring those DMD modes whose growth rate changes sign. On the other hand, as anticipated, retaining only those terms appearing in the DMD expansion (1) that exhibit zero (or very small) growth rate, the final attractor may be obtained from transient data. This latter idea can be used to strongly decrease the CPU time required to compute the final attractor, avoiding continuing the numerical integration until final attractor is reached, which may require an extremely large timespan, specially near bifurcation points. Transient dynamics are usually found in the initial stage of a numerical simulation, where a large number of modes (either physical or spurious) appear simultaneously.Let us illustrate these two applications by considering the cylinder wake, namely, the wake behind a circular cylinder immersed in a free stream with constant, spatially uniform velocity, which is a basic fluid dynamics problem [55, 56]. This problem depends on a unique nondimensional parameter, the Reynolds number Re=UD/ν, where U is the free stream velocity, D is the cylinder diameter, and ν is the dynamic viscosity of the fluid, assumed to be incompressible. Thus, the nondimensional velocity and pressure obey the continuity and Navier-Stokes equations: (32)∇·v=0,(33)∂v∂t+v·∇v=-∇p+1ReΔv,which are usually integrated imposing v is equal to the free stream velocity at the inlet boundary, no-slip boundary conditions at the boundary of the body, appropriate (nonreflecting) boundary conditions at the sides and exit of the computational domain, and periodic boundary conditions in the spanwise direction, with a period Lz. The wake is two-dimensional and steady for small Re, but at Re≃46 it exhibits a first unsteady instability (Hopf bifurcation [57, 58]) that produces a still two-dimensional but unsteady, periodic von Karman vortex street flow. The flow remains orbitally stable up to Re≃190, where it suffers a secondary bifurcation (Floquet multiplier=1) and becomes three-dimensional [59] for spanwise period Lz>3D (with D the cylinder diameter). The resulting mode is a long-wave, periodic, flow known asmode A. For larger values of Re, at Re=259, a second periodic short-wave mode appears if Lz>1.5D that is known asmode B. Modes A and B are illustrated in Figure 10. These modes are usually obtained via Floquet stability analysis of the basic two-dimensional periodic solution [59–61], but they can also be obtained via HODMD applied to the whole dynamics obtained via numerical simulations [62]. It is relevant to note that for higher values of the Reynolds number, at Re ≃380, a new instability takes place that is associated with a pair of complex-conjugate Floquet multipliers that give a family ofquasi-periodic solutions [61]. These modes were also identified in [62].Figure 10
Streamwise velocity contours for modes A (left) and B (right) obtained at Re=280, with and Lz=6.99.As for temporal extrapolation via the HODMD method, this idea has been used [11] precisely in the three-dimensional cylinder wake, with Re=220 and Lz=4. The numerical data was obtained using the above-mentioned numerical solver Nek5000 [54]. Temporal HODMD extrapolation divided by five the CPU time required computing the final attractor with a reasonable accuracy; see [11] for further details.Additional applications involving numerically generated data include the analysis of the turbulent wake of a cross flow wind turbine [63] and an off-shore wind turbine [64], as well as the analysis of wakes interaction [65].Concerning analyses based on experimental data, a good example, linked to the field of renewable energies, is the postprocess of light detection and ranging (LiDAR) [66, 67] measurements. LiDAR is an experimental method for the remote measurement of the line-of-sight component of wind speed, usually employed in the wind energy industry with different goals, including the maximization of the power generated by the wind turbine. Measurements are made at a few (say, M) points upstream the wind turbine, based on detection of the Doppler shift for light backscattered from natural aerosols transported by the wind in the atmosphere. LiDAR measurements offer M time-dependent signals that present high noise levels (15-20%). Thus, the spatial complexity M is small, making HODMD a suitable tool to analyze the data and detect the dominant frequencies. A great difficulty is that LiDAR experiments only measure at distances between 20 m and 200 m from the measurement device, which is a restriction in certain situations. Figure 11 shows six planes of LiDAR measurements located at distances between 33 m and 201 m upstream the wind turbine. HODMD is used to analyze such data and to predict the wind velocity at a seventh distance (228 m), which is well beyond the LiDAR device measurement range. The idea used in [68] is to extrapolate the spatial modes in the DMD expansion, guessing the data at the seventh distance with relative errors ~ 2%. The computational cost for predicting these measurements is negligible, making it possible to easily update the model in real-time (leading to a fast, purely data-driven ROM).Figure 11
Measurements upstream the wind turbine in the LiDAR experimental campaign. Data are available in the six planes that are closer to the wind turbine, while the data on the seventh plane are to be predicted.Additional applications using experimental data where the HODMD method have proved to give good results including nonisothermal flows around square objects in a channel [69] and the analysis of aeroelastic data in aircraft flight tests [70]. The aim in the latter application is to obtain the natural aeroelastic frequencies, and the associated damping rates, of the various aircraft parts (e.g., wings or fuselage) from data obtained at some sensors along the aircraft when this is subject to various types of impulsive excitation.


## 4. Spatiotemporal Koopman Decomposition

As seen in the previous sections, HODMD is suited for analyzing complex dynamics and provides accurate signal descriptions in terms of dominant frequencies and DMD modes. Nevertheless, the spatial components of the data analyzed may bring some relevant information that, combined with the temporal information, helps in the analysis and interpretation of the data. The main benefit of this idea, which gives rise to the spatiotemporal Koopman decomposition (STKD) method presented below, is that it further simplifies the spatiotemporal description of complex data in terms of TWs. The STKD method leads to data approximations of the type of (11), which for the sake of clarity is rewritten here as (34)vx,y,t=∑m,namnumnyeνm+iκmx+δn+iωnt,for x1≤x≤x1+X and t1≤t≤t1+T. This expansion can be seen as a linear combination ofmonochromatic TWs, with phase velocities cmn=-ωn/κm, which may be organized in various types of more complex standing or propagative patterns [28]. For a given finite accuracy, it is frequent that only some combinations of m and n need to be retained that, moreover, are grouped along straight lines in ω versus κ plane (called the ω-κdiagram below), which may be either horizontal or oblique and correspond to eitherstanding or traveling patterns, respectively.Pure TWs are a particular case, in which the number of spatial and temporal modes coincide, the amplitude matrixanm is diagonal, δn=νm=0, and ωn=-cκm. Then, (34) can be expressed as (35)vx-ct,y≃∑n=1Nanmunmyeiκmx-ct.In this case, the relevant points in the ω-κ diagram are in a straight line passing through the origin. If c=0, then the pattern is a pure SW. When the points in the ω-κ diagram can be organized in families of horizontal or oblique, parallel straight lines, then the pattern is a modulated SW or TW, respectively. Combinations of these two give more general patterns; see the examples presented below.In some cases, the dynamics present oscillatory behavior in more than one spatial direction. Concentrating for simplicity in the case in which two spatial directionsx and y are involved, the STKD method provides an expansion of the type (36)vx,y,z,t≃∑m1,m2,nam1m2num1m2nzeνm11+iκm11x+νm22+iκm22y+δn+iωnt,which also give TWs. For instance, if the temporal and spatial growth rates are all equal to zero and the wavenumbers and frequencies are such that (37)ωnκm112+κm222=-c,κm22=dκm11,then pattern (36) represents a pure TW with phase velocity c that propagates in the direction of the lines dy-x= constant. As in the simpler case in which only one spatial direction is involved, pure SWs correspond to the case c=0, and more general patterns are obtained combining these basic patterns. The extension to the case in which the three spatial variables are involved is straight forward.


### 4.1. The STKD Method

Let us now present the STKD method, considering the discrete counterpart of (34), with the discrete values of x and t defined as xj=x1+jΔx and tk=t1+kΔt, for j=1,…,J and k=1,…,K, with Δx=X/J and Δt=T/K; as in the temporal HODMD method, the continuous decomposition is obtained from the discrete decomposition by just allowing xj and tk to take continuous values. In addition, for simplicity, we consider a scalar state variable in the strictly one-dimensional case, in which no dependence on the transversal variable y is present. Namely, we substitute (34) by (38)vxj,tk=∑m,namnumneνm+iκmxj+δn+iωntk,where the modes umn are complex numbers with unit absolute value. The method proceeds in two steps.Step 1 (dimension reduction).
As in HODMD, truncated SVD with a tunable toleranceε1 (which determines the number of retained SVD modes, M) is applied to the J×K snapshot matrix (15), but in the present case the reduced spatial and temporal modes are scaled differently. In other words, (17) is replaced by (39)V1K≃PΣQ⊤≡X^1J⊤T^1K,where M×J and M×K matrices X1J and T1K are defined as (40)X^1J=ΣP⊤,T^1K=ΣQ⊤.
These matrices are called thereduced spatial and temporal snapshot matrices, respectively, and their columns, x^j and t^k, thereduced spatial and temporal snapshots, respectively. Comparison between (40) and (17) shows that while the SVD singular values were used in the temporal HODMD method to scale the reduced temporal snapshots, they are equidistributed now to rescale the reduced spatial and temporal snapshots, which are both scaled with the square root of the singular values.Step 2 (computation of the STKD expansion).
The higher DMD-d algorithm described in Section2 is applied to the reduced spatial and temporal snapshots, using appropriate indexes dx and dt, respectively, which leads to (41)x^j≃∑m=1Ma^mxq^mxeνm+iκmxj,forj=1,…,J,(42)t^k≃∑n=1Na^ntq^nteδn+iωntkfork=1,…,K.The numbers of retained modes, M and M, are thespectral spatial and temporal complexities, respectively, and are selected in terms of the spatial and temporal amplitudes, a^mx and a^nt, respectively, as we did in Section 2. Namely, M and N are the smallest values of the indexes m and n such that (43)a^M+1xa^1x<ε2,a^N+1ta^1t<ε2.for some tunable tolerance ε2.
Now, substituting (41)-(42) into the columns of the reduced snapshot matrices X^1J and T^1K appearing in (39) yields the STKD expansion (38), with (44)amn=q^mn>0,qmn=q^mnq^mn,
where(45)q^mn=a^mxa^ntq^mx⊤q^nt.Note that the modes qmn are complex numbers exhibiting unit absolute value, as required, and that the resulting expansion (38) exhibits M×N spatiotemporal modes. However, this number of modes can be decreased by eliminating those modes such that (46)amn<ε2,where the threshold ε2 is tunable.The method described above can be called aparallel method because the spatial and temporal variables, x and t, are treated in the same way. The extension of the method to the case in which more than one spatial variable is involved (see, e.g., (34) and (36)) is straight forward, except that the snapshot matrix becomes a snapshot tensor and HOSVD needs to be used instead of SVD in Step 1 (as in Section 3.2). For the sake of brevity, we do not describe this extension here; see [28] for a detailed development and illustration of this extension.Instead of the parallel method described above, a simpler but less consistentsequential method may be used to obtain the discrete STKD expansion (38) as follows. First, the temporal DMD-d algorithm described in Section 2, with an appropriate index dt, is applied to the snapshot matrix, which leads to (47)vxj,tk≃∑n=1Nantunjteδn+iωntk.As a second step, we apply the spatial DMD-d algorithm, with an appropriate index dx, to the N×J snapshot matrix W1J whose (n,j) element is antunjt, which leads to (48)antunj≃∑m=1Mamxumnxeνm+iκmxj.Substituting this into (47) and setting (49)amn=amxumnx,umn=amxumnxamxumnxyields the discrete expansion (38), as required. However, as in the parallel method, the number of terms in this expansion is M×N, which can be decreased by ignoring those modes that satisfy (46).Two remarks about the sequential method are now in order:(i)
The order in which the variablesx and t are treated above can be interchanged, which generally gives slightly different results.(ii)
As with the parallel method, when more than one spatial variable is present, the snapshot matrix becomes an snapshot tensor and standard SVD must be replaced by HOSVD in the application of the HODMD method.Because the parallel method is more consistent, specially when treating more than one spatial variable, it will be this method that will be used in the remainder of this section, unless otherwise stated.As with the HODMD, the STKD method can be used to analyze experimental data. The noisy artifacts that may cover the actual physical data can be cleaned using the proper tolerances and applying the method iteratively (as in the example presented in Section3.2). Similarly, the STKD method is suitable for analyzing transient dynamics (as in the example shown in Section 3.3).


### 4.2. Illustration of the STKD Method in a Toy Model

Let us now apply these two methods in the following toy model, already considered in [29], defined as (50)ux,t=0.5+sin⁡x2cos⁡κ~1x-ω~1t+0.5cos⁡κ~2x-ω~2t,with κ~1=2, κ~2=10, ω~1=2π/45, and ω~2=10. It is a simple (but demanding because ω~1 and ω~2 are disparate from each other; see Figure 12-left) model that can be cast into the form (34), with νm=δn=0, involving 12 wavenumbers, namely, ±κ~1,±κ~2,±(κ~1±1),±(κ~2±1), and 4 frequencies, namely, ±ω~1 and ±ω~2. Thus the spectral spatial and temporal complexities are M=12 and N=4, respectively. The pattern is temporally quasi-periodic because ω~1 and ω~2 are incommensurable, but spatially periodic, with a period equal to 2π.Figure 12
Color map of the toy model (50) (left plot) and its restriction to the sampled spatiotemporal square where the STKD method is applied (right plot).The STKD method is applied in the sampled square0≤x≤10,0≤t≤5 (see Figure 12, right), where x and t are discretized using 50 and 25 points, respectively. The tunable parameters of the method are ε1=10-8, ε2=10-7, dx=1, and dt=10. The obtained ω-κ diagram is given in Figure 13, which shows that the pattern is a modulated TW, obtained as a linear superposition of two pure TWs, whose phase velocities are c1=ω1/κ1=π/45 and c2=ω2/κ2=1/10. The method provides the involved 12 wavenumbers and 4 frequencies with zero machine accuracy, while the RRMS error (see (16)) of the reconstruction is ~10-14. It is remarkable that these extremely good results are obtained using a sampled spatial interval that is somewhat larger to the spatial period, but a temporal sampled interval that is much smaller than the largest temporal period, which is 2π/ω~1=45. Also note that the extrapolation from the sampled time interval 0≤t≤5 (Figure 12, right) to the larger time interval 0≤t≤100 (Figure 12, left) still produces a good reconstruction, whose RRMS error, as defined in (16), is ~10-14.Figure 13
ω - κ diagram of the toy model (50).The sequential method has also been applied to analyze this toy model, obtaining results that are similar to those described above.


## 4.1. The STKD Method

Let us now present the STKD method, considering the discrete counterpart of (34), with the discrete values of x and t defined as xj=x1+jΔx and tk=t1+kΔt, for j=1,…,J and k=1,…,K, with Δx=X/J and Δt=T/K; as in the temporal HODMD method, the continuous decomposition is obtained from the discrete decomposition by just allowing xj and tk to take continuous values. In addition, for simplicity, we consider a scalar state variable in the strictly one-dimensional case, in which no dependence on the transversal variable y is present. Namely, we substitute (34) by (38)vxj,tk=∑m,namnumneνm+iκmxj+δn+iωntk,where the modes umn are complex numbers with unit absolute value. The method proceeds in two steps.Step 1 (dimension reduction).
As in HODMD, truncated SVD with a tunable toleranceε1 (which determines the number of retained SVD modes, M) is applied to the J×K snapshot matrix (15), but in the present case the reduced spatial and temporal modes are scaled differently. In other words, (17) is replaced by (39)V1K≃PΣQ⊤≡X^1J⊤T^1K,where M×J and M×K matrices X1J and T1K are defined as (40)X^1J=ΣP⊤,T^1K=ΣQ⊤.
These matrices are called thereduced spatial and temporal snapshot matrices, respectively, and their columns, x^j and t^k, thereduced spatial and temporal snapshots, respectively. Comparison between (40) and (17) shows that while the SVD singular values were used in the temporal HODMD method to scale the reduced temporal snapshots, they are equidistributed now to rescale the reduced spatial and temporal snapshots, which are both scaled with the square root of the singular values.Step 2 (computation of the STKD expansion).
The higher DMD-d algorithm described in Section2 is applied to the reduced spatial and temporal snapshots, using appropriate indexes dx and dt, respectively, which leads to (41)x^j≃∑m=1Ma^mxq^mxeνm+iκmxj,forj=1,…,J,(42)t^k≃∑n=1Na^ntq^nteδn+iωntkfork=1,…,K.The numbers of retained modes, M and M, are thespectral spatial and temporal complexities, respectively, and are selected in terms of the spatial and temporal amplitudes, a^mx and a^nt, respectively, as we did in Section 2. Namely, M and N are the smallest values of the indexes m and n such that (43)a^M+1xa^1x<ε2,a^N+1ta^1t<ε2.for some tunable tolerance ε2.
Now, substituting (41)-(42) into the columns of the reduced snapshot matrices X^1J and T^1K appearing in (39) yields the STKD expansion (38), with (44)amn=q^mn>0,qmn=q^mnq^mn,
where(45)q^mn=a^mxa^ntq^mx⊤q^nt.Note that the modes qmn are complex numbers exhibiting unit absolute value, as required, and that the resulting expansion (38) exhibits M×N spatiotemporal modes. However, this number of modes can be decreased by eliminating those modes such that (46)amn<ε2,where the threshold ε2 is tunable.The method described above can be called aparallel method because the spatial and temporal variables, x and t, are treated in the same way. The extension of the method to the case in which more than one spatial variable is involved (see, e.g., (34) and (36)) is straight forward, except that the snapshot matrix becomes a snapshot tensor and HOSVD needs to be used instead of SVD in Step 1 (as in Section 3.2). For the sake of brevity, we do not describe this extension here; see [28] for a detailed development and illustration of this extension.Instead of the parallel method described above, a simpler but less consistentsequential method may be used to obtain the discrete STKD expansion (38) as follows. First, the temporal DMD-d algorithm described in Section 2, with an appropriate index dt, is applied to the snapshot matrix, which leads to (47)vxj,tk≃∑n=1Nantunjteδn+iωntk.As a second step, we apply the spatial DMD-d algorithm, with an appropriate index dx, to the N×J snapshot matrix W1J whose (n,j) element is antunjt, which leads to (48)antunj≃∑m=1Mamxumnxeνm+iκmxj.Substituting this into (47) and setting (49)amn=amxumnx,umn=amxumnxamxumnxyields the discrete expansion (38), as required. However, as in the parallel method, the number of terms in this expansion is M×N, which can be decreased by ignoring those modes that satisfy (46).Two remarks about the sequential method are now in order:(i)
The order in which the variablesx and t are treated above can be interchanged, which generally gives slightly different results.(ii)
As with the parallel method, when more than one spatial variable is present, the snapshot matrix becomes an snapshot tensor and standard SVD must be replaced by HOSVD in the application of the HODMD method.Because the parallel method is more consistent, specially when treating more than one spatial variable, it will be this method that will be used in the remainder of this section, unless otherwise stated.As with the HODMD, the STKD method can be used to analyze experimental data. The noisy artifacts that may cover the actual physical data can be cleaned using the proper tolerances and applying the method iteratively (as in the example presented in Section3.2). Similarly, the STKD method is suitable for analyzing transient dynamics (as in the example shown in Section 3.3).


## 4.2. Illustration of the STKD Method in a Toy Model

Let us now apply these two methods in the following toy model, already considered in [29], defined as (50)ux,t=0.5+sin⁡x2cos⁡κ~1x-ω~1t+0.5cos⁡κ~2x-ω~2t,with κ~1=2, κ~2=10, ω~1=2π/45, and ω~2=10. It is a simple (but demanding because ω~1 and ω~2 are disparate from each other; see Figure 12-left) model that can be cast into the form (34), with νm=δn=0, involving 12 wavenumbers, namely, ±κ~1,±κ~2,±(κ~1±1),±(κ~2±1), and 4 frequencies, namely, ±ω~1 and ±ω~2. Thus the spectral spatial and temporal complexities are M=12 and N=4, respectively. The pattern is temporally quasi-periodic because ω~1 and ω~2 are incommensurable, but spatially periodic, with a period equal to 2π.Figure 12
Color map of the toy model (50) (left plot) and its restriction to the sampled spatiotemporal square where the STKD method is applied (right plot).The STKD method is applied in the sampled square0≤x≤10,0≤t≤5 (see Figure 12, right), where x and t are discretized using 50 and 25 points, respectively. The tunable parameters of the method are ε1=10-8, ε2=10-7, dx=1, and dt=10. The obtained ω-κ diagram is given in Figure 13, which shows that the pattern is a modulated TW, obtained as a linear superposition of two pure TWs, whose phase velocities are c1=ω1/κ1=π/45 and c2=ω2/κ2=1/10. The method provides the involved 12 wavenumbers and 4 frequencies with zero machine accuracy, while the RRMS error (see (16)) of the reconstruction is ~10-14. It is remarkable that these extremely good results are obtained using a sampled spatial interval that is somewhat larger to the spatial period, but a temporal sampled interval that is much smaller than the largest temporal period, which is 2π/ω~1=45. Also note that the extrapolation from the sampled time interval 0≤t≤5 (Figure 12, right) to the larger time interval 0≤t≤100 (Figure 12, left) still produces a good reconstruction, whose RRMS error, as defined in (16), is ~10-14.Figure 13
ω - κ diagram of the toy model (50).The sequential method has also been applied to analyze this toy model, obtaining results that are similar to those described above.


## 5. Applications of STKD

Let us now illustrate in some applications the ability of the STKD method to uncover the spatiotemporal structure of complex dynamics. Among these applications, we consider the CGLE (already analyzed with the temporal HODMD method in Section3.1) and the dynamics of the thermal convection in a rotating spherical shell subject to gravitational force.


### 5.1. Application of the STKD Method to the CGLE

For the CGLE, the STKD counterpart of (27) is (51)Ax,t=eiω0t∑mnamnumneiκmx+ωnt,where the modes umn are now unit complex numbers. Concentrating on the two cases analyzed in Section 3.1 and Figures 3 and 4, we use the same snapshots and tunable parameters used in Section 3.1, including the index dt, which is set as being equal to the index d used in Section 3.1. The index dx, instead, is set as dx=1. In both cases, the RRMS reconstruction errors, as defined in (16), are comparable to their counterparts obtained in Section 3.1 using the temporal HODMD method, retaining 25 and 71 spatiotemporal modes for low and high accuracy, respectively, in the case considered in Figure 3 and 17 and 49 modes for low and high accuracy, respectively, in the case considered in Figure 4. The ω-κ diagrams for these two cases are given in Figure 14. The left plot in this figure shows that points can be conformed in either horizontal or oblique straight lines, which clearly indicates that the pattern can be considered as either a modulated SW or a modulated TW, consistent with what is clearly seen in Figure 3. As anticipated in Section 3.1, the value of ω0 in the expansion (51) equals 11.6, while the frequencies ωn appearing in the expansion in the right hand side of (51) are positive and negative harmonics of the fundamental frequency ω1=118.7, which is associated with the wavenumber κ1=π. Thus, this pattern is slowly modulated (because ω0 is small compared to ω1) and when it is seen as a TW, the phase velocity is c=±ω1/κ1=±37.8.Figure 14
ω - κ diagrams resulting from the STKD method applied to the CGLE for the cases considered in Figure 3 (left plot) and Figure 4 (right plot). As in Figure 3, left, and Figure 4, left, crosses and circles label results obtained with high and low accuracies, respectively.Concerning the right plot in Figure14, the points in this plot are conformed in the indicated oblique straight lines, which means that the pattern is a modulated TW. As in Section 3.1, the frequency ω0 appearing in (51) equals -21.2, while the various frequencies appearing in the expansion in the right hand side of this equations are all combinations of the fundamental frequencies ω1=156.6 and ω2=22.2, whose associated wavenumbers are κ1=π and κ2=2π, respectively; the associated phase velocities are c1=-ω1/κ1=-49.8 and c2=-ω2/κ2=-3.53, which are indicated in Figures 3 and 4, middle and right plots.In summary, the STKD method computes quite well the frequencies and wavenumbers involved in SWs and TWs for the CGLE and permits straightforward identification of the associated phase velocities.


### 5.2. STKD in a Three-Dimensional Spherical Shell with Radial Thermal Convection, Subject to Radial Gravity

Let us consider the three-dimensional thermal convection in a rotating spherical shell, heated from the inside, subject to a radial gravitational force. This problem is relevant in geophysical and astrophysical fluid dynamics [71], in connection with the transport of mass and energy in the upper atmospheres of stars and giant planets, where it is specially interesting the identification of periodic and quasi-periodic rotating waves promoted by convective instabilities.The continuity, momentum, and energy equations governing the dynamics of the flow, using the Boussinesq approximation in a rotating (with angular velocityΩ) frame linked to the spherical shell, are (52)∇·v,(53)∂tv+v·∇v+2E-1k×v=-∇Π+∇2v+Θr,(54)σ∂tΘ+v·∇Θ=∇2Θ+Rη1-η-2r-3r·v.Here, v, Π, and Θ are the velocity vector, the pressure, and the perturbation of the temperature from the quiescent purely convective state, respectively, r=(x,y,z) is the position vector referred to a Cartesian coordinate system, with the origin at the center of the sphere and the z axis along the axis of rotation, r=|r|, and k is the unit vector along the z axis. The boundary conditions are homogeneous Dirichlet for both the velocity and temperature at the inner and outer spheres.The problem depends on threenondimensional parameters, thePrandtl numberσ=ν/κ (where ν and κ are the kinematic viscosity and thermal conductivity, respectively), theEckman numberE=ν/(Ωd2) (where d=ro∗-ri∗ is the radial gap between the spheres), theRayleigh numberR=γαΔTd4/(νκ) (where α is the thermal expansion coefficient, γri∗ is the imposed radial gravity, and ΔT is increment in temperature), and theinner to outer radii ratioη=ri∗/ro∗. This problem is invariant under the symmetry group SO(2)×Z2, generated by rotations about the z axis and up-down reflection on the equatorial plane. Invariance under rotation is essential to obtain TWs (which are rotating waves along the azimuthal direction in the present context), most of which are also (instantaneously) invariant under the up/down reflection symmetry.The computational domain is composed by 25, 128, and 32 grid points in the radial, colatitudinal, and azimuthal directions, respectively. Thenumerical solver used to solve this problem is described in [72–75], where the problem is reformulated by taking the curl of the momentum equation (which eliminates the pressure) and writing the solenoidal velocity field (which eliminates the continuity equation) as (55)v=∇×Ψr+∇×∇×Φr,where Ψ and Φ are thetoroidal and poloidal scalar potentials, respectively, first introduced by Chandrasekhar [76]. The numerical data used in the analysis below have been provided to us by Professors Marta Net and Joan Sánchez. These data correspond to E=10-4, η=0.35, σ=0.1, and two values of the Rayleigh number, namely, R=8·105 and R=1.1·106. In both cases, the STKD method is applied to the temperature field, which gives the expansion (56)Θrj1,θj2,ϕ,t≃∑m,n=1M,NamnΘmnj1j2eikmϕ+ωnt,where rj1 and θj2 denote the mesh points in the radial and colatitudinal directions, respectively, and ϕ denotes the azimuthal coordinate. As in the previous applications of the STKD method, two sets of values of the tunable parameters are considered, namely, (ε1,ε2)=(10-4,5·10-3) and (10-5,5·10-4), which will be referred to as low and high accuracy, respectively. The index d may be selected as d=1 for both the spatial and temporal applications of the HODMD method that are needed in the second step of the STKD method.ForR=8·105, the STKD method is applied considering 100 equispaced snapshots in a time interval of length 0.3, selected after sufficiently large timespan to eliminate transient behavior. The RRMS reconstruction error, as defined in (16), is 1.6·10-3 and 1.5·10-5 for low and high accuracies, respectively, retaining 5 and 9 spatiotemporal modes, respectively. The ω-κ diagram is as given in Figure 15, which shows that for both low and high accuracy; the relevant points (whose frequencies and wavenumbers are multiples of ω1=45.5 and κ1=8, respectively) are in a straight line passing through the origin. As explained in Section 4, this means that the pattern is a periodic, pure TW; namely, it is steady in a rotating frame with phase speed, c=-ω1/κ1=-5.69. Such a steady solution is considered in Figure 16, where the fundamental wavenumber κ1=8 is clearly seen in the equatorial section ϕ=π/2.Figure 15
Theω-κ diagram (with circles and crosses indicating low and high accuracy) for the attractor at R=8⋅105.Figure 16
A snapshot showing the contour levels ofΘ at the indicated surfaces for R=8⋅105.Concerning the caseR=1.1·106, in which the dynamics are more complex than in the former case, the STKD method is applied considering 500 equispaced snapshots in a time interval of length 1.4, again selected after sufficiently large timespan to eliminate transient behavior. Now, the RRMS reconstruction error, as defined in (16), is 8·10-3 and 1.4·10-4 for low and high accuracies, respectively, retaining 17 and 47 spatiotemporal modes, respectively. The ω-κ diagram is as given in Figure 17, where it can be seen that the relevant points are aligned along oblique straight lines. The fundamental frequencies and wavenumbers are (ω1,κ1)=(47.6,8) and (-239.8,6), while the remaining points are linear combinations (with integer coefficients) of these. Thus, the pattern is a modulated TW, which is periodic in a frame moving with either of the phase velocities c1=ω1/κ1=-5.95 or c2=ω2/κ2=40.Figure 17
Counterpart of Figure15R=1.1⋅106.Since the counterpart of the space-time diagrams presented for the CGLE in Section3.1 would lead to four-dimensional plots that are not possible in the present case, the character of the pattern as a modulated TW is illustrated in Figure 18 considering the restricted space-time diagrams (along the azimuthal coordinate) in two equatorial circles. As can be seen in these plots, the modulated TW character (with the already identified phase velocities, indicated here with straight lines) of the pattern is clearly seen.Figure 18
The quasi-periodic TW forR=1.1⋅106: color maps of Θ along the azimuthal coordinate at the indicated equatorial circles.


### 5.3. Additional Applications

The STKD method has already been used in the analysis of TWs appearing in the following:(i)
Off-shore wind turbines [64], whose structure is of great technological interest(ii)
The three-dimensional cylinder wake [77], which permits identifying well the TWs associated with the already mentioned periodic modes A and B, and also with some quasi-periodic modes. These modes play a fundamental role in the associated dynamics [61]Additional ongoing applications include the following:(i)
Interacting wakes of bluff bodies, which generally gives complex quasi-periodic phenomena, including at least the frequencies and wavenumbers associated with the individual wakes(ii)
Several turbomachinery configurations(iii)
The reconstruction the 3D flow from various (PIV) 2D measurements, taking advantage of the fact that if the 3D flow is given by (34), then the whole unsteady 3D flow field may be reconstructed if the (steady) 2D modes umn(y) are available. These can be computed by fitting with unsteady 2D measurements. Note that this is a demanding application of great scientific and technological interest(iv)
The identification of the large scale dynamics in fully turbulent flows from either experimental (PIV) or numerical data resulting from direct numerical simulation. Again, this is a demanding application of great scientific interest


## 5.1. Application of the STKD Method to the CGLE

For the CGLE, the STKD counterpart of (27) is (51)Ax,t=eiω0t∑mnamnumneiκmx+ωnt,where the modes umn are now unit complex numbers. Concentrating on the two cases analyzed in Section 3.1 and Figures 3 and 4, we use the same snapshots and tunable parameters used in Section 3.1, including the index dt, which is set as being equal to the index d used in Section 3.1. The index dx, instead, is set as dx=1. In both cases, the RRMS reconstruction errors, as defined in (16), are comparable to their counterparts obtained in Section 3.1 using the temporal HODMD method, retaining 25 and 71 spatiotemporal modes for low and high accuracy, respectively, in the case considered in Figure 3 and 17 and 49 modes for low and high accuracy, respectively, in the case considered in Figure 4. The ω-κ diagrams for these two cases are given in Figure 14. The left plot in this figure shows that points can be conformed in either horizontal or oblique straight lines, which clearly indicates that the pattern can be considered as either a modulated SW or a modulated TW, consistent with what is clearly seen in Figure 3. As anticipated in Section 3.1, the value of ω0 in the expansion (51) equals 11.6, while the frequencies ωn appearing in the expansion in the right hand side of (51) are positive and negative harmonics of the fundamental frequency ω1=118.7, which is associated with the wavenumber κ1=π. Thus, this pattern is slowly modulated (because ω0 is small compared to ω1) and when it is seen as a TW, the phase velocity is c=±ω1/κ1=±37.8.Figure 14
ω - κ diagrams resulting from the STKD method applied to the CGLE for the cases considered in Figure 3 (left plot) and Figure 4 (right plot). As in Figure 3, left, and Figure 4, left, crosses and circles label results obtained with high and low accuracies, respectively.Concerning the right plot in Figure14, the points in this plot are conformed in the indicated oblique straight lines, which means that the pattern is a modulated TW. As in Section 3.1, the frequency ω0 appearing in (51) equals -21.2, while the various frequencies appearing in the expansion in the right hand side of this equations are all combinations of the fundamental frequencies ω1=156.6 and ω2=22.2, whose associated wavenumbers are κ1=π and κ2=2π, respectively; the associated phase velocities are c1=-ω1/κ1=-49.8 and c2=-ω2/κ2=-3.53, which are indicated in Figures 3 and 4, middle and right plots.In summary, the STKD method computes quite well the frequencies and wavenumbers involved in SWs and TWs for the CGLE and permits straightforward identification of the associated phase velocities.


## 5.2. STKD in a Three-Dimensional Spherical Shell with Radial Thermal Convection, Subject to Radial Gravity

Let us consider the three-dimensional thermal convection in a rotating spherical shell, heated from the inside, subject to a radial gravitational force. This problem is relevant in geophysical and astrophysical fluid dynamics [71], in connection with the transport of mass and energy in the upper atmospheres of stars and giant planets, where it is specially interesting the identification of periodic and quasi-periodic rotating waves promoted by convective instabilities.The continuity, momentum, and energy equations governing the dynamics of the flow, using the Boussinesq approximation in a rotating (with angular velocityΩ) frame linked to the spherical shell, are (52)∇·v,(53)∂tv+v·∇v+2E-1k×v=-∇Π+∇2v+Θr,(54)σ∂tΘ+v·∇Θ=∇2Θ+Rη1-η-2r-3r·v.Here, v, Π, and Θ are the velocity vector, the pressure, and the perturbation of the temperature from the quiescent purely convective state, respectively, r=(x,y,z) is the position vector referred to a Cartesian coordinate system, with the origin at the center of the sphere and the z axis along the axis of rotation, r=|r|, and k is the unit vector along the z axis. The boundary conditions are homogeneous Dirichlet for both the velocity and temperature at the inner and outer spheres.The problem depends on threenondimensional parameters, thePrandtl numberσ=ν/κ (where ν and κ are the kinematic viscosity and thermal conductivity, respectively), theEckman numberE=ν/(Ωd2) (where d=ro∗-ri∗ is the radial gap between the spheres), theRayleigh numberR=γαΔTd4/(νκ) (where α is the thermal expansion coefficient, γri∗ is the imposed radial gravity, and ΔT is increment in temperature), and theinner to outer radii ratioη=ri∗/ro∗. This problem is invariant under the symmetry group SO(2)×Z2, generated by rotations about the z axis and up-down reflection on the equatorial plane. Invariance under rotation is essential to obtain TWs (which are rotating waves along the azimuthal direction in the present context), most of which are also (instantaneously) invariant under the up/down reflection symmetry.The computational domain is composed by 25, 128, and 32 grid points in the radial, colatitudinal, and azimuthal directions, respectively. Thenumerical solver used to solve this problem is described in [72–75], where the problem is reformulated by taking the curl of the momentum equation (which eliminates the pressure) and writing the solenoidal velocity field (which eliminates the continuity equation) as (55)v=∇×Ψr+∇×∇×Φr,where Ψ and Φ are thetoroidal and poloidal scalar potentials, respectively, first introduced by Chandrasekhar [76]. The numerical data used in the analysis below have been provided to us by Professors Marta Net and Joan Sánchez. These data correspond to E=10-4, η=0.35, σ=0.1, and two values of the Rayleigh number, namely, R=8·105 and R=1.1·106. In both cases, the STKD method is applied to the temperature field, which gives the expansion (56)Θrj1,θj2,ϕ,t≃∑m,n=1M,NamnΘmnj1j2eikmϕ+ωnt,where rj1 and θj2 denote the mesh points in the radial and colatitudinal directions, respectively, and ϕ denotes the azimuthal coordinate. As in the previous applications of the STKD method, two sets of values of the tunable parameters are considered, namely, (ε1,ε2)=(10-4,5·10-3) and (10-5,5·10-4), which will be referred to as low and high accuracy, respectively. The index d may be selected as d=1 for both the spatial and temporal applications of the HODMD method that are needed in the second step of the STKD method.ForR=8·105, the STKD method is applied considering 100 equispaced snapshots in a time interval of length 0.3, selected after sufficiently large timespan to eliminate transient behavior. The RRMS reconstruction error, as defined in (16), is 1.6·10-3 and 1.5·10-5 for low and high accuracies, respectively, retaining 5 and 9 spatiotemporal modes, respectively. The ω-κ diagram is as given in Figure 15, which shows that for both low and high accuracy; the relevant points (whose frequencies and wavenumbers are multiples of ω1=45.5 and κ1=8, respectively) are in a straight line passing through the origin. As explained in Section 4, this means that the pattern is a periodic, pure TW; namely, it is steady in a rotating frame with phase speed, c=-ω1/κ1=-5.69. Such a steady solution is considered in Figure 16, where the fundamental wavenumber κ1=8 is clearly seen in the equatorial section ϕ=π/2.Figure 15
Theω-κ diagram (with circles and crosses indicating low and high accuracy) for the attractor at R=8⋅105.Figure 16
A snapshot showing the contour levels ofΘ at the indicated surfaces for R=8⋅105.Concerning the caseR=1.1·106, in which the dynamics are more complex than in the former case, the STKD method is applied considering 500 equispaced snapshots in a time interval of length 1.4, again selected after sufficiently large timespan to eliminate transient behavior. Now, the RRMS reconstruction error, as defined in (16), is 8·10-3 and 1.4·10-4 for low and high accuracies, respectively, retaining 17 and 47 spatiotemporal modes, respectively. The ω-κ diagram is as given in Figure 17, where it can be seen that the relevant points are aligned along oblique straight lines. The fundamental frequencies and wavenumbers are (ω1,κ1)=(47.6,8) and (-239.8,6), while the remaining points are linear combinations (with integer coefficients) of these. Thus, the pattern is a modulated TW, which is periodic in a frame moving with either of the phase velocities c1=ω1/κ1=-5.95 or c2=ω2/κ2=40.Figure 17
Counterpart of Figure15R=1.1⋅106.Since the counterpart of the space-time diagrams presented for the CGLE in Section3.1 would lead to four-dimensional plots that are not possible in the present case, the character of the pattern as a modulated TW is illustrated in Figure 18 considering the restricted space-time diagrams (along the azimuthal coordinate) in two equatorial circles. As can be seen in these plots, the modulated TW character (with the already identified phase velocities, indicated here with straight lines) of the pattern is clearly seen.Figure 18
The quasi-periodic TW forR=1.1⋅106: color maps of Θ along the azimuthal coordinate at the indicated equatorial circles.


## 5.3. Additional Applications

The STKD method has already been used in the analysis of TWs appearing in the following:(i)
Off-shore wind turbines [64], whose structure is of great technological interest(ii)
The three-dimensional cylinder wake [77], which permits identifying well the TWs associated with the already mentioned periodic modes A and B, and also with some quasi-periodic modes. These modes play a fundamental role in the associated dynamics [61]Additional ongoing applications include the following:(i)
Interacting wakes of bluff bodies, which generally gives complex quasi-periodic phenomena, including at least the frequencies and wavenumbers associated with the individual wakes(ii)
Several turbomachinery configurations(iii)
The reconstruction the 3D flow from various (PIV) 2D measurements, taking advantage of the fact that if the 3D flow is given by (34), then the whole unsteady 3D flow field may be reconstructed if the (steady) 2D modes umn(y) are available. These can be computed by fitting with unsteady 2D measurements. Note that this is a demanding application of great scientific and technological interest(iv)
The identification of the large scale dynamics in fully turbulent flows from either experimental (PIV) or numerical data resulting from direct numerical simulation. Again, this is a demanding application of great scientific interest


## 6. Conclusions

A review has been presented on the purely data-driven HODMD and STKD methods and their multiple applications. These methods decompose either temporal or spatiotemporal experimental or numerical signals as a combination of modes. In particular, HODMD is an extension of standard DMD that has general validity, namely it provides good results (identifying the relevant temporal growth rates and frequencies) in cases in which standard DMD fails. STKD also computes the relevant spatial growth rates and wavenumbers, and providesω-κ diagrams that permit identifying TWs or SWs in the given data. These methods allow for two main tasks:(i)
Uncovering the physics underlying the given data. In particular, in the case of experimental data, these methods are able to filter noise quite efficiently and extract patterns that were completely masked by experimental noise.(ii)
Constructing ROMs whose online operation is very fast. This is because the HODMD and STKD reconstructions involve semianalytic descriptions that only require a limited amount of operations.The HODMD and STKD methods have been illustrated in various applications of basic scientific interest, including the CGLE, the ZNMF jet, the thermal convection in a spherical rotating shell, and the wake of a circular cylinder. Several, more engineering oriented applications dealing with, e.g., aircraft flight tests and wind turbine problems, have also been considered.



---

*Source: 6920783-2018-12-12.xml*