# A robust dynamic mode decomposition methodology for an airfoil undergoing
transonic shock buffet

Andre Weiner111Post-doctoral researcher, flow modeling and control group,
<EMAIL_ADDRESS>and Richard Semaan.222Group leader, flow modeling
and control group<EMAIL_ADDRESS>Technische Universität
Braunschweig, Institute of Fluid Mechanics, Braunschweig, 38108, Germany

###### Abstract

Dynamic mode decomposition (DMD) is a data-driven technique widely used to
analyze and model fluid problems including transonic buffet flows. Despite its
strengths, DMD is known to suffer from sensitivities to the selected settings
and the characteristics of the employed data. In this work, we closely examine
the aforementioned sensitivities, identify possible pitfalls, and provide best
practices to robustly perform DMD on a flow exhibiting transonic shock buffet.
Specifically, we assess various DMD variants and test their sensitivity to the
POD rank truncation, the sampling rate, and the time window. A critical
enabler to our analysis is a new presentation of the DMD algorithm as a
modular framework consisting of five distinct steps. The test results show a
high sensitivity to rank truncation and a low sensitivity to the time window
selection when a reasonable time period is used. The tests also highlight the
existing dangers of aliasing, when the sampling rate is too low. Finally, a
list of practical recommendations and guidelines on how to accurately and
robustly perform DMD on a transonic buffet flow is provided.

## Nomenclature

$a$ = | local speed of sound
---|---
$a_{\infty}$ = | freestream speed of sound
$\mathbf{A}$ = | discrete linear operator
$\mathbf{\mathcal{A}}$ = | continuous linear operator
$\tilde{\mathbf{A}}$ = | reduced discrete linear operator
$\mathbf{b}$ = | DMD mode amplitudes
$c$ = | chord length
$c_{l}$ = | lift coefficient
$c_{p}$ = | pressure coefficient
$Co$ = | Courant number
$\mathbf{D}_{\mathbf{b}}$ = | diagonal matrix based on $\mathbf{b}$
$e$ = | total energy density
$f$ = | frequency
$f_{b}$ = | buffet frequency
$f_{s}$ = | snapshot sampling frequency
$\tilde{f}$ = | normalized frequency
$I_{i}$ = | importance of mode $i$
$k$ = | kinematic energy density
$\mathbf{L}$ = | low-rank contributions to $\mathbf{M}$
$\mathbf{M}$ = | data matrix containing all snapshots
$\mathrm{Ma}$ = | local Mach number
$\mathrm{Ma}_{\infty}$ = | freestream Mach number
$N$ = | number of snapshots
$N_{\mathrm{cells}}$ = | number of control volumes
$N_{\mathrm{cycle}}$ = | number of buffet cycles
$r$ = | matrix rank
$r_{\mathrm{opt}}$ = | optimal matrix rank
$Re_{\infty}$ = | freestream Reynolds number
$\mathbf{S}$ = | sparse contributions to $\mathbf{M}$
$t$ = | time
$\tilde{t}$ = | normalized time
$\hat{u}$ = | internal energy density
$\mathbf{u}$ = | velocity vector
$u_{\tau}$ = | friction velocity
$\mathbf{U}$ = | matrix of left singular vectors
$U_{\infty}$ = | freestream velocity
$\mathbf{V}$ = | matrix of right singular vectors
$\mathbf{V}_{\mathbf{\omega}}$ = | Vandermode matrix based on $\mathbf{\omega}$
$w$ = | weighting factor based on cell geometry
$\mathbf{W}$ = | matrix of eigenvectors of $\tilde{\mathbf{A}}$
$\tilde{x}$ = |  $x$-coordinate divided by chord length
$\mathbf{x}_{n}$ = | snapshot of state vector at time step $n$
$\mathbf{X}$ = | data matrix of all but the last snapshots
$\mathbf{X}_{r}$ = | truncated generic matrix $\mathbf{X}$
$\tilde{y}$ = |  $y$-coordinate divided by chord length
$y^{+}$ = | normalized coordinate normal to surface
$\mathbf{Y}$ = | data matrix with all but the first snapshots
$\alpha$ = | angle of attack
$\gamma$ = | heat capacity ratio
$\gamma_{0}$ = | sparsity control parameter in the spDMD
$\Delta t$ = | time step
$\Delta\mathbf{X}$ = | error in $\mathbf{X}$
$\Delta\mathbf{Y}$ = | error in $\mathbf{Y}$
$\lambda$ = | eigenvalue of $\mathbf{A}$
$\lambda_{0}$ = | sparsity control parameter in the robust SVD
$\mathbf{\Lambda}$ = | matrix of eigenvalues of $\mathbf{A}$
$\mu_{x}$ = | temporal mean of generic $x$
$\rho$ = | density
$\sigma$ = | singular value
$\sigma_{x}$ = | temporal standard deviation of generic $x$
$\mathbf{\Sigma}$ = | matrix of singular values
$\tilde{\tau}$ = | normalized correlation time lag
$\nu$ = | kinematic viscosity
$\phi$ = | DMD phase angle
$\mathbf{\Phi}$ = | DMD modes, eigenvectors of $\mathbf{A}$
$\mathbf{\Phi}_{\mathbf{b}}$ = |  $\mathbf{\Phi}$ with columns normalized by $\mathbf{b}$
$\mathbf{\omega}$ = | trace of $\mathbf{\Omega}$
$\mathbf{\Omega}$ = | matrix of eigenvalues of $\mathbf{\mathcal{A}}$
$\Omega$ = | domain
DES = | Detached Eddy Simulation
DDES = | Delayed DES
DMD = | Dynamic Mode Decomposition
IDDES = | Improved DDES
iPSP = | instationary Pressure Sensitive Paint
LES = | Large Eddy Simulation
mrDMD = | multi resolution DMD
PIV = | Particle Image Velocimetry
POD = | Proper Orthogonal Decomposition
PSD = | Power Spectral Density
SA = | Spalart-Allmaras
SST = | Shear Stress Transport
SVD = | Singular Value Decomposition
spDMD = | sparsity-promoting DMD
TDMD = | Total least-squares DMD
RANS = | Reynolds Averaged Navier Stokes
UDMD = | Unitary DMD
URANS = | Unsteady RANS
ZDES = | Zonal DES

##  1 Introduction

The transonic shock buffet is a flow phenomenon characterized by self-
sustained oscillations caused by interactions between the shock wave and the
boundary layer over an airfoil in transonic flow conditions. A good
understanding of the shock buffet enables predictions of the flight envelope
and, possibly, a mitigation of the phenomenon by means of flow control.
Numerous experimental and numerical studies have been conducted to deepen our
understanding of the buffet mechanism. A comprehensive overview of transonic
shock buffet over airfoils and wings is provided by Giannelis et al. [1]. One
of the open questions regarding the understanding of two-dimensional transonic
buffet is the nature of the interaction between shock motion, boundary layer
separation, trailing edge vortex shedding, and acoustic waves propagating in
the upstream direction. A better understanding of the underlying flow physics
requires accurate and robust analyses. One of the most widely used spectral
analysis and reduced-order modeling techniques is the dynamic mode
decomposition (DMD). DMD has been repeatedly applied to transonic buffet flows
[2, 3, 4, 5, 6, 7, 8]. However, DMD is known to suffer from sensitivities to
the selected settings and the data characteristics. In this work, we benchmark
various DMD variants and identify a robust methodology that yields physical
and repeatable results.

The DMD is a purely data-driven modal decomposition technique [9, 10]. In
contrast to similar techniques, such as proper orthogonal decomposition (POD),
DMD associates each spatial mode with exactly one frequency and growth rate.
The coherence of specific frequencies is often more informative for
understanding the flow dynamics than a ranking based on second-order flow
statistics [11], especially if acoustic waves need to be identified. Despite
its strengths, the DMD core algorithm is sensitive to noise and transients in
the data [12] and requires careful parameter tuning, e.g., for the number of
POD modes or regularization parameters. Here, we use the term noise
generically for components in the data that are unwanted, e.g., because they
have little to no coherence. In this context, turbulent fluctuations in scale-
resolved simulations are also considered noise. To improve DMD’s accuracy and
robustness, a variety of extensions to the core algorithm have been developed
[13]. In the following, we review various DMD implementations applied to
transonic shock buffet flows.

Differences between the DMD implementations may be grouped according to (i)
the type of data used, (ii) the definition of the state vector, and (iii) the
employed DMD algorithm and its settings. The type of data is either numerical
or experimental, which affects the amount of information available. While
simulations provide pressure, velocity, and other field variables across the
entire flow domain, experiments typically yield only a single surface (e.g.,
via instationary pressure sensitive paint, iPSP [6]), or field (e.g., via
particle image velocimetry, PIV [7]) variable. Moreover, simulations provide
typically higher temporal and spatial resolution, while experiments provide
much longer time sequences. Table 2 illustrates the typical differences in
selected studies between the type of data employed for DMD of a flow
experiencing transonic buffet. The ratio of sampling frequency $f_{s}$ and
buffet frequency $f_{b}$ yields the number of snapshots per buffet cycle,
whereas $N_{\mathrm{cycle}}$ is the overall number of recorded buffet cycles.

Table 2: Summary of selected literature resources illustrating the differences between the type of data employed for DMD of a flow experiencing transonic buffet. The ratio of sampling frequency $f_{s}$ and buffet frequency $f_{b}$ yields the number of snapshots per buffet cycle. $N_{\mathrm{cycle}}$ is the overall number of recorded buffet cycles. authors | $f_{s}/f_{b}$ | $N_{\mathrm{cycles}}$ | method
---|---|---|---
Ohmichi et al. [2] | 100 | 20 | 3D ZDES
Poplingher and Raveh [3] | 10-1000 | 8 | 2D URANS
Massini et al. [6] | 30 | 700 | iPSP
Hoffmann et al. [7] | 23 | 256 | PIV

The state vector may consist of volume, slice, or surface data of one or
multiple field variables. Simulations provide a larger degree of freedom in
terms of possible state vectors. However, more complex state vectors in
combination with high-resolution meshes can quickly lead to intractable memory
requirements. Ohmichi et al. [2] used density, velocity, and pressure in a
volume about a swept wing to define the state. The authors also tested reduced
state vectors containing only velocity or only pressure and reported a
comparable dominant buffet frequency among all cases. Other numerical works
[3, 4, 5] analyzed only the pressure field. To better approximate the integral
operation of the inner product and reduce the mesh dependency of the DMD
results, Ohmichi et al. [2] weighted the state vector with the square root of
the cell volume. Such weighting was not reported by other authors [3, 4, 5].

The final grouping of DMD applied to transonic buffet flows relates to the
type of algorithm employed and its settings. The first setting common to all
DMD variants is the number of POD modes on which to project the linear
operator. While Ohmichi et al. [2] truncated the POD basis after 51 modes
without a clear justification, Das et al. [4] chose the number of POD modes
based on the reconstruction error. The latter criterion may appear more
rigorous at first sight. However, a closer inspection reveals that the rank
selection problem is now simply shifted to satisfy a user-defined
reconstruction error. Other authors [3, 6, 7, 5] appear to have used the full
POD basis, which is predetermined by the number of available snapshots.
However, the larger the POD basis, the more difficult the mode selection
becomes. Moreover, truncation of the POD basis is also important to reduce
noise in the data. To filter out important modes, some authors [2, 7] resorted
to sparsity-promoting DMD (spDMD) introduced by Jovanović et al. [14] or to a
variant thereof [15]. Poplingher and Raveh [3] employed the integral criterion
by Kou and Zhang [16], whereas Massini et al. [6] weighted the amplitudes with
the eigenvalues as suggested by Tu et al. [17]. The disadvantages of the spDMD
are its elevated cost that stems from the iterative solution procedure and the
additional parameter tuning of the sparsity control parameter. It is worth
mentioning the recent multi-resolution DMD (mrDMD) implementation of Das et
al. [5], where DMD was sequentially applied to increasingly smaller time
windows, from which low-frequency dynamics are gradually removed [18].
However, little insight was gained from the additional effort.

A common objective among most studies applying DMD on transonic buffet flows
is their main focus on identifying the dominant buffet frequency and its
harmonics. Although important, the buffet mode is certainly not the only
relevant one in the transonic buffet flow. Among the above-mentioned
investigations, only Hoffmann et al. [7] computed dominant modes at
frequencies considerably higher than that of the buffet, which the authors
attribute to vortex shedding at the trailing edge. Interestingly, the authors
also detected waves with a frequency similar to that of the vortex shedding
along a streamline above the airfoil. The waves travel roughly at the speed of
sound from the trailing towards the leading edge. However, the acoustic waves
were detected using spatial-temporal correlation coefficients and not through
DMD. A direct interaction between vortex shedding, acoustic waves, and shock
motion by means of DMD was therefore not established, even though DMD is, in
principle, a suitable tool to investigate such connections.

The plethora of methods and our own analysis leading up to this work [19]
demonstrate the sensitivity of the dynamic mode decomposition to the state
number definition, the rank truncation, the sampling rate, the time window,
and the employed DMD variant. This sensitivity makes the physical
interpretation of results and the quantitative comparison between studies
extremely challenging. In this work, we closely examine the aforementioned
sensitivities, identify possible pitfalls, and provide best practices to
perform DMD on flows exhibiting transonic shock buffet. For the sake of
reproducibility, we publish the complete work including simulation cases,
datasets, processing pipelines, and visualizations in a dedicated
repository333https://github.com/FlowModelingControl/naca0012_shock_buffet.

The remainder of this article is structured as follows. In section 2 the
simulation setup and data acquisition are briefly described. In section 3, we
recast DMD as a modular framework that consists of five steps with various
options. In section 4, we compare the outcome of various DMD variants, derive
best practices, and apply the recommended procedure to a slice of data. Our
final assessment and recommendations are provided in section 5.

##  2 Numerical setup

The current investigations are based on numerical simulation data. Within the
project, multiple two- and three-dimensional simulations with different mesh
resolutions and settings were conducted and validated [19]. However, the
current DMD analysis relies only on one three-dimensional numerical simulation
setup. In this section, we briefly describe the numerical approach, the
validation process, and the data acquisition. For more details, the reader is
referred to Weiner and Semaan [19].

###  2.1 Simulation approach

The simulations are performed using the rhoCentralFoam solver from the
OpenFOAM-v2012444https://www.openfoam.com/ library. OpenFOAM implements a
finite volume approach for unstructured meshes with a cell-centered, co-
located variable arrangement. Face fluxes are reconstructed based on the two
neighboring cell-centered values of each face. rhoCentralFoam solves the
conservation equations for mass, momentum, and energy by operator splitting
[20]. In the predictor step, the convection of conserved quantities, i.e.,
mass density, momentum density, and energy density, are computed explicitly.
In the corrector step for momentum and energy, the diffusive terms are
discretized and solved implicitly for the primitive variables, i.e., energy
and velocity. The convective transport is computed based on the interpolation
scheme by Kurganov et al. [21]. To interpolate the cell-centered values of the
conserved variables to the faces, van Leer’s limiter [22] is employed. Linear
interpolation is used for all remaining face interpolations.

The solver supports generic RANS, LES, or hybrid RANS-LES approaches. Attempts
to predict transonic buffet with standard SA and $k$-$\omega$-SST RANS method
at a critical angle of attack were not successful [19]. The solution converged
to a steady state, which is a commonly-reported problem when employing RANS
[19, 23, 3]. Instead of modifying model constants or increasing the angle of
attack to trigger the buffet, we resort to hybrid SA-IDDES simulations. Our
SA-IDDES implementation is based on the work by Gritskevich et al. [24]. As
detailed in [19], the IDDES simulations are capable of capturing the flow
physics including the buffet phenomenon.

The airfoil geometry is a simple subsonic NACA-0012 with a blunt trailing edge
of height $2.5\times 10^{-3}c$, where $c$ denotes the chord length. The choice
for the popular NACA-0012 is motivated by our desire to use a public domain
geometry for which experimental data are available. The airfoil geometry along
with the computational mesh is presented in figure 1. The block-structured
C-type mesh is generated with the blockMesh utility. It consists of a total of
$N_{\mathrm{cells}}\approx 4.9\times 10^{6}$ cells, and a wall-normal
resolution of $y^{+}=\frac{yu_{\tau}}{\nu}=3.73$, where $u_{\tau}$ is the
friction velocity and $\nu$ is the kinematic viscosity. Considering the
airfoil’s nose as the origin, the lower, upper, and upstream boundaries are
located at a distance of $10c$, whereas the downstream boundary is $30c$ away.
The mesh has a spanwise width of $0.2c$, and most surface elements have an
aspect ratio close to unity.

(a) full mesh

(b) airfoil vicinity

(c) top view of the airfoil surface mesh

Figure 1: Various views of the computational mesh employed in the three-
dimensional simulation.

Since the convective transport step is solved explicitly, the Courant number
is limited to $\mathrm{Co}_{\mathrm{max}}=0.25$. For larger values in the
range $0.3<\mathrm{Co}<1$, the numerical solution is stable but yields
spurious pressure values about the leading edge. The dimensionless time step
employed in the simulations is $\Delta\tilde{t}_{\mathrm{sim}}\approx
10^{-5}$, where the dimensionless time is defined as:

$\tilde{t}=\frac{ta_{\infty}}{c}=\frac{tU_{\infty}}{c\textrm{Ma}_{\infty}}\ ,$
(1)

and $a_{\infty}$, $U_{\infty}$, and $\mathrm{Ma}_{\infty}$ are the freestream
speed of sound, velocity, and Mach number, respectively. All simulations are
conducted at $\mathrm{Re}=U_{\infty}c/\nu=10^{7}$ and
$\textrm{Ma}=U_{\infty}/a_{\infty}=0.75$.

###  2.2 Validation

This section provides an excerpt of the more extensive validation procedure
presented in Weiner and Semaan [19]. The validation mainly relies on
experimental data from McDevitt and Okuno [25]. Figure 2 shows a comparison
between experimental and numerical pressure coefficient $c_{p}$ distributions
at a pre-onset angle of attack of $\alpha=2^{\circ}$. In the figure,
$\tilde{x}=x/c$ is a non-dimensionalized coordinate. Similar normalization is
applied to the other coordinates yielding $\tilde{y}$ and $\tilde{z}$. As the
figure shows, the numerical and experimental distributions compare favorably.
We note the discrepancy in the shock location, which is predicted
approximately $0.05c$ upstream compared to the experiment.

Figure 2: Comparison between the numerical and experimental pressure
coefficient distributions of the NACA-0012 at $\alpha=2^{\circ}$ and
$\textrm{Ma}_{\infty}=0.75$. The experimental data is from McDevitt and Okuno
[25].

At buffet conditions of $\alpha=4^{\circ}$, a comparison of the pressure
coefficient distribution is not possible, because of the lack of reference
experimental data. Therefore, we rely on the buffet frequency to validate the
simulation results at these conditions. Our results yield a buffet frequency
of $\tilde{f}_{b}=\frac{2\pi cf}{U_{\infty}}\approx 0.65$ compared to
$\tilde{f}_{b}=0.47$ from the experiment. This discrepancy is within the
expected range of numerically-obtained results from the literature considering
that no additional tuning of turbulence model parameters was conducted [26,
19, 16].

###  2.3 Data acquisition

Because we are interested in the transonic buffet dynamics, all results
presented in section 4 are based on the simulation at $\alpha=4^{\circ}$.
Snapshots of density, velocity, pressure, and local Mach number are saved at
an interval of $\Delta\tilde{t}\approx 5.5\times 10^{-2}$. The time window
used for the analysis, excluding the initial numerical transient, spans over
$\Delta\tilde{T}_{\mathrm{win}}\approx 64$ corresponding to 5 full buffet
cycles. Consequently, about 230 snapshots per cycle are available. This very
high sampling rate is chosen with the intention to analyze the influence of
the sampling rate on the DMD spectrum.

To illustrate various scenarios, we analyze three different types of data: (i)
surface pressure, (ii) a slice from the $\tilde{x}$-$\tilde{y}$-plane at
$\tilde{z}=0$, and (iii) volume data. For the surface data, we only consider
the pressure on the upper surface, since tests including both upper and lower
surfaces did not yield any significant differences. For both the slice and
volume data, only a subdomain spanning $-1\leq\tilde{x}\leq 4$ and
$-1\leq\tilde{y}\leq 1$ is retained. This masking operation reduces the
overall data size and excludes unconverged downstream-convected structures
from the analysis. For brevity, we only present results from a data slice.
Results from the surface pressure and the volume data are provided in the
repository.

##  3 Dynamic mode decomposition

Since its inception, many DMD variants have been proposed to improve its
accuracy, robustness, efficiency, and applicability; see reference [13] for a
comprehensive review. In this study, we are interested in identifying a DMD
“recipe” that yields robust results when subjected to noise and changes in the
rank truncation. Navigating the plethora of existing DMD variants is a
daunting task. Therefore, we recast the dynamic mode decomposition as a
modular five-step framework comprising (a) the definition of data matrices,
(b) the optional de-noising, (c) the operator identification, (d) the
computation of the mode amplitudes, and (e) the selection of dynamically
meaningful modes. Figure 3 illustrates this framework, wherein multiple DMD
variants are denoted. In principle, any option in a step can be combined with
any option from another step. However, not all combinations are sensible. In
the current section, we focus on sensible and computationally tractable
combinations.

However, before we detail the new DMD framework, we briefly review the
singular value decomposition (SVD), which plays an integral part in every DMD
variant. The SVD is a frequently-used matrix factorization technique that
yields the following decomposition

$\mathbf{X}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\ast},$ (2)

where $\mathbf{U}$ and $\mathbf{V}$ are unitary matrices with orthonormal
columns, and $\mathbf{\Sigma}$ is a diagonal matrix (considering the economy
SVD). Here, ∗ denotes the complex conjugate transpose. One of the most useful
properties of SVD is that it provides an optimal low-rank approximation of a
matrix $\mathbf{X}$, which can be mathematically expressed as:

$\underset{\mathbf{X}_{r}}{\mathrm{argmin}}\left|\left|\mathbf{X}-\mathbf{X}_{r}\right|\right|_{F}=\mathbf{X}_{r}\mathbf{\Sigma}_{r}\mathbf{V}^{\ast}_{r}\quad\textbf{subject
to}\quad\mathrm{rank}(\mathbf{X}_{r})=r,$ (3)

where $||\cdot||_{F}$ is the Frobenius norm, and $r$ is the rank [27]: The
rank $r$ is a user-defined parameter (a hyperparameter). Selecting a suitable
value for $r$ is not always an easy task, especially for turbulent flows. To
facilitate the rank truncation in the presence of noise, Gavish and Donoho
[28] suggested a simple set of equations to explicitly compute the optimal
rank as:

$r_{\mathrm{opt}}=\underset{i}{\mathrm{argmin}}|\sigma_{i}-\kappa(\beta)\sigma_{\mathrm{med}}|,$
(4)

where $\sigma_{\mathrm{med}}$ is the median singular value, $\beta$ is the
ratio of rows to columns of $\mathbf{X}$, and the function $\kappa(\beta)$ is
given by:

$\kappa(\beta)=0.56\beta^{3}-0.95\beta^{2}+1.82\beta+1.43.$ (5)

$\mathbf{X}=\left[\begin{array}[]{cccc}|&|&&|\\\ \mathbf{x}_{0}&\mathbf{x}_{1}&...&\mathbf{x}_{N-2}\\\ |&|&&|\\\ \end{array}\right]\quad\mathbf{Y}=\left[\begin{array}[]{cccc}|&|&&|\\\ \mathbf{x}_{1}&\mathbf{x}_{2}&...&\mathbf{x}_{N-1}\\\ |&|&&|\\\ \end{array}\right]$Data matrices total least-squares [29] |  robust [30]
---|---
$\underset{\mathbf{Y}+\mathbf{\Delta Y}=\mathbf{A}\left(\mathbf{X}+\mathbf{\Delta X}\right)}{\mathrm{argmin}}\left|\left|\left[\mathbf{\Delta X},\mathbf{\Delta Y}\right]^{T}\right|\right|_{F}$ | $\underset{\mathbf{L}+\mathbf{S}=\mathbf{M}}{\mathrm{argmin}}\left(\left|\left|\mathbf{L}||_{\ast}+\lambda_{0}||\mathbf{S}\right|\right|_{1}\right)$
Pre-processing/de-noising standard [10] |  unitary [31] |  optimized [32]
---|---|---
$\underset{\mathbf{A}}{\mathrm{argmin}}\left|\left|\mathbf{Y}-\mathbf{AX}\right|\right|_{F}$ | $\underset{\mathbf{A}^{\ast}\mathbf{A}=\mathbf{I}}{\mathrm{argmin}}\left|\left|\mathbf{Y}-\mathbf{AX}\right|\right|_{F}$ | $\underset{\mathbf{\omega},\mathbf{\Phi}_{b}}{\mathrm{argmin}}\left|\left|\mathbf{X}-\mathbf{\Phi}_{b}\mathbf{V}_{\mathbf{\omega}}\right|\right|_{F}$
Operator standard [10] |  optimal [14] |  sparsity-promoting [14] |
---|---|---|---
$\mathbf{b}=\mathbf{\Phi}^{\dagger}\mathbf{x}_{0}$ | $\underset{\mathbf{b}}{\mathrm{argmin}}\left|\left|\mathbf{X}-\mathbf{\Phi}\mathbf{D}_{\mathbf{b}}\mathbf{V}_{\mathbf{\omega}}\right|\right|_{F}$ | $\underset{\mathbf{b}}{\mathrm{argmin}}\left(\left|\left|\mathbf{X}-\mathbf{\Phi}\mathbf{D}_{\mathbf{b}}\mathbf{V}_{\mathbf{\omega}}\right|\right|_{F}+\gamma_{0}\left|\left|\mathbf{b}\right|\right|_{1}\right)$ |
Mode amplitudes amplitude [10] |  $\lambda$-weighted [17] |  integral contr. [16]
---|---|---
$I_{i}=|b_{i}|$ | $I_{i}=|b_{i}\lambda_{i}^{N-1}|$ | $I_{i}=\sum\limits_{j=0}^{N-1}|b_{i}\lambda_{i}^{j}|$
Mode selection$(a)$$(b)$$(c)$$(d)$$(e)$ Figure 3: Schematic illustration of
the DMD modular framework consisting of five distinct steps. Each step
includes various choices, which lead to different DMD variants. We note that
step (b) is optional and is therefore marked with a dashed line.

###  3.1 Data matrices

The first step in the modular DMD framework is the construction of two data
matrices $\mathbf{X}$ and $\mathbf{Y}$ comprising snapshots of the system’s
state at fixed time intervals. The snapshots in $\mathbf{Y}$ are shifted by
one timestep with respect to $\mathbf{X}$. Because SVD and DMD are data-driven
methods, one could theoretically use any flow variable or a combination of
multiple variables in the domain (or a subdomain) to define the state vector
$\mathbf{x}$ [13]. However, it is typically advantageous to use state vectors
that yield an interpretable and robust modal decomposition. In compressible
flows, the state vector could consist of density $\rho$, pressure $p$,
velocity $\mathbf{u}=\left(u_{x},u_{y},u_{z}\right)^{T}$, or other derived
field variables. It is implicitly assumed that all flow states are functions
of space. The inner product of two states $\mathbf{x}_{1}$ and
$\mathbf{x}_{2}$ is then defined as:

$\langle\mathbf{x}_{1},\mathbf{x}_{2}\rangle=\int_{\Omega}\mathbf{x}_{1}\cdot\mathbf{x}_{2}\>\mathrm{d}v,$
(6)

where $\Omega$ is the spatial domain under consideration. The resultant norm
is:

$||\mathbf{x}||^{2}=\int_{\Omega}\mathbf{x}\cdot\mathbf{x}\>\mathrm{d}v.$ (7)

The norm plays an essential role in the operator identification step. For
example, if we choose $\mathbf{x}=\left(u_{x},u_{y},u_{z}\right)^{T}$, the
squared norm yields twice the flow’s kinetic energy in $\Omega$. If the flow
is incompressible, the kinetic energy is a preserved quantity and we could
enforce this constraint in the DMD operator identification [31].

For compressible flows, where the state variables are coupled, Rowley et al.
[33] suggested a state vector that combines the velocity vector and the speed
of sound as:

$\mathbf{x}=\left(\sqrt{\frac{2}{(\gamma(\gamma-1))}}\>a,u_{x},u_{y},u_{z}\right)^{T},$
(8)

where $\gamma$ is the adiabatic index and $a$ is the local speed of sound. The
squared norm of definition (8) yields twice the flow’s stagnation energy [33],
which is expected to be approximately conserved in the type of flows
considered here. More details on the derivation of the compressible inner
product are provided in appendix Details about the compressible inner product.

When dealing with snapshot data originating from a finite-volume solution, the
state vector is finite-dimensional, and the inner product 6 is approximated by
a sum over cell-centered values multiplied by the cell volume. If the mesh is
uniform, the volume-weighting can be neglected. However, in the present study,
the mesh around the airfoil is highly non-uniform, and neglecting the cell
volume would introduce a strong mesh dependency in the DMD results. The
weighting is implemented in practice by multiplying each entry of the state
vector by the square root of the corresponding cell volume. The dot-product of
two weighted state vectors then yields the desired approximation of the volume
integral.

Besides the physically-motivated construction of the state vector, it is also
possible to enrich the state with time delays by stacking multiple consecutive
snapshots into extended column vectors [34, 35]. However, each additional
snapshot increases the size of the final data matrices, which can be
prohibitive for large datasets. The problem can be mitigated significantly by
a dimensionality reduction step preceding the time delay embedding [35].
Nonetheless, we do not consider time delay embeddings in this study.

###  3.2 Pre-processing/de-noising

The second optional step addresses the often-encountered small-scale broadband
dynamics that stem from turbulence. To mitigate the sensitivity of DMD
analysis to the resolved turbulence, it might be sensible to apply a de-
noising technique to the data before proceeding with the subsequent steps. In
the total least-squares DMD (TDMD) [29], the data are considered to have some
noise in $\mathbf{X}$ and $\mathbf{Y}$, denoted by $\Delta\mathbf{X}$ and
$\Delta\mathbf{Y}$, which one attempts to minimize as,

$\underset{\mathbf{Y}+\mathbf{\Delta
Y}=\mathbf{A}\left(\mathbf{X}+\mathbf{\Delta
X}\right)}{\mathrm{argmin}}\left|\left|\left[\mathbf{\Delta X},\mathbf{\Delta
Y}\right]^{T}\right|\right|_{F}\,.$ (9)

While the linear operator $\mathbf{A}$ appears in the TDMD constraint, Hemati
et al. [29] present the solution of the total least-squares problem as a two-
step procedure involving an initial subspace projection step followed by the
operator identification. We adopt this viewpoint and consider the subspace
projection as an optional pre-processing step before the operator
identification.

An alternative approach for de-noising applies a robust SVD on the data matrix
decomposing it into low-rank $\mathbf{L}$ and sparse contributions
$\mathbf{S}$ by solving the following optimization problem [30]:

$\underset{\mathbf{L}+\mathbf{S}=\mathbf{M}}{\mathrm{argmin}}\left(\left|\left|\mathbf{L}||_{\ast}+\lambda_{0}||\mathbf{S}\right|\right|_{1}\right).$
(10)

The operator identification is subsequently performed on the low-rank
contributions $\mathbf{L}$. While the TDMD approach requires only explicit
computations, the robust SVD method involves a relatively costly iterative
solution [30], which is why we exclude it from our current assessment study.

###  3.3 Operator identification

The third step identifies a linear operator $\mathbf{A}$ that maps the state
$\mathbf{x}_{n}$ to $\mathbf{x}_{n+1}$, such that:

$\mathbf{x}_{n+1}=\mathbf{A}\mathbf{x}_{n}.$ (11)

Since the exact linear operator will typically vary between snapshot pairs,
$\mathbf{A}$ is defined in the least-squares sense as:

$\underset{\mathbf{A}}{\mathrm{argmin}}\left|\left|\mathbf{Y}-\mathbf{AX}\right|\right|_{F}.$
(12)

The solution of problem (12) without any further constraints on $\mathbf{A}$
yields the most common operator definition:

$\mathbf{A}\coloneqq\mathbf{Y}\mathbf{X}^{\dagger},$ (13)

where $\dagger$ denotes the Moore-Penrose inverse.

The operator’s eigen-decomposition, namely
$\mathbf{A}=\mathbf{\Phi\Lambda\Phi}^{-1}$, has several useful properties. For
example, it can be used to make future predictions based on the relation:

$\mathbf{x}_{n}=\mathbf{\Phi\Lambda}^{n}\mathbf{b},$ (14)

where $\mathbf{\Phi}$ denotes the DMD modes and
$\mathbf{b}=\mathbf{\Phi}^{-1}\mathbf{x}_{0}$ the mode amplitudes [10]. In
fluid flow problems, the state vector is typically extremely large such that
the eigendecomposition of $\mathbf{A}$ has to be approximated. The most widely
used procedure invokes a low-rank approximation of $\mathbf{A}$ by means of an
SVD of the first snapshots matrix,
$\mathbf{X}\approx\mathbf{U}_{r}\mathbf{\Sigma}_{r}\mathbf{V}_{r}^{\ast}$ and
the definition of a reduced operator

$\tilde{\mathbf{A}}=\mathbf{U}_{r}^{\ast}\mathbf{AU}_{r}=\mathbf{U}_{r}^{\ast}\mathbf{X}^{\prime}\mathbf{V}_{r}\mathbf{\Sigma}_{r}^{-1}\,,$
(15)

which has the same non-zero eigenvalues as $\mathbf{A}$. The modes are then
reconstructed following the formulation suggested by Tu et al. [17] yielding
the so-called exact DMD modes:

$\mathbf{\Phi}=\mathbf{Y}\mathbf{V}_{r}\mathbf{\Sigma}_{r}^{-1}\mathbf{W},$
(16)

where the matrix $\mathbf{W}$ contains the eigenvectors of
$\tilde{\mathbf{A}}$.

This standard operator identification is the most widely used one.
Nonetheless, it still lacks accuracy and robustness for a range of
applications. To address this shortcoming, many variants have been proposed.
One variant is the optimized DMD [32], where the problem is cast as an
optimization problem in terms of the scaled eigenvector matrix
$\mathbf{\Phi}_{b}$ and the Vandermode matrix $\mathbf{V}_{\mathbf{\omega}}$

$\underset{\mathbf{\omega},\mathbf{\Phi}_{b}}{\mathrm{argmin}}\left|\left|\mathbf{X}-\mathbf{\Phi}_{b}\mathbf{V}_{\mathbf{\omega}}\right|\right|_{F}.$
(17)

The Vadermode matrix is based on the trace of the eigenvalue matrix, i.e.,
$\mathbf{\omega}=\mathrm{tr}\left(\mathbf{\Omega}\right)$; see appendix Time-
continuous linear operator for more details. The reformulated optimization
problem comes at the cost of a much more involved iterative solution procedure
[32]. Even though the optimized DMD is usually more accurate [32, 36], the
number of potentially important modes may increase, since the symmetry of
complex-conjugate mode pairs is broken [36]. Moreover, simple numerical
experiments show that the resulting eigenvalues may be very different from
that of the original dynamical system [12]. Due to the computational and
numerical complexity of the optimized DMD, we do not consider this variant in
our investigations.

One variant we do consider is from the physics-informed DMD family proposed by
Baddoo et al. [31], where a number of physically-motivated constraints for the
operator are introduced. One constraint we adopt in this study is self-
adjointness $\mathbf{A}^{\ast}\mathbf{A}=\mathbf{I}$, which assumes
conservation of the state vector’s norm. It is mathematically expressed as:

$\underset{\mathbf{A}^{\ast}\mathbf{A}=\mathbf{I}}{\mathrm{argmin}}\left|\left|\mathbf{Y}-\mathbf{AX}\right|\right|_{F}.$
(18)

This constraint is only sensible if the state’s norm is actually conserved (at
least approximately). We examine this assumption in the following.

Figure 4 presents the temporal evolution of various state vector norms
normalized by their corresponding temporal mean. The state vector norms for
field variables are presented with and without cell volume weighting. Figure 4
(a) shows the state vector norm time sequence that is based on the surface
pressure coefficient $c_{p}$ on the suction side. As the figure shows, the
norm undergoes large fluctuations of more than $10\%$ caused by the shock
motion and the boundary layer separation. These large fluctuations clearly
indicate that the norm of this state vector is not conserved sufficiently.

In contrast to the surface pressure, field variable norms fluctuate less. The
density-based state vector norms presented in figure 4 (b) for the unweighted
and weighted state vary $\approx 4.5\%$ and $\approx 1.3\%$, respectively. The
strong reduction resulting from weighting demonstrates its importance in
reducing mesh dependency.

The state vector norm time sequences based on the two- and three-dimensional
velocity (figure 4 (c)) yield normalized kinetic energy fluctuations of less
than $\approx 1.3\%$. Differences between the two- and three-dimensional state
vector norms are minimal since the spanwise velocity is small compared to the
other two components and thus contributes little to the flow’s kinetic energy.
Finally, we present the normalized state vector norms based on the local speed
of sound solely, as well as in combination with the velocity vector in figure
4 (d). The latter definition corresponds to the state vector in equation (8).
This physically-motivated state vector fluctuates less than $\approx 0.3\%$.

In summary, state conservation may be assumed for all herein-considered
definitions except for the surface pressure. Particularly suited is the state
vector that conserves the stagnation energy (equation (8)). This assessment is
in agreement with the results presented in section 4.

Figure 4: Time sequence of various state vector norms $||\mathbf{x}||$
normalized with their corresponding temporal mean norm $||\mu_{\mathbf{x}}||$.
The legends indicate the field variables comprising the state vector. The
keyword weighted indicates a weighting with the square root of the cell
volume. Except for the surface pressure distributions $c_{p}$, all field data
are extracted from a planar slice; see section 2.

###  3.4 Mode amplitudes

The fourth step concerns the computation of the mode amplitude vector
$\mathbf{b}$. The standard relation to computing the mode amplitudes results
from the solution of the initial value problem (29) [10]:

$\mathbf{b}=\mathbf{\Phi}^{\dagger}\mathbf{x}_{0}\ .$ (19)

However, since $\mathbf{A}$ is defined in the least-squares sense, there might
be better options to compute the amplitudes, such as with optimal DMD [14],
where $\mathbf{b}$ is optimized for a given eigendecomposition:

$\underset{\mathbf{b}}{\mathrm{argmin}}\left|\left|\mathbf{X}-\mathbf{\Phi}\mathbf{D}_{\mathbf{b}}\mathbf{V}_{\mathbf{\omega}}\right|\right|_{F}\
.$ (20)

Another approach, the sparsity-promoting DMD (spDMD) [14], adds an $L_{1}$
norm to the optimization problem:

$\underset{\mathbf{b}}{\mathrm{argmin}}\left(\left|\left|\mathbf{X}-\mathbf{\Phi}\mathbf{D}_{\mathbf{b}}\mathbf{V}_{\mathbf{\omega}}\right|\right|_{F}+\gamma_{0}\left|\left|\mathbf{b}\right|\right|_{1}\right)\,.$
(21)

The additional complexity of the spDMD requires an iterative solution
procedure. We performed a priori tests with spDMD and found no advantage over
optimal DMD when combined with the improved selection criteria presented in
the next section, which is in agreement with the observations in [16].
Moreover, the spDMD introduces an additional hyperparameter $\gamma_{0}$ to
control the sparsity, which requires tuning.

###  3.5 Mode selection

The fifth and final step in the framework addresses the mode selection. Since
the DMD modes are of unit length and the eigenvalues reflect the dynamic
behavior (growing, decaying, periodicity), the magnitude of the amplitudes
$b_{i}$:

$I_{i}=|b_{i}|\,.$ (22)

may be used to sort the DMD modes. However, noise in the data often yields
modes associated with large amplitudes and $|\lambda_{i}|<1$, leading to
quickly-vanishing modes contributing very little to the full reconstruction.
To account for the dampening, several authors [26, 6] rank the modes based on
$\lambda_{i}^{N-1}$-weighted amplitudes as:

$I_{i}=|b_{i}\lambda_{i}^{N-1}|\,.$ (23)

A more refined version of this idea was introduced by Kou and Zang [16], where
the eigenvalue-weighted amplitudes are summed up over all $N$ time steps:

$I_{i}=\sum\limits_{j=0}^{N-1}|b_{i}\lambda_{i}^{j}|\,.$ (24)

We refer to this criterion as the integral mode selection.

###  3.6 Abbreviations of DMD variants

The various options in each step create a vast space of possible combinations.
To distinguish between the DMD variants we have tested, we use the
abbreviations summarized in table 3. For example, “opt. TDMD, int.” refers to
the DMD variant with total least-squares pre-processing, the standard operator
identification, optimal mode coefficients, and mode ranking based on the
integral mode selection criterion. We emphasize that we restrict our tests to
sensible and computationally-tractable combinations.

Table 3: List of abbreviations of DMD options according to the framework illustrated in figure 3. Abbreviation | Explanation
---|---
DMD | standard operator definition (12)
TDMD | total least-squares pre-processing (9)
UDMD | unitary linear operator (18)
opt. | optimal mode amplitudes (20)
amp. | mode selection based on $|b_{i}|$ (22)
int. | mode selection based on the integral criterion (24)

##  4 Results

In this section, we test different DMD variants with different settings on the
transonic buffet flow. We assess each algorithm based on its robustness and
its accuracy in identifying the flow dynamics. Robustness tests include rank
sensitivity, sampling rate sensitivity, and time-window sensitivity. Using the
test results, we identify a robust recipe and apply it to a planar slice
through the domain. If not indicated otherwise, the state vector in all tests
is based on the stagnation energy. We emphasize that we performed the same
analyses on surface and volume data, but these results are omitted here for
brevity. The interested reader is referred to the accompanying repository. All
DMD analyses are conducted with the Python library flowTorch [37].

Before we proceed with the tests, it is important to identify the relevant
reference frequencies against which we benchmark our DMD results. For this
purpose, we analyze the flow velocity magnitude $U=||\mathbf{u}||$ at four
characteristic probe locations. The exact probe locations are depicted in
figure 5, which presents the normalized flow velocity magnitude’s mean
$\mu_{U}/U_{\infty}$ and standard deviation $\sigma_{U}/U_{\infty}$. The four
sensors are positioned both heuristically and based on the standard deviation
distribution.

Figure 5: Temporal mean and standard deviation of the flow velocity magnitude
$U=||\mathbf{u}||$ normalized with $U_{\infty}$. Markers indicate the probe
locations.

Figure 6 presents the power spectral density (PSD) of the four sensor signals.
All four spectra exhibit a peak at $\tilde{f}_{b}=\frac{2\pi
cf}{U_{\infty}}\approx 0.65$, which is the buffet frequency. The prominence of
the buffet frequency is a clear indicator of the shock buffet dynamics
dominance across the flow field. Besides the peak at $\tilde{f}_{b}\approx
0.65$, three sensor signals show harmonics of $\tilde{f}_{b}$. Also of
relevance are the dynamics of the wake shedding, which are best captured with
sensor 2 located near the airfoil’s trailing edge. The power spectral density
shows multiple peaks associated with the wake shedding dynamics within the
range $8\tilde{f}_{b}<\tilde{f}<25\tilde{f}_{b}$.

Figure 6: Power spectral density (PSD) of the velocity magnitude
$U=||\mathbf{u}||$ at four selected probe locations. The index $i$ of probe
$U_{i}$ corresponds to the marker label in figure 5. Also illustrated are
frequencies of interest presented in nondimensional form.

###  4.1 Rank sensitivity

In the first of a series of tests, we assess the capability of different DMD
variants to robustly identify the buffet frequency under rank variation.
Figure 7 shows close-ups of the stability circle for twelve DMD variants. Each
subfigure presents the dominant eigenvalue for a rank variation spanning
$5<r<500$. The red star indicates the mean eigenvalue over the tested rank
range, whereas the red ellipse quantifies the extent of the scatter. The
ellipse is oriented along the two principal components of the eigenvalues,
with the half-axes spanning twice the standard deviation along each axis.
Subfigures with a non-visible ellipse indicate a scatter beyond the shown
range.

For each rank setting, the dominant eigenvalue is identified following two
mode selection criteria. Hence, the spectral contribution $I_{i}$ of
eigenvalue $i$ is either the largest magnitude of amplitude $b_{i}$ (amp.) or
of the integral criterion (int.). We reasonably assume that the most dominant
eigenvalue with a positive imaginary part is that of the buffet. We recall
that $r$ sets the dimensionality of the POD basis. As $r$ increases, smaller
flow structures or potential noise are included, which, ideally, should not
affect the dominant buffet eigenvalue. Hence, the markers in figure 7 should
cluster in a small region as the rank is varied. However, as the figure shows,
many variants exhibit a scatter. The scatter is reduced, i.e., the robustness
is increased, when the mode selection based on the integral criterion is
employed rather than that based on the amplitude. Another significant
robustification is achieved when utilizing optimal mode amplitudes. These
observations hold for all other state vectors tested in this work.

Figure 7: Most dominant eigenvalues according to $I_{i}$ for the stagnation
energy state vector over a rank range of $5<r<500$. The red star indicates the
mean eigenvalue over the rank range. The red ellipse is oriented along the two
principal components of the eigenvalues, and the half-axes span twice the
standard deviation along each axis.

Both DMD and TDMD perform similarly, but the scatter within the dominant
cluster is slightly larger for the TDMD. In comparison, UDMD delivers
perfectly consistent results, especially when combined with optimal
amplitudes. The statement holds for the other tested state vectors. Only for
the surface pressure, UDMD is not able to identify the buffet eigenvalue
consistently, which is not surprising considering the state conservation
analysis presented in section 3.3.

The above-mentioned observations are also mirrored in figure 8, which presents
the frequency associated with the largest eigenvalue amplitude
$\tilde{f}_{\mathrm{max}}$ against the rank for various DMD variants. To
improve visibility, markers below $0.9\tilde{f}_{\mathrm{max}}$ and above
$1.1\bar{f}_{\mathrm{max}}$ are pinned to the plot’s lower and upper boundary.
$r_{\mathrm{opt}}$ denotes the optimal rank according to equation (4). We
utilize the buffet frequency resulting from the sensor spectral analysis as a
reference. In addition to the previous observations, figure 8 provides the
first quantitative description of the rank sensitivity. Both DMD and TDMD are
rank-sensitive, but TDMD yields more outliers and more scatter about the
reference frequency. Again, the optimal mode amplitude and integral mode
selection approaches exhibit the most robust and accurate results. Among the
three operator variants, the standard DMD matches the buffet frequency best if
the POD basis is sufficiently truncated, e.g, $r<400$. In both DMD and TDMD
variants, the optimal rank $r_{\mathrm{opt}}$ according to equation (3)
appears to provide a good compromise between robustness and flow resolution.
UDMD is the least rank-sensitive algorithm, regardless of the mode amplitude
and mode selection approach. However, its $\tilde{f}_{\mathrm{max}}$ is
slightly lower than the reference $\tilde{f}_{b}$. It is not entirely clear
why UDMD delivers a frequency shift of the dominant mode.

Figure 8: Frequency associated with the largest eigenvalue amplitude
$\tilde{f}_{\mathrm{max}}$ depicted against the rank for various DMD variants.
To improve visibility, markers below $0.9\tilde{f}_{\mathrm{max}}$ and above
$1.1\bar{f}_{\mathrm{max}}$ are pinned to the plot’s lower and upper boundary.
$r_{\mathrm{opt}}$ denotes the optimal rank according to equation (3).

The strong robustness to rank truncation, and the ability to identify the
buffet mode frequency might give the impression of an overall superiority of
UDMD. However, in our analysis, we have observed strange artifacts appearing
in the buffet mode distributions (as well as other modes) when the rank is
changed. On the other hand, the DMD buffet mode remains almost identical over
a wide range of rank truncations. Again, it is not entirely clear why UDMD
exhibits these artifacts. We refer the reader to appendix UDMD buffet mode for
a comparison between the DMD and UDMD mode distributions. For these reasons,
we mostly employ henceforth the standard DMD variant for operator selection.

The rank sensitivity investigation is not complete without examining the full
spectrum. Figure 9 shows full DMD spectra for ranks $5<r<250$ using the
optimal DMD variant with integral mode selection. The markers’ size is scaled
by the integral criterion, and the ten most dominant modes are colored red.
The probability density (PD) to the right of the scatter plot visualizes the
regularity of occurrence of each frequency. Ideally, the red markers follow a
straight horizontal line, which corresponds to a small rank sensitivity of the
dominant modes. As can be seen, this robust DMD variant consistently
identifies not only the buffet frequency but also several of its harmonics
over a large range of ranks. For higher frequencies, the variability increases
with the rank but remains acceptable. It is worth noting the identified modes
at $\tilde{f}\approx 9$ and $\tilde{f}\approx 12$, which we shall revisit
later in the manuscript.

Figure 9: DMD spectra over the rank range $5<r<250$ for the state vector based
on the stagnation energy. The marker size is scaled by $I_{i}$ (integral
criterion) and the top 10 markers are colored in red.

The rank sensitivity is also influenced by the choice of the state vector.
Figure 10 shows DMD spectra for the density-based state vector. Compared to
figure 9, the rank sensitivity is much larger, especially for frequencies
$\tilde{f}>2$. Moreover, no relevant modes in the vortex shedding regime about
$\tilde{f}\approx 10$ are detected. Thus, physically-motivated state vectors
are better suited for the DMD analysis, even if the unitary DMD variant is not
employed.

Figure 10: DMD spectra for $5<r<250$ and the density-based state vector. The
marker size is scaled by $I_{i}$ (integral criterion) and the top 10 markers
are colored in red.

###  4.2 Sampling rate sensitivity

In the second series of tests, we investigate the DMD sensitivity to the
sampling rate $\tilde{f}_{s}$. Here, three DMD variants are examined, and the
relevant mode frequencies are identified according to the integral criterion.
We gradually decrease the sampling rate from the maximum value
$\tilde{f}_{s}/\tilde{f}_{b}\approx 230$, halving the number of snapshots at
each step. The rank is selected automatically following Gavish and Donoho
[28].

Figure 11: Ten most relevant DMD mode frequencies for a range of sampling
rates. The dashed line demarcates the Nyquist frequency. All three DMD
variants employ the integral criterion.

Figure 11 shows the identified ten most relevant mode frequencies for a range
of sampling rates. The dashed line demarcates the Nyquist frequency. As the
figure shows, all DMD variants identify the buffet frequency reliably, even if
as few as three snapshots per buffet cycle are used. However, for sampling
rates below $\tilde{f}_{s}/\tilde{f}_{b}\approx 100$, the upper end of the
spectrum gets clipped. The same threshold value was also reported by
Poplingher and Raveh [3] in their DMD analysis. This clipping is associated
with aliasing, where high-frequency modes associated with wake phenomena, are
now aliased towards lower frequencies. This is observed by the frequency
clustering near $f_{\mathrm{Nyquist}}$ and by the high frequencies modes now
appearing at lower frequencies. These observations are similar across all the
tested state vectors and DMD variants. We note the difficulty in avoiding such
aliasing in experimental settings. Considering the values listed in table 2
and the current state-of-the-art field measurement techniques, the current
achievable sampling rates are barely enough to capture the shedding frequency.

###  4.3 Time window sensitivity

As a final test, we examine DMD’s sensitivity to the time window selection,
which consists of sensitivity to the starting time and the time window length.
These tests are particularly relevant to numerical simulations, where the user
has to decide when and how long to sample.

For the starting time analysis, we apply DMD on a series of sequentially-
shifted sequences, each sequence being shifted by a small time step. We have
tested three different time window lengths spanning between two and four
buffet cycles. In figure 12, we present the results for three buffet cycles
long sequence. The spectra are visualized analogously to the ones in section
4.1, but with the abscissa now showing the starting time $\tilde{t}_{0}$
instead of the rank. All spectra are computed on a planar slice through the
flow field using optimal DMD. The marker size in the figure is scaled by
$I_{i}$ using the integral criterion, and the top ten markers are colored in
red.

Figure 12: Opt. DMD spectra for a shifting time window of constant size
covering three buffet cycles. The marker size is scaled by $I_{i}$ (integral
criterion) and the top ten markers are colored in red. The state vector is the
stagnation energy-conserving definition in equation (8).

As figure 12 shows, the starting time has little impact on the DMD spectrum.
Most relevant frequencies, including the high frequencies, remain consistent
during the time shift. We note that a perfectly non-varying spectrum with
constant top frequencies over a shifting $\tilde{t}_{0}$ is not expected,
because the IDDES simulation at this Reynolds number yields buffet cycle-to-
cycle variations. On the other hand, more consistent spectra can be expected
with URANS data, because of its more regular time sequences.

For an appropriate rank selection and state vector definition, such as the one
used in figure 12, the DMD results are also insensitive to the time window
length. The low sensitivity to the total sampling time is in agreement with
previously-reported DMD results from numerical and experimental investigations
[3, 2, 6]. DMD spectra using other state vector definitions and rank
truncations exhibit more variations (not shown) than what we observe in figure
12. More details are provided in the supplementary repository.

In summary, for an appropriate rank and state vector definition, the choice of
a time window has little influence on the DMD results. We note that longer
time sequences might be necessary if the cycle-to-cycle differences become
stronger, which might be the case in experiments.

###  4.4 Application

Based on the test results, we apply DMD on a planar slice through the flow
field employing the following settings:

  * •

State vector: the state vector is the stagnation energy-conserving definition
(8), excluding $u_{z}$ because of its small contribution. The cell-centered
values are weighted by the square root of the cell volume. All available
snapshots are used, allowing to detect frequencies up to
$\tilde{f}_{max}\approx 76$.

  * •

Operator: the standard operator definition (13) is employed without unitary
constraint.

  * •

Rank:  the size of the SVD basis is selected according to equation (3).

  * •

Mode amplitudes: the optimization problem (20) is solved to obtain the optimal
mode amplitudes.

  * •

Mode selection: the modes’ importance is evaluated based on the integral
criterion (24).

Similarly successful implementations on surface and volumetric data using the
same settings are provided in the supplementary repository.

Figure 13 presents the resulting DMD spectrum that exhibits a clear peak at
$\tilde{f}_{b}\approx 0.63$, which is in agreement with the frequency
extracted from the lift coefficient and the velocity probes
$\tilde{f}_{b}\approx 0.65$. Besides the main buffet mode, the spectrum shows
at least two harmonics to $\tilde{f}_{b}$, as well as modes within the
frequency range $12<\tilde{f}<13$ that are associated with vortex shedding. No
dominant modes were detected beyond the vortex shedding frequency.

Figure 13: DMD spectrum. The dashed line denotes the buffet frequency
$\tilde{f}_{b}\approx 0.65$ that was extracted from the power spectral density
of the velocity probes. For better visibility, the spectrum is limited to
frequencies less than $\tilde{f}=25$.

Figure 14 presents the real part of two selected modes: the buffet mode, and a
wake mode. The $u_{x}$-component of the buffet mode shows activity about the
shock, within the separated boundary layer, and in the rear part above the
boundary layer, which is in agreement with Das et al. [5]. Downstream of the
shock, the $u_{y}$-component exhibits a region of high values where the
$u_{x}$-component has a minimum, indicating a large periodic displacement of
fluid normal to the airfoil’s surface. The $a$-component of the buffet mode
appears mostly like a mixture of the other two components.

Figure 14: Selected DMD modes associated with shock buffet (top row) and wake
phenomena (bottom row). The figures only show the real part.

The second mode at $\tilde{f}=12.58$ captures flow separation in the boundary
layer and the wake. The structures with alternating intensity correspond to
the expected vortex shedding. Moreover, vertical lamellar structures below and
above the airfoil can be observed, particularly in the $u_{x}$\- and
$a$-components. The nature of these waves is elucidated through the
examination of the phase angle and the wave propagation speed. Figure 15
presents the phase angle distribution, which is computed as

$\phi=\tan^{-1}\frac{\Im(u_{x})}{\Re(u_{x})}\,$

where $\Im(u_{x})$ and $\Re(u_{x})$ denote the imaginary and real components
of the mode. The white contour lines demark the phase at $\phi=0$. The waves’
propagation direction can be inferred from the phase angle gradient [3].
Outside the wake and above and below the airfoil, the slope in the $\tilde{x}$
direction is positive, indicating that the waves travel upstream. On the other
hand, the phase angle slope in $\tilde{x}$ within the separated wake is
negative, indicating wave propagation downstream.

Figure 15: Phase angle distribution of a shedding mode at $\tilde{f}=12.58$.
The white contour lines demark zero phase angle.

Besides the propagation direction, the phase angle gradient can also be used
to compute the wave propagation speed [3]. However, in this study, we opt for
using the spatio-temporal correlation coefficient [7] to infer the wave speed
along a line. Figure 16 presents the spatio-temporal correlation coefficient
of the velocity magnitude from DMD-reconstructed snapshots using a single mode
at $\tilde{f}=12.58$. The correlation is evaluated along $0.7\leq\tilde{y}\leq
0.8$ at $\tilde{x}=0.5$. The propagation velocity can be elucidated from the
slope of the ridge lines [7]. As expected, the velocity is negative (i.e.,
pointing upstream) and very close to the speed of sound in the absolute frame
of reference (the free stream velocity is set such that $M_{\infty}=0.75$).
Hence, these waves at $\tilde{f}=12.58$ are the well-known acoustic waves that
originate in the shear layer and the trailing edge and propagate upstream [38,
1]. They are clearly slaved to the shear layer dynamics.

Figure 16: Spatio-temporal correlation coefficient of the velocity magnitude
from DMD reconstructed snapshots using a single mode at $\tilde{f}=12.58$. The
correlation is evaluated along $0.7\leq\tilde{y}\leq 0.8$ at $\tilde{x}=0.5$.
The dashed line demarks the wave propagation speed.

Before we conclude this section, we would like to point out the simplicity
involved in computing the wave propagation speed. While Feldhusen-Hoffmann et
al. [7] applied a band-pass filter to the PIV data before computing the
spatio-temporal correlation, such filtering is not necessary for a
reconstructed DMD mode, since it is associated with a single frequency.
Omitting the filtering step is a clear advantage since no prior knowledge of
the wave frequency range is required.

##  5 Conclusions and recommendations

In this study, we benchmark various DMD variants on simulation data of
transonic shock buffet flow with the objective to identify a robust and
accurate methodology. The analysis is enabled by the introduction of a simple
modular framework of the DMD algorithm. The framework consists of five
distinct steps, where each step includes various choices that lead to
different DMD variants.

The analysis is conducted on a scale-resolving IDDES simulation of a NACA-0012
under transonic buffet conditions. The IDDES simulation with the Spalart-
Almaras turbulence model is performed using rhoCentralFoam solver from
OpenFOAM-v2012. The simulation is validated against the experimental results
of McDevitt and Okuna [25].

The conducted tests examine the various DMD variants’ sensitivity to rank
truncation, sampling rate, and time window selection. The tests yield a DMD
variant with corresponding settings that provide robust and accurate
identification of the relevant flow dynamics. The suggested workflow can
detect the shock-boundary-layer interaction reliably and also identifies
upstream-propagating acoustic waves that originate in the shear layer and at
the trailing edge. The waves’ frequency is directly linked to the vortex
shedding at the trailing edge. Similar results in terms of the DMD spectra and
flow phenomena encoded in the modes could be obtained based on slice, surface,
and volume data, which is in agreement with the observation by Ohmichi et al.
[2] and supports the possibility to compare numerical and experimental data
from different sources.

Another objective of this study is to provide recommendations on how to best
conduct dynamic mode decomposition on a transonic buffet flow. Based on the
test results, the following general recommendations are inferred:

  1. 1.

When possible, employ a physically-motivated state vector, such as that
defined in equation (8). This state vector is physically interpretable and
enables using the so-called unitary DMD. Clearly, this recommendation is
challenging for experimental data, where only velocity or pressure fields are
typically acquired.

  2. 2.

Weigh the cell-centered values by the square root of the cell volume (or
surface area in case of two dimensions). The weighing improves the inner
product operator calculation, reduces mesh dependency, and yields more
accurate results.

  3. 3.

Extra care should be taken when truncating the SVD, which has shown the
highest sensitivity on the DMD results. For the current flow, an automatic
rank truncation according to equation (3) was deemed appropriate. However,
this criterion might not be ideal for other flows or data types. Hence, we
recommend a repeated examination of the rank sensitivity for each new problem
to be examined. The rank should be selected based on an iterative procedure
that yields a robust and accurate outcome. We caution against using strong
truncation that filters out the low energy-containing acoustic waves.

  4. 4.

Sample your results at a high-enough rate to avoid aliasing. Following the
Nyquist criteria, the sampling rate should be two to three times the maximum
frequency of interest.

  5. 5.

While for a robust DMD variant and settings the time window has little
influence, a reasonable time sequence should be nonetheless employed. A
minimum of two buffet cycles should be used for the analysis, where the time
window spans an integer number of cycles. For a short sequence, time windows
that span a non-integer number of cycles risk skewing the results.

  6. 6.

We recommend using the standard operator definition according to equation
(13). This definition is much more computationally tractable than the
optimized operator definition and more consistent than the so-called unitary
definition. Whereas UDMD does identify the buffet frequency robustly, the
corresponding mode occasionally exhibits strange artifacts. We posit UDMD as a
very promising variant, which requires further investigations to uncover its
full potential.

  7. 7.

The optimal mode amplitude (equation (20)) should be employed. It is more
accurate than the standard definition and less expensive than the sparsity-
promoting approach.

  8. 8.

Use the integral criterion according to equation (24) to select the modes. The
approach is best suited to handle quickly-vanishing modes that contribute very
little to the full reconstruction.

The provided workflow including the numerical setup, the datasets, and all
routines for data processing, modal analysis, and visualization are fully
available in complementary repository.

We hypothesize that this recipe is also applicable to other fluid flow
problems. The authors are currently pursuing this postulation. Other possible
future investigations include the exploration of other DMD variants that
‘wrap’ around the five-step framework we have introduced. Examples include
bagging DMD [39], multiscale DMD [40], or higher-order DMD [35]. Higher-order
DMD could better capture the nonlinearities in the data while keeping the
memory requirements tractable. However, these ‘wrapper’ algorithms are also
more computationally demanding and deliver results that are more complex to
interpret.

## Appendix

### Details about the compressible inner product

In this section, we provide some details on the choice of compressible inner
products introduced by Rowley et al. [33]. The starting point is the total
energy per unit mass, which is defined as the sum of internal energy $\hat{u}$
and kinetic energy $k$:

$e=\hat{u}+k.$ (25)

The kinetic energy per unit mass in Cartesian coordinates is given by:

$k=0.5\left(u_{x}^{2}+u_{y}^{2}+u_{z}^{2}\right).$ (26)

Assuming ideal gas and isentropic flow, the internal energy may be expressed
in terms of the local speed of sound $a$ [33]:

$\hat{u}=\frac{a^{2}}{\gamma\left(\gamma-1\right)},$ (27)

where $\gamma$ is the adiabatic index. Substituting equations (26) and (27) in
(25), weighting by the density $\rho$, and integrating over the domain
$\Omega$ yields:

$\int_{\Omega}\rho\left[\frac{a^{2}}{\gamma\left(\gamma-1\right)}+0.5\left(u_{x}^{2}+u_{y}^{2}+u_{z}^{2}\right)\right]\mathrm{d}v\approx\mathrm{const.}$
(28)

The integral is only approximately constant because the flow is not perfectly
isentropic in the case of transonic shock buffet and there might be some
netflux of energy across the boundary of $\Omega$.

### Time-continuous linear operator

In the time-continuous domain, the problem of advancing the state vector in
time is posed as:

$\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}=\mathbf{\mathcal{A}x},$ (29)

where $\mathbf{\mathcal{A}}$ is the continuous counterpart of $\mathbf{A}$.
Both operators are related as:

$\mathbf{A}=\exp{\left(\mathbf{\mathcal{A}}\Delta t\right)},$ (30)

and similarly, the relation between the eigenvalues is:

$\mathbf{\Lambda}=\exp{\left(\mathbf{\Omega}\Delta t\right)},$ (31)

where the diagonal of $\mathbf{\Omega}$ contains the eigenvalues of
$\mathbf{\mathcal{A}}$. Equation (29) has an analytical solution [10]:

$\mathbf{x}(t)=\mathbf{\Phi}\mathrm{exp}\left(\mathbf{\Omega}t\right)\mathbf{b}.$
(32)

Noticing that the eigenvalues $\omega_{i}$ in $\mathbf{\Omega}$ are usually
complex numbers, the exponential term may be split up and reformulated using
Euler’s formula:

$\mathrm{exp}\left(\omega_{i}t\right)=\mathrm{exp}\left(\Re(\omega_{i})t\right)\left[\mathrm{cos}\left(\Im(\omega_{i})t\right)+j\mathrm{sin}\left(\Im(\omega_{i})t\right)\right].$
(33)

Therefore, the $i$th eigenvalue’s real part $\Re(\omega_{i})$ is interpreted
as growth rate of mode $i$, and the imaginary part $\Im(\omega_{i})$ yields
the mode’s characteristic frequency $2\pi f_{i}=\Im(\omega_{i})$.

### UDMD buffet mode

Figure 17 shows the buffet modes resulting from UDMD for three different SVD
ranks. While the eigenvalue corresponding to the buffet mode remains nearly
identical, the mode undergoes significant changes. There is no obvious
tendency for these changes regarding the selected rank. In contrast, the
buffet mode resulting from the standard DMD operator remains almost identical
over a wide range of $r$.

Figure 17: UDMD and DMD buffet modes for varying rank truncation.

## Acknowledgments

The authors gratefully acknowledge financial support by the German Research
Foundation (DFG) received within the research unit FOR 2895 “Unsteady flow and
interaction phenomena at high-speed stall conditions”. We also acknowledge the
discussions with Steve Brunton, Benjamin Herrmann, and Peter Baddoo.

## References

  * Giannelis et al. [2017] Giannelis, N. F., Vio, G. A., and Levinski, O., “A review of recent developments in the understanding of transonic shock buffet,” _Progress in Aerospace Sciences_ , Vol. 92, 2017, pp. 39–84.  10.1016/j.paerosci.2017.05.004.
  * Ohmichi et al. [2018] Ohmichi, Y., Ishida, T., and Hashimoto, A., “Modal Decomposition Analysis of Three-Dimensional Transonic Buffet Phenomenon on a Swept Wing,” _AIAA Journal_ , Vol. 56, No. 10, 2018, pp. 3938–3950.  10.2514/1.J056855.
  * Poplingher et al. [2019] Poplingher, L., Raveh, D. E., and Dowell, E. H., “Modal Analysis of Transonic Shock Buffet on 2D Airfoil,” _AIAA Journal_ , Vol. 57, No. 7, 2019, pp. 2851–2866.  10.2514/1.J057893.
  * Das et al. [2020] Das, A., Carrese, R., Marzocca, P., and Levinski, O., _Flow Diagnostics of Transonic Shock Buffet Phenomenon of a Supercritical Airfoil using Dynamic Mode Decomposition_ , 2020.  10.2514/6.2020-1988.
  * Das et al. [2022] Das, A., Marzocca, P., Levinski, O., and Das, R., _Investigation onto Deep Transonic Buffet Condition of a Supercritical Airfoil using Multiresolution Dynamic Mode Decomposition_ , 2022.  10.2514/6.2022-1956.
  * Masini et al. [2020] Masini, L., Timme, S., and Peace, A. J., “Analysis of a civil aircraft wing transonic shock buffet experiment,” _Journal of Fluid Mechanics_ , Vol. 884, 2020, p. A1.  10.1017/jfm.2019.906.
  * Feldhusen-Hoffmann et al. [2021] Feldhusen-Hoffmann, A., Lagemann, C., Loosen, S., Meysonnat, P., Klaas, M., and Schröder, W., “Analysis of transonic buffet using dynamic mode decomposition,” _Experiments in Fluids_ , Vol. 62, No. 4, 2021, p. 66.  10.1007/s00348-020-03111-5.
  * Zhao et al. [2020] Zhao, Y., Dai, Z., Tian, Y., and Xiong, Y., “Flow characteristics around airfoils near transonic buffet onset conditions,” _Chinese Journal of Aeronautics_ , Vol. 33, No. 5, 2020, pp. 1405–1420.  10.1016/j.cja.2019.12.022.
  * Schmid [2010] Schmid, P. J., “Dynamic mode decomposition of numerical and experimental data,” _Journal of Fluid Mechanics_ , Vol. 656, 2010, p. 5–28.  10.1017/S0022112010001217.
  * Kutz et al. [2016a] Kutz, J. N., Brunton, S. L., Brunton, B. W., and Proctor, J. L., _Dynamic Mode Decomposition_ , Society for Industrial and Applied Mathematics, Philadelphia, PA, 2016a.  10.1137/1.9781611974508.
  * Schmid et al. [2011] Schmid, P. J., Li, L., Juniper, M. P., and Pust, O., “Applications of the dynamic mode decomposition,” _Theoretical and Computational Fluid Dynamics_ , Vol. 25, No. 1, 2011, pp. 249–259.  10.1007/s00162-010-0203-9.
  * Wu et al. [2021] Wu, Z., Brunton, S. L., and Revzen, S., “Challenges in dynamic mode decomposition,” _Journal of The Royal Society Interface_ , Vol. 18, No. 185, 2021, p. 20210686.  10.1098/rsif.2021.0686.
  * Schmid [2022] Schmid, P. J., “Dynamic Mode Decomposition and Its Variants,” _Annual Review of Fluid Mechanics_ , Vol. 54, No. 1, 2022, pp. 225–254.  10.1146/annurev-fluid-030121-015835.
  * Jovanović et al. [2014] Jovanović, M. R., Schmid, P. J., and Nichols, J. W., “Sparsity-promoting dynamic mode decomposition,” _Physics of Fluids_ , Vol. 26, No. 2, 2014, p. 024103.  10.1063/1.4863670.
  * Ohmichi [2017] Ohmichi, Y., “Preconditioned dynamic mode decomposition and mode selection algorithms for large datasets using incremental proper orthogonal decomposition,” _AIP Advances_ , Vol. 7, No. 7, 2017, p. 075318.  10.1063/1.4996024.
  * Kou and Zhang [2017] Kou, J., and Zhang, W., “An improved criterion to select dominant modes from dynamic mode decomposition,” _European Journal of Mechanics - B/Fluids_ , Vol. 62, 2017, pp. 109–129.  10.1016/j.euromechflu.2016.11.015.
  * Tu et al. [2014] Tu, J. H., Rowley, C. W., Luchtenburg, D. M., Brunton, S. L., and Kutz, J. N., “On dynamic mode decomposition: Theory and applications,” _Journal of Computational Dynamics_ , Vol. 1, No. 2, 2014, pp. 391–421.  10.3934/jcd.2014.1.391.
  * Kutz et al. [2016b] Kutz, J. N., Fu, X., and Brunton, S. L., “Multiresolution Dynamic Mode Decomposition,” _SIAM Journal on Applied Dynamical Systems_ , Vol. 15, No. 2, 2016b, pp. 713–735.  10.1137/15M1023543.
  * Weiner and Semaan [2022] Weiner, A., and Semaan, R., _Simulation and modal analysis of transonic shock buffets on a NACA-0012 airfoil_ , 2022.  10.2514/6.2022-2591.
  * Greenshields et al. [2010] Greenshields, C. J., Weller, H. G., Gasparini, L., and Reese, J. M., “Implementation of semi-discrete, non-staggered central schemes in a colocated, polyhedral, finite volume framework, for high-speed viscous flows,” _International Journal for Numerical Methods in Fluids_ , Vol. 63, No. 1, 2010, pp. 1–21.  10.1002/fld.2069.
  * Kurganov et al. [2001] Kurganov, A., Noelle, S., and Petrova, G., “Semidiscrete Central-Upwind Schemes for Hyperbolic Conservation Laws and Hamilton–Jacobi Equations,” _SIAM Journal on Scientific Computing_ , Vol. 23, No. 3, 2001, pp. 707–740.  10.1137/S1064827500373413.
  * van Leer [1974] van Leer, B., “Towards the ultimate conservative difference scheme. II. Monotonicity and conservation combined in a second-order scheme,” _Journal of Computational Physics_ , Vol. 14, No. 4, 1974, pp. 361–370.  10.1016/0021-9991(74)90019-9.
  * Zimmermann et al. [2018] Zimmermann, D.-M., Mayer, R., Lutz, T., and Krämer, E., “Impact of Model Parameters of SALSA Turbulence Model on Transonic Buffet Prediction,” _AIAA Journal_ , Vol. 56, No. 2, 2018, pp. 874–877.  10.2514/1.J056193.
  * Gritskevich et al. [2012] Gritskevich, M. S., Garbaruk, A., Schütze, J., and Menter, F. R., “Development of DDES and IDDES Formulations for the k-$\omega$ Shear Stress Transport Model,” _Flow, Turbulence and Combustion_ , Vol. 88, 2012, pp. 431–449.  10.1007/s10494-011-9378-4.
  * MeDevitt and Okuno [1985] MeDevitt, J. B., and Okuno, A. F., “Static and dynamic pressure measurements on a NACA 0012 airfoil in the Ames High Reynolds Number Facility,” _NASA TR NASA-TP-2485_ , 1985.
  * Iovnovich and Raveh [2012] Iovnovich, M., and Raveh, D. E., “Reynolds-Averaged Navier-Stokes Study of the Shock-Buffet Instability Mechanism,” _AIAA Journal_ , Vol. 50, No. 4, 2012, pp. 880–890.  10.2514/1.J051329.
  * Brunton and Kutz [2019] Brunton, S. L., and Kutz, J. N., _Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control_ , Cambridge University Press, 2019.  10.1017/9781108380690.
  * Gavish and Donoho [2014] Gavish, M., and Donoho, D. L., “The Optimal Hard Threshold for Singular Values is $4/\sqrt{3}$,” _IEEE Transactions on Information Theory_ , Vol. 60, No. 8, 2014, pp. 5040–5053.  10.1109/TIT.2014.2323359.
  * Hemati et al. [2017] Hemati, M. S., Rowley, C. W., Deem, E. A., and Cattafesta, L. N., “De-biasing the dynamic mode decomposition for applied Koopman spectral analysis of noisy datasets,” _Theoretical and Computational Fluid Dynamics_ , Vol. 31, No. 4, 2017, pp. 349–368.  10.1007/s00162-017-0432-2.
  * Scherl et al. [2020] Scherl, I., Strom, B., Shang, J. K., Williams, O., Polagye, B. L., and Brunton, S. L., “Robust principal component analysis for modal decomposition of corrupt fluid flows,” _Phys. Rev. Fluids_ , Vol. 5, 2020, p. 054401.  10.1103/PhysRevFluids.5.054401.
  * Baddoo et al. [2021] Baddoo, P. J., Herrmann, B., McKeon, B. J., Kutz, J. N., and Brunton, S. L., “Physics-informed dynamic mode decomposition (piDMD),” , 2021.  10.48550/ARXIV.2112.04307.
  * Askham and Kutz [2018] Askham, T., and Kutz, J. N., “Variable Projection Methods for an Optimized Dynamic Mode Decomposition,” _SIAM Journal on Applied Dynamical Systems_ , Vol. 17, No. 1, 2018, pp. 380–416.  10.1137/M1124176.
  * Rowley et al. [2004] Rowley, C. W., Colonius, T., and Murray, R. M., “Model reduction for compressible flows using POD and Galerkin projection,” _Physica D: Nonlinear Phenomena_ , Vol. 189, No. 1, 2004, pp. 115–129.  10.1016/j.physd.2003.03.001.
  * Brunton et al. [2016] Brunton, B. W., Johnson, L. A., Ojemann, J. G., and Kutz, J. N., “Extracting spatial–temporal coherent patterns in large-scale neural recordings using dynamic mode decomposition,” _Journal of Neuroscience Methods_ , Vol. 258, 2016, pp. 1–15.  10.1016/j.jneumeth.2015.10.010.
  * Le Clainche and Vega [2017] Le Clainche, S., and Vega, J. M., “Higher Order Dynamic Mode Decomposition,” _SIAM Journal on Applied Dynamical Systems_ , Vol. 16, No. 2, 2017, pp. 882–925.  10.1137/15M1054924.
  * Kaptanoglu et al. [2020] Kaptanoglu, A. A., Morgan, K. D., Hansen, C. J., and Brunton, S. L., “Characterizing magnetized plasmas with dynamic mode decomposition,” _Physics of Plasmas_ , Vol. 27, No. 3, 2020, p. 032108.  10.1063/1.5138932.
  * Weiner and Semaan [2021] Weiner, A., and Semaan, R., “flowTorch - a Python library for analysis and reduced-order modeling of fluid flows,” _Journal of Open Source Software_ , Vol. 6, No. 68, 2021, p. 3860.  10.21105/joss.03860.
  * Lee [1990] Lee, B. H. K., “Oscillatory shock motion caused by transonic shock boundary-layer interaction,” _AIAA Journal_ , Vol. 28, No. 5, 1990, pp. 942–944.  10.2514/3.25144.
  * Sashidhar and Kutz [2022] Sashidhar, D., and Kutz, J. N., “Bagging, optimized dynamic mode decomposition for robust, stable forecasting with spatial and temporal uncertainty quantification,” _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , Vol. 380, No. 2229, 2022, p. 20210199.  10.1098/rsta.2021.0199.
  * Dylewsky et al. [2019] Dylewsky, D., Tao, M., and Kutz, J. N., “Dynamic mode decomposition for multiscale nonlinear physics,” _Phys. Rev. E_ , Vol. 99, 2019, p. 063311.  10.1103/PhysRevE.99.063311.