# Harnessing the instability mechanisms in airfoil flow for the data-driven
forecasting of extreme events

Benedikt Barthel, Themistoklis Sapsis
Department of Mechanical Engineering,
Massachusetts Institute of Technology,
77 Massachusetts Ave., Cambridge, MA 02139 Corresponding author:
<EMAIL_ADDRESS>
###### Abstract

This work addresses the data-driven forecasting of extreme events in the flow
over a static airfoil. For certain Reynolds numbers and flow configurations,
airfoils are subject to sporadic high amplitude fluctuations in the
aerodynamic forces. These extreme excursions may be seen as prototypical
examples of the kind of unsteady and intermittent dynamics relevant to the
flow around airfoils and wings in a variety of laboratory and real-world
applications. Here we investigate the instability mechanisms at the heart of
these extreme events, and how knowledge thereof may be harnessed for efficient
data driven forecasting. Through a wavelet and spectral analysis of the flow
we find that the extreme events arise due to the instability of a specific
frequency component distinct from the vortex shedding mode. During these
events this extreme event manifold draws energy from the energetically
dominant vortex shedding flow and undergoes an abrupt inverse cascade of
energy transfer from small to large scales. We also investigate the spatial
dependence of the temporal correlation and mutual information between the
surface pressure and the aerodynamic forces, with the aim of identifying
regions of the airfoil amenable to sparse sensing and the efficient
forecasting of extremes. Building on previous work on predictive machine
learning models, we show that relying solely on the mutual information for
optimal sensor placement fails to improve model prediction over uniform or
random sensor placement. However, we show that by isolating the extreme event
frequency component offline through a wavelet transform we are able to
circumvent the requirement for a recursive long-short term memory (LSTM)
network – resulting in a significant reduction in computational complexity
over the previous state of the art. Using the wavelet pre-processed data in
conjunction with an extreme event-tailored loss function we find that our
model is capable of forecasting extreme events using only three pressure
sensors. Furthermore, we find our model to be robust to sensor location –
showing promise for the use of our model in dynamically varying applications.

##  1 Introduction

Many engineering systems are subject to rare high-amplitude fluctuations
commonly referred to as extreme events (Sapsis, 2021). While here we focus
primarily on fluid-structure interactions, such events occur in a wide variety
of systems ranging from climate systems to stock markets. Although rare,
events like gusts or rogue waves have a disproportionate affect on the fatigue
life of aircraft, naval vessels, or marine infrastructure. Due to their rare
nature, the prediction of these events is inherently challenging, especially
as they often occur in complex systems were the physical mechanisms are
unknown (Farazmand and Sapsis, 2019b). Various authors have ventured to
address this problem through strategies such as optimal sampling (Mohamad and
Sapsis, 2018; Sapsis, 2020b; Blanchard and Sapsis, 2021) or training
strategies which preferentially amplify rare events (Guth and Sapsis, 2019; Qi
and Majda, 2020; Rudy and Sapsis, 2022).

Two dimensional airfoil flow is one of the canonical test cases for the
dynamics of fluid-structure interaction and has been the subject of extensive
study for decades (Lissaman, 1983; Kim et al., 2012). However, the advent of
machine learning and data driven techniques has unlocked new lenses to study
this and other classical problems in the field of fluid dynamics (Brenner et
al., 2019; Fukami et al., 2020; Brunton et al., 2020). In particular, several
authors including Gomez et al. (2019); Maulik et al. (2020); Rudy and Sapsis
(2022) have proposed neural network models for the reconstruction of the flow
from sparse measurements. Data driven prediction’s based the airfoil surface
pressure have been of particular interest due to it’s practical measurability.
Most aerospace applications will require predictions made from sparse arrays
of sensors, and thus effective strategies for the optimal distribution of
pressure sensors is of critical importance. Hou et al. (2019); Le Provost et
al. (2020) used both data assimilation and convolutional and recursive neural
networks for the context of prediction of the leading edge suction parameter
(LESP), while Rudy and Sapsis (2022) investigated a range of neural network
models for the prediction of the drag coefficient.

Oscillator flows such as airfoil flows are generally insensitive to noise and
exhibit multiple characteristic time scales (Williamson, 1996; Symon et al.,
2018). This makes them an ideal candidate for the study of multi-scale slow-
fast type extreme events – as classified by Farazmand and Sapsis (2019a). From
this perspective, the Reynolds number regime
$\mathcal{O}\sim\left(10^{4}\right)$ is of particular interest. This regime
lies between steady laminar and fully turbulent regimes, and is especially
susceptible to highly nontrivial dynamics which depend significantly on angle
of attack and Reynolds number (Wang et al., 2014; Gopalakrishnan Meena et al.,
2018; Menon and Mittal, 2020). One tool for the study of such systems is the
continuous wavelet transform, which in a manner analogous to the short time
Fourier transform, quantifies the time varying strength of a signal’s
frequency components (Addison, 2016; Mojahed et al., 2021). Wavelet analysis
has previously been used in the context of extreme events by Cousins and
Sapsis (2014, 2016b); Bayındır (2016) and Srirangarajan et al. (2013) for the
detection of rogue waves and turbulent bursts in pipe flow respectively. The
majority of these studies have focused on spatial wavelet transforms, while
here we perform a temporal analysis. Additionally, the combination of wavelet
analysis with machine learning models using the type of output weighted
strategies discussed above remains largely unexplored.

This work conducts an investigation of the physical mechanisms driving the
extreme bursting events observed in the flow over a two-dimensional airfoil at
$R=17,500$ at constant angle of attack. We build on previous work by Rudy and
Sapsis (2022) who studied the data driven reconstruction of this flow using a
range of neural network architectures. We exploit the findings of our analysis
to design extreme event tailored – also referred to as output-weighted
(Sapsis, 2020b) – data processing and training strategies for the efficient
data-driven prediction of extreme events from optimal sparse sampling of the
surface pressure. For a broad discussion of output-weighted strategies for
neural networks applied to this and other systems see Rudy and Sapsis (2021).

The rest of the paper is organized as follows. In §2 we describe the problem
under investigation. In §3 we perform a statistical analysis of the data and
describe the physical mechanisms driving the extreme events, and how these
manifest in the observed data. In §4 we discuss the limitations of offline
optimal sensing algorithms. The main results of this work are then presented
in §5, and we provide some discussion of our findings in §6.

##  2 Problem Description

We consider a 2D direct numerical simulation of an incompressible flow around
a NACA 4412 airfoil at an angle of attack $\alpha=5^{\circ}$ and a cord length
based Reynolds number $R=17,500$. The flow is governed by the Navier-Stokes
and continuity equations,

$\frac{\partial\mathbf{u}}{\partial
t}+\mathbf{u}\cdot\nabla\mathbf{u}-\frac{1}{R}\nabla^{2}\mathbf{u}+\nabla
p=0,$ (1) $\nabla\cdot\mathbf{u}=0,$ (2)

where $\mathbf{u}\equiv[u(x,y,t),v(x,y,t)]$ is the velocity, $p(x,y,t)$ is the
pressure field, $t$ is time, and $x,y$ are the spatial dimensions parallel and
perpendicular to the free stream respectively. The simulation is carried out
using the open source spectral element code Nek5000 developed by Fischer et
al. (2008) with 4368 elements at spectral order 7 and a convective outflow
boundary condition (Dong, 2015). We use the same data set as Rudy and Sapsis
(2022) who report that further refinement of the numerical grid did not
meaningfully alter the results. At this Reynolds number the flow is
susceptible to intermittent, yet non-periodic turbulent bursts which manifest
as high amplitude fluctuations in the drag coefficient. Therefore we focus on
two observables: the surface pressure and the drag coefficient. The former is
a practically measurable quantity and will serve as the input to our model,
while the latter encodes the extreme events and will serve as the model
output.

Throughout this work we define $s$ as a generalized measure of arc length
measured clockwise from the leading edge (as shown in figure 1). For example,
$s\in[0,0.5)$ refers to the upper surface of the airfoil and $s\in[0.5,1)$
refers to under side. The surface pressure is saved at 100 equally spaced
locations around the airfoil surface. A visualization of the airfoil flow, the
simulation grid and the arc length measure are summarized in figure 1, and we
refer the interested reader to Rudy and Sapsis (2022) for a more detailed
discussion of the numerical method.

The aerodynamic forces are computed using skin friction and surface pressure
according to

$\mathbf{F}(t)=\oint_{s}\left(\bm{\tau}(s,t)+\mathbf{n}P(s,t)\right)ds=L(t)\hat{\mathbf{e}}_{x}+D(t)\hat{\mathbf{e}}_{y}.$
(3)

Here $t$ is time, and $x$ and $y$ represent the directions parallel and normal
to the free stream respectively. The lift and drag coefficients are then
defined as

$C_{L}(t)\equiv\frac{L(t)}{\rho U_{\infty}^{2}c},~{}\ ~{}\ ~{}\ ~{}\ ~{}\ ~{}\
~{}\ ~{}\ ~{}\ C_{D}(t)\equiv\frac{D(t)}{\rho U_{\infty}^{2}c}.$ (4)

To distinguish the extreme events from the background vortex shedding we apply
a Gaussian smoothing operation to the time series of the drag coefficient to
extract the non-periodic behaviour,

$q(t)\equiv\left(K*C_{D}\right)(t).$ (5)

where $K(t^{\prime})\propto exp\left(-(t^{\prime}/2f_{v})^{2}\right)$ is a
Gaussian smoothing kernel. $f_{v}=1.44$ is the most energetic frequency and
corresponds to periodic vortex shedding. Moving forward we simply refer to
$q(t)$ as the drag. In addition to the raw pressure signal, we also consider a
version of the pressure with the same Gaussian filter applied,

$\tilde{P}(s,t)\equiv\left(K*P(s,t)\right)(s,t).$ (6)

We refer to $P(s,t)$ and $\tilde{P}(s,t)$ as the raw and filtered pressures
respectively. In general we will consider the pressure measured at subset of
discrete sensor locations, and thus treat the surface pressure as a vector
valued quantity $\mathbf{P}(t)\in\mathbb{R}^{n}$, where $n$ is the number of
sensors. An illustrative example of the drag as well as the raw and filtered
pressures is shown in figure 2.

This flow was previously investigated by Rudy and Sapsis (2022) using a deep
long-short term memory (LSTM) network. Those authors considered a variety of
input observables, and found that the extreme events could be predicted from a
range of different observables including full and reduced order descriptions
of the flow field as well as surface pressure. This suggest that the extreme
events are a result of an underlying physical instability inherent in the
governing equations. The primary focus of this work is to identify and exploit
this mechanism for optimal sensing and forecasting. Practically, we aim to
predict future extreme events observed in the drag from sparse measurements of
the surface pressure as efficiently as possible. In other words we seek a data
driven map

$\mathbf{P}(t)\rightarrow q(t+\tau)$ (7)

for maximum lead time $\tau$, with minimal
$\operatorname{dim}\left(\mathbf{P}\right)$, and at minimal computational
cost.

Figure 1: From top left: computational grid, snapshot of vorticity and airfoil
geometry, with arclength measure $s$. Dots represent every other sensor
location. Image adapted with permission from Rudy and Sapsis (2022). Figure 2:
Illustrative example of the data. From top to bottom: filtered drag
coefficient, raw pressure signal, and filtered pressure signal. Pressure data
is taken at a single representative sensor location $25\%$ of the way along
the upper surface.

##  3 Statistical Analysis and Physical Mechanisms of Extremes

In order to gain a deeper insight into the dynamics of the flow and mechanisms
driving the extreme events, we first perform a detailed statistical analysis
of the data. We first analyze the surface pressure, as this will serve as the
basis of our modeling efforts, then in §3.3 we analyze the vorticity field –
both globally and locally along the airfoil surface – to further probe the
extreme event dynamics. Before presenting our results, we first review some
definitions we use throughout the following sections. For a signal $x(t)$ with
discrete values $x_{i}$ and distribution $X$, we define the mean $\mu_{x}$,
variance $\sigma^{2}_{x}$, and the probability density function $f_{X}(x)$.
For two signals $x(t)$ and $y(t)$ the covariance is defined as

$\sigma_{xy}\equiv\operatorname{cov}(X,Y)=\frac{1}{n-1}\sum_{i=1}^{n}\left(x_{i}-\mu_{x}\right)\left(y_{i}-\mu_{y}\right).$
(8)

To further quantify the connection between two signals we also make use the
mutual information (MI) defined as

$I(X,Y)\equiv\int_{y}\int_{x}f_{X,Y}(x,y)\log\left(\frac{f_{X,Y}(x,y)}{f_{X}(x)f_{Y}(y)}\right)dxdy,$
(9)

where $f_{X,Y}$ is the joint probability density function of $X$ and $Y$. The
MI is the Kullback–Leibler (KL) divergence between the joint probability
distribution and the product of the marginal probability distributions – it
quantifies the error in the assumption that two distributions $X$, and $Y$ are
uncorrelated. We also propose the “extreme event conditioned mutual
information”, defined as the mutual information integrated only over values of
the output greater than two standard deviations from the mean: $y>2\sigma_{y}$
– all values of the input $X$ are included –

$I_{EE}(X,Y)\equiv\int_{x}\int_{y>2\sigma_{y}}f_{X,Y}(x,y)\log\left(\frac{f_{X,Y}(x,y)}{f_{X}(x)p_{Y}(y)}\right)dxdy.$
(10)

We choose a cut-off of two standard deviations, however we found that the
results were not sensitive to changes of $\pm\sigma$. For the results
presented in this section the probability density functions in (9) and (10)
are approximated using Monte Carlo estimation using 50,000 samples and the
relevant integrals are then carried out using trapezoidal integration.

###  3.1 Mutual Information Structure

To investigate the spatial dynamics of the surface pressure we compute the
covariance and mutual information matrices:
$\operatorname{cov}\left(P(s,t),P(s^{\prime},t)\right)$ and
$I\left(P(s,t),P(s^{\prime},t)\right)$ for both the raw pressure $P$ and the
filtered pressure $\tilde{P}$. These quantify the information shared between
different locations along the airfoil. The covariance matrices and mutual
information matrices are shown in figures 3(a) and 3(b) respectively. The left
plot shows the raw pressure signal and the right shows the filtered pressure.
We notice that the results for the raw and filtered pressure are qualitatively
similar, and thus the following discussion applies to both.

These results reveal three distinct regions. First, the underside of the
airfoil, $0.5<s<1.0$. This region displays a high degree of mutual information
and strong correlation. Second, the front section of the upper surface,
$0<s<0.3$. This region exhibits similar features as the underside: strong
mutual information and correlation, however, in this region the mutual
information drops off much more quickly with separation between the sensor
locations. These results imply that these regions are amenable to sparse
sensor distribution, since any additional sensor is unlikely to contribute new
information. Note also that due to the airfoil having a non-zero angle of
attack we see strong negative correlation between the upper and lower surface
pressures. However, there is little mutual information between the upper and
lower surfaces. This implies that measurements on one surface do not
necessarily provide information about the other. The exception to this is the
rear part of the upper surface, $0.3<s<0.5$. In this region there is little to
no mutual information and significant variation in the correlation. As a
result, this region likely requires relatively higher sensor density. We note
that the transition point between the first two regions, $s=0.3$, coincides
with the point of flow separation (see figure 1). Therefore the increased
disorder observed in section 3 is likely due to to the complexity and
increased unsteadiness of the flow in this region.

We also compute the standard and extreme event mutual information between the
surface pressure and the drag coefficient: $I\left(P(s,t),q(t+\tau)\right)$
and $I_{EE}\left(P(s,t),q(t+\tau)\right)$. These are plotted in figure 4 for
the raw and filtered pressure signals for a range of $\tau$. As we are
interested in the spatial variation of these quantities, to ease comparison we
normalize each by its maximum value. In all cases we do not observe strong
dependence on the lead time $\tau$. As with the intra-pressure mutual
information we see strong spatial dependence in the region $0.3<s<0.5$. For
both the raw and filtered pressure signal, the extreme event mutual
information is (locally) peaked in this region. On the other hand, for the
standard mutual information this region lies in the trough of the spatial
distribution. This suggests that the mechanisms driving the extreme events are
strongest in the separation region. However, due to their rarity, this
connection is not reflected in the standard mutual information profile. Next,
we analyze this extreme event mechanism, and its connection to the extreme
event mutual information in more detail.

(a)

(b)

Figure 3: Pressure covariance matrix cov$\left(P(s),P(s^{\prime})\right)$ (a)
and mutual information matrix $I\left(P(s),P(s^{\prime})\right)$ (b) . Raw
pressure, $P(s,t)$ (left), filtered pressure, $\tilde{P}(s,t)$ (right). Here
$s$ is the generalized arc length around the airfoil measured clockwise from
the leading edge (Fig 1). Figure 4: Mutual information (left) and extreme
event mutual information (right) between raw (top row) and filtered (bottom
row) pressure signal along the airfoil surface and the output drag coefficient
for $\tau=0,3,7$ (blue circles, red triangles, green crosses) as a function of
arc length $s$. All curves are normalized by their maximum value.

###  3.2 Extreme Event Mechanisms

Extreme events occurring in dynamical systems are known to arise due to a
variety of factors – not all of which are fully understood. One class of
dynamical system known to give rise to extreme events are slow-fast multi-
scale systems (Farazmand and Sapsis, 2019b). In such cases, the system evolves
on two or more manifolds which have significant separation of characteristic
time scales. At most times, the system evolves along the slower manifold.
Occasionally, the trajectory may encounter an instability of this slow
manifold, resulting in the trajectory rapidly approaching the fast manifold.
Once the unstable region has passed the system relaxes back to the slow
manifold. Such phenomena are often observed as sporadic high amplitude bursts
(Farazmand and Sapsis, 2019b).

Airfoil flow is an example of such a multi-scale system. Such systems have
multi-peaked spectral content – or in other words they have multiple
characteristic frequencies. At this Reynolds number there are two high
frequencies (fast manifolds), the vortex shedding frequency, $f_{v}=1.44$
corresponding to the energetically dominant oscillatory flow, and a second
frequency corresponding to the extreme event manifold, $f_{e}=0.4$. Figure 5
shows the standard and pre-multiplied temporal Fourier power spectrum of the
filtered surface pressure defined as

$\mathcal{P}(s,f)\equiv\int P(s,t)e^{-ift}dt,$ (11)

and

$\mathcal{P}_{pm}(s,f)\equiv f\mathcal{P}(s,f),$ (12)

respectively. The latter is useful for visualizing higher frequency content as
it de-emphasizes the slow dynamics ($f\rightarrow 0$). In the standard power
spectrum there is a clear maximum close to $f=0$, corresponding to the slow
dynamics. The extreme event frequency (manifold) is also evident in the plain
spectrum, but is best seen in pre-multiplied spectrum, which exhibits a clear
peak around $f=0.4$. We show the spectrum of the filtered pressure as the
vortex shedding frequency at $f=1.44$ is much stronger than either the slow or
extreme event dynamics and obscures these when included in the spectrum. Note
also the increased magnitude in the region $0.3<s<0.5$ consistent with the
results of §3.1.

The connection between this frequency and the extreme events is best
interpreted through the wavelet transform. The wavelet transform allows for
the visualization of the time varying strength of a signal’s frequency
content. The wavelet transform has been used in the past to identify extreme
events by for example Srirangarajan et al. (2013) to detect bursts in pipe
flow and Cousins and Sapsis (2014, 2016a) for the early detection of rogue
waves. The continuous wavelet transform (CWT) of a signal $x(t)$ is defined as

$\hat{X}(f,t)\equiv\mathcal{W}\left(x(t)\right)=\sqrt{\frac{f}{f_{c}}}\int^{\infty}_{-\infty}x(s)\psi\left(f\frac{s-t}{f_{c}}\right)ds.$
(13)

Here $\psi(t)$, is the wavelet function, and $f_{c}$ is the wavelet specific
center frequency. The wavelet function is not unique, but must satisfy several
conditions including finite energy and localized support (Addison, 2016). Here
we use the Morlet wavelet,

$\psi(t)=e^{-t^{2}/2}\cos(5t).$ (14)

Moving forward we refer to the wavelet transform of the pressure signal as
$\hat{P}(s,f,t)$ where $f$ is the frequency. The wavelet transform of the
pressure signal at $s=0.34$ is shown in the upper panel of figure 6. This
location corresponds to the peak in the spatial power spectrum in figure 5.
Figure 6 clearly shows the bursts of energy at $f=f_{e}=0.4$.

We define the extreme event indicator $\bm{\gamma}$ as the wavelet coefficient
which maximizes the spectrogram of the filtered pressure signal, i.e. for
$f=f_{e}$

$\gamma(s,t)\equiv|\hat{P}(s,f_{e},t)|.$ (15)

In the lower panel of figure 6 we show the clear correlation between
$\bm{\gamma}$ and the extreme drag events. This connection is even further
highlighted in figure 7 where we compare the spatial dependence of the norm

$\|\gamma(s)\|=\sqrt{\int|\gamma(s,t)|^{2}dt},$ (16)

to the standard and extreme mutual information profiles (also shown in figure
4). As we previously found that these mutual information profiles do not
depend significantly on $\tau$, we show only the distribution for $\tau=0$.
Notice that the norm of the wavelet coefficient peaks in the same region of
the airfoil, $0.3<s<0.5$, as the extreme event mutual information between the
pressure signal and the drag. This suggests that the extreme events observed
in this flow are indeed of the multi-scale system type, and that the
connection (as quantified by the extreme event mutual information) is a
reflection of the strength of the extreme event frequency. In §5 we show that
preprocessing the pressure signal to extract this extreme event manifold
allows for a drastic reduction in the model complexity and data required for
accurate forecasting.

Figure 5: Standard (left) and pre-multiplied (right) power spectrum of
filtered surface pressure.  Figure 6: Absolute value of wavelet transformed
pressure signal at $s=0.34$ (upper panel). Extreme event indicator
$\bm{\gamma}=|\hat{P}|_{f=f_{e}}$ at same location (red triangles) compared to
drag coefficient (black line) (lower panel).  Figure 7: Spatial distribution
of norm of wavelet coefficient evaluated at extreme event frequency (black
line) compared to standard and extreme event mutual information between raw
and filtered pressure signal and drag (blue circles). Raw pressure (top row),
filtered pressure (bottom row), standard mutual information (left column),
extreme event mutual information (right column). Blue lines are the same as in
figure 4.

###  3.3 Flowfield Analysis

In §3.1 and §3.2 we identified the unstable extreme event manifold and its
connection to the separation region of the airfoil. Here we zoom out and
analyze the full boundary layer to better understand the mechanisms at the
heart of the instability and the subsequent extreme drag events. To this end
we follow the time evolution of the vorticity field,

$\omega(x,y,t)\equiv\frac{\partial v}{\partial x}-\frac{\partial u}{\partial
y}$ (17)

over the course of a single extreme event from $t=911-931$, with a peak at
$t=921$. Figure 8 shows 8 snapshots of the vorticity field over this time
interval – for clarity we focus on the region near the boundary layer. The
corresponding values of the instantaneous drag coefficient are shown in the
top left panel of the same figure. The red markers in the latter represent the
time instances of the 8 vorticity snapshots. We see clear evidence of boundary
layer separation/disorder during the time instance corresponding to the peak
in drag coefficient. Although we show only a single extreme event here, this
disordered behaviour of the boundary layer was observed during the peak of all
extreme events.

The transient dynamics observed in the full vorticity field are subtle, and
thus to better understand these transient dynamics we compute the wavelet
transform of the entire vorticity field,

$\hat{\Omega}(x,y,t,f)=\mathcal{W}\left(\omega(x,y,t)\right).$ (18)

This allows us to investigate the component of the vorticity evolving with the
extreme event frequency identified in §3.2 – which is not the dominant
energetic contributor to the full field, and is thus liable to be obscured in
the snapshots in figure 8.

####  3.3.1 Global Dynamics

To analyze the dynamics of the boundary layer we first consider the wavelet
transform of the full flowfield – focusing on the the extreme event frequency,

$\hat{\Omega}_{e}(x,y,t)\equiv\hat{\Omega}(x,y,t,f)|_{f=f_{e}=0.4},$ (19)

and the vortex shedding frequency,

$\hat{\Omega}_{v}(x,y,t)\equiv\hat{\Omega}(x,y,t,f)|_{f=f_{v}=1.44}.$ (20)

The wavelet component associated with the extreme event frequency (19) and the
vortex shedding frequency (19) are plotted in figures 9 and 10 respectively
for the same time instances as in figure 8. In the former, we see clear
evidence of a coherent structure with relatively small characteristic spatial
length scale which undergoes a transient instability resulting in a temporary
loss of coherence during the extreme event before recovering as the drag
coefficient returns to its nominal state. The vortex shedding mode has a much
larger characteristic length scale – on the order of the vortical structures
seen in figure 8 – and does not appear to undergo any significant changes
during the extreme event.

To further illustrate the dynamics of these two frequency components we
compute the temporal correlation function

$r_{\alpha}(t_{0},t)\equiv\left|\int\int\hat{\Omega}^{*}_{\alpha}(x,y,t_{0})\hat{\Omega}_{\alpha}(x,y,t)dxdy\right|,$
(21)

where ∗ denotes the complex conjugate and $\alpha=e,v$ and the integration is
performed over the entire domain. For the special case where $t_{0}=t$ this is
equivalent to the $L_{2}$ norm of the wavelet mode

$\|\hat{\Omega}_{\alpha}\|^{2}(t)\equiv\int\int|\hat{\Omega}_{\alpha}(x,y,t)|^{2}dxdy.$
(22)

These metrics respectively quantify the temporal evolution of the shape
(length scale) and magnitude of the vorticity at a specific temporal
frequency.

The correlation (21) is plotted in the upper panel of figure 11 for both
$\hat{\Omega}_{e}$ (blue circles) and $\hat{\Omega}_{v}$ (red triangles). We
fix $t_{0}=911$ as a representative snapshot corresponding to the vorticity
structure during the quiescent periods – however, any quiescent time instance
could be used. We clearly observe a systematic and drastic loss of coherence
in the extreme event mode during the spikes in the drag coefficient.
Inspection of figure 9 suggests that this is at least in part due to an
increase in the dominant spatial length scale. The coherence of the vortex
shedding mode actually fluctuates – with what further analysis reveals to be
at close to the extreme event frequency – with an amplitude that slightly
increases during the extreme events, however no drastic loss of coherence is
observed. The significance of this fluctuation is not immediately clear, but
it suggests some interaction between the two frequency components – the
investigation of which is the focus of ongoing research.

The central and lower panel of figure 11 show the time evolution of the norm
(22) of the extreme event mode (blue circles) and vortex shedding mode (red
triangles). The central panel shows the large difference in magnitude between
these two frequency components – the vortex shedding mode generally contains
an order of magnitude more energy than the extreme event mode. The lower panel
compares these norms (normalized by their value at $t_{0}$) to the drag
coefficient (black). For reference, in the lower panel we also plot the
normalized surface pressure extreme event wavelet coefficient,
$\gamma(s,t)_{s=0.34}$ (15) (green squares). As in the previous sections, the
location $s=0.34$ is chosen as it is located within the separation region.
While as previously noted, the magnitude of $\gamma$ peaks in sync with the
extreme drag events, the global norm of the extreme event component of the
vorticity, $\hat{\Omega}_{e}$, drops in magnitude during the same time
intervals. In contrast to both of these, the dominant vortex shedding mode,
(20), is significantly more stable and exhibits a much smaller relative drop
in magnitude during the extreme events.

####  3.3.2 Dynamics on the Airfoil Surface

To better understand the instability of the extreme event mode and the
associated transfer of energy we also compute the spatial Fourier transform of
the temporally wavelet transformed vorticity evaluated at the airfoil surface.
In other words we compute

$\tilde{\Omega}_{s}(t,f,k_{s})\equiv\left|\mathcal{F}\left(\hat{\Omega}(x(s),y(s),t,f)\right)\right|=\left|\int\hat{\Omega}(x(s),y(s),t,f)e^{-ik_{s}s}ds\right|$
(23)

where $k_{s}$ is the spatial wavenumber with respect to the arclength $s$
along the airfoil surface defined in §2 and $x(s)$ and $y(s)$ are the
coordinates of that surface. This reduces the spatial dimensions from two to
one, and thus allows us to visualize the transfer of energy between various
spatial and temporal scales as a function of time. The isocontours of this
quantity are plotted in figure 12 over a time horizon covering two extreme
events – the drag coefficient is also plotted for reference. This plot
succinctly summarizes the observations discussed above. First we see that
during the quiescent periods, the energy of the extreme event mode, which
actually seems to meander slightly about $f=0.4$, is concentrated at a
wavenumber $k_{s}\approx 20$. Then during the extreme events, energy is drawn
from the higher frequency vortex shedding mode, which serves as an energy
reservoir, leading to the instability of the extreme event mode which abruptly
transfer its energy to a lower wave number $k_{s}\approx 10$. This appears in
figure 12 as the “pinching off” of the isocontours during the spikes in the
drag. This is the manifestation of loss of coherence and increase in spatial
scale of the extreme event mode observed in figure 9 and quantified in the
upper panel of 11. We note in closing that this phenomenon of the extreme
events evolving on a manifold distinct from the energetically dominant one has
been observed in other systems such as for example Kolmogorov flow (Farazmand
and Sapsis, 2017; Wan et al., 2018). In that case the energy of the flow is
dominated by a specific triad of spatial wavenumbers, however projecting the
flow onto this triad fails to predict the extreme energy dissipation events
observed in that flow. Similarly, in the case of the airfoil flow considered
here, extracting only the dominant vortex shedding dynamics would miss the
extreme event dynamics entirely. However, our findings are contrary to the far
more common phenomenon of instabilities transferring energy from large to
small scales.

Figure 8: Drag coefficient (top left) and temporal snapshots of vorticity
field traversing a single extreme event. Red markers in top left plot indicate
snapshot time instances.  Figure 9: Temporal snapshots of the vorticity
wavelet coefficient evaluated at the extreme event frequency
$\hat{\Omega}_{e}(x,y,t)$. Same time snapshots as figure 8. Note the loss of
coherence and increase in spatial scales during the extreme event. Figure 10:
Temporal snapshots of the vorticity wavelet coefficient evaluated at the
vortex shedding frequency $\hat{\Omega}_{v}(x,y,t)$. Same time snapshots as
figures 8 and 9. No significany spatial changes are observed during the
extreme event, indicating that vortex shedding mode acts as energy reservoir.
Figure 11: Top panel: normalized correlation function (21) for $t_{0}=911$ for
$\hat{\Omega}_{e}$ (blue circles) and $\hat{\Omega}_{v}$(red triangles).
Center panel: vorticity wavelet coefficient norms $\|\hat{\Omega}_{e}\|^{2}$
(blue circles), $\|\hat{\Omega}_{v}\|^{2}$ (red triangles). Lower panel:
normalized vorticity wavelet coefficient norms (same markers as central panel)
and normalized surface pressure wavelet coefficient corresponding to extreme
event frequency, $\gamma(s,t)_{s=0.34}$ (15) (green squares). In the upper and
lower panel the drag coefficient $q(t)$ is shown in black.  Figure 12:
Isocontours of vorticity amplitude on the airfoil surface as a function of
temporal frequency, $f$, spatial wavenumber, $k_{s}$, and time, $t$. Blue
arrows represent direction of energy transfer.

##  4 Offline Sparse Sensor Placement

The analysis of mutual information structure described in §3 indicates that
certain sections of the airfoil are statistically more informative of the drag
coefficient. To test the practical implications of this discovery we first
propose an offline strategy to optimally select sensor locations. Such an
algorithm does not require actually training the neural network. Therefore, it
can be thought of as a prepossessing step which allows us to optimally design
the network architecture prior to training. At each iteration, the sampling
algorithm, which is outlined in algorithm 1, selects the next best sensor
location by maximizing a cost function referred to as an acquisition function.
Throughout this work we use the term “acquisition function” strictly in
connection with such a sampling strategy, and the term “cost function” to
refer to the cost function used to train a given model.

For a given application, the choice of acquisition function is not obvious,
see for example Chaloner and Verdinelli (1995); Sapsis (2020a); Yang et al.
(2021). In this framework sensor locations are selected sequentially, and
therefore we seek locations which are maximally informative of the drag
coefficient and minimally redundant with respect to the previously placed
sensors. Thus, we propose the following two acquisition functions based on the
previously defined standard and extreme event mutual information

$a^{1}_{j+1}(s,\textbf{s}^{*}_{j})=\frac{I\left(P(s),q\right)}{\frac{1}{j}\sum_{k=1}^{j}I\left(P(s),P(s_{k}^{*})\right)},$
(24)
$a^{2}_{j+1}(s,\textbf{s}^{*}_{j})=\frac{I_{EE}\left(P(s),q\right)}{\frac{1}{j}\sum_{k=1}^{j}I\left(P(s),P(s_{k}^{*})\right)}.$
(25)

The numerator – the mutual information between the pressure signal with the
drag coefficient – rewards predictive capability. The denominator – the
average of the intra-pressure sensor mutual information – penalizes
redundancy. This second condition ensures that sensors are not placed in
locations which do not contribute information not already encoded in
previously placed sensors.

There is no unique way to quantify the connection between a prospective sensor
location and the previously placed sensors, and the average used here is only
one option. Therefore, we also considered a second acquisition function where
the arithmetic mean in the denominator of (24) and (25) is replaced with a
geometric mean, but we did not observe significantly different results. A more
exhaustive study of candidate functions is beyond the scope of this work, and
so for the sake of brevity we restrict ourselves to (24) and (25). Going
forward we refer to any results obtained using this algorithm as offline-
mutual-information ($OMI_{N}$) where $N$ is the number of sensors.

Algorithm 1 Optimal Placement of $N$ sensors

while $j<N$ do
$s^{*}_{j+1}=\operatorname{argmax}a_{j+1}(s,\textbf{s}^{*}_{j})$ update
$\mathbf{s}_{j+1}^{*}=[\mathbf{s}^{*}_{j},s^{*}_{j+1}]$

end while

###  4.1 Results: Sensor Placement

We apply algorithm 1 with acquisition functions (24) and (25) to our data set
to compute the first six optimal sensor locations. Because the results of
figure 3(b) indicate that the general behaviour of the mutual information is
not dependent on the lead time $\tau$, we fix $\tau=0$. Additionally, in order
to facilitate comparison with Rudy and Sapsis (2022), we consider only the raw
pressure signal. The acquisition function landscape for at each iteration is
plotted in figure 13. The optimal senor locations after each iteration are
then summarized in figure 14.

The globally optimal sensor locations are simply the points of maximum
standard and extreme event mutual information – these are located at
approximately $s=0.15$ for (24) and $s=0.3$ for (25) respectively. However,
inspection of figure 13 indicates that in the latter case the acquisition
function landscape does not display any significant variation along the
airfoil, calling into question the viability of (25) as a practical metric for
optimal sensor placement. For the standard mutual information case, (24), at
iteration 2-4 the acquisition function exhibits multiple local maxima of
roughly equal value in the region $0.3<s<0.5$. These multiple peaks are
sequentially “picked off” throughout the iterations 2-4. This phenomenon is
also observed, but to a slightly lesser extent, in the extreme event mutual
information case (25). The similarity of these local maxima mean that these 4
sensor locations should be thought of as an “optimal grouping” rather than a
strict ranking, since measurement noise or numerical errors could affect their
ordering. However, we note that adding small amounts of noise did not
significantly impact the qualitative features of the results. Figure 14
highlights that while the ordering of the sensors varies between the standard
and extreme event versions of the model, the final distribution of the optimal
sensors is qualitatively very similar.

The under side of the airfoil $(0.5<s<1)$ is completely ignored by the
algorithm until iteration 5 for the standard case and iteration 6 for the
extreme event case, where a strong maximum is observed just downstream of the
leading edge. This solitary underside sensor near the leading edge is
consistent with the mean pressure profile observed in the flow over an
inclined airfoil. The mean pressure gradient (w.r.t. arc length) is generally
significant along the upper surface but relatively weak along the lower
surface. Therefore, a single sensor can capture a significant amount of the
information of the pressure field along the lower surface, since once the jump
in pressure across the leading edge is established there is not much more to
be gained from further probing the pressure along the underside of the
airfoil.

Figure 13: OMI acquisition function landscape for sensors 1-6 from left to
right, top to bottom. Black stars indicates maximum point. Standard mutual
information (blue squares), extreme event mutual information (red triangles).

(a)

(b)

Figure 14: Optimal sensor locations after each iteration $N$. Standard mutual
information (left), extreme event mutual information (right) (a).
Visualization of optimal sensor locations on airfoil, standard mutual
information (blue), extreme event mutual information (red) (b). Sensor
locations are labeled by $N$.

###  4.2 Results: Evaluation

To test the efficacy of the proposed offline sensor placement algorithm we
train the same LSTM network described in Rudy and Sapsis (2022) using the
first 5 optimal sensor locations predicted using algorithm 1 with acquisition
function (24). Due to the similarities of the sensor locations predicted by
(24) and (25) and the high computational cost of training the network we omit
the predictions of (25) from this analysis. We compare those results to those
in Rudy and Sapsis (2022) using 50 sensors spaced equally around the airfoil.
The network architecture is

$\rightarrow FC32\rightarrow LSTM32\rightarrow LSTM32\rightarrow
FC32\rightarrow FC16\rightarrow FC8\rightarrow FC4\rightarrow FC1\rightarrow$
(26)

where FC and LSTM stand for ‘fully-connected’ and ‘long-short term memory’
respectively, and the swish activation function is applied between each layer
(Ramachandran et al., 2017). In order to isolate the effects of the sensor
placement, we make no changes to the architecture or other than the input
dimension and utilize the same training strategies as Rudy and Sapsis (2022).
Training is conducted using $70\%$ of the data, with the remaining 30 $\%$
split evenly between validation and testing. The model was trained over 3
restarts using 140 history points and a mean square error loss function until
the validation error failed to decrease for 10 epochs – no regularization was
used. The interested reader is referred to Rudy and Sapsis (2022) for a more
detailed description of the network architecture and training strategy.

We compare three different models: the reference case from
rudy_prediction_2022 using 50 sensor locations, the $OMI_{5}$ model, as well
as a second reference case using 5 uniformly spaced sensors – all three use
the raw pressure data as an input. The last case is included to verify that
any potential benefit of our algorithm is actually due to the algorithmic
placement of the sensors and not simply a reflection of oversampling by
rudy_prediction_2022. The three models are summarized in table 1.

To compare the predictive capabilities of the models we compute both the mean
square error (MSE) of the model prediction as well as the maximum adjusted
area under the precision-recall curve – a metric introduced by
guth_machine_2019 which quantifies the accuracy of extreme event prediction.
The area under the precision-recall curve is then defined as

$\alpha(\chi)=\int_{0}^{1}S(R,\chi)dR,$ (27)

where the event rate $\chi$ is defined as the probability that the output
exceeds some threshold, the precision, $S$, the ratio of correct event
predictions to total event predictions, and the recall, $R$, is the ratio of
correct event predictions to the actual number of events. The maximum adjusted
value is then defined as

$\alpha^{*}=\operatorname{max}_{\chi}\left(\alpha(\chi)-\chi\right).$ (28)

When the value of $\alpha^{*}$ is large (approaches unity) the model is very
good at predicting rare events, alternatively, when the value approaches zero
a model does no better than a guess based on the aggregate frequency of
extreme events.

Figure 15 compares the mean absolute error, MAE, and $\alpha^{*}$ of the
various models for a range of lead times $\tau$. As expected, for all cases
MAE increases and $\alpha^{*}$ decreases with $\tau$ – it is more difficult to
predict the far future. Comparing the models, we first observe that there is
no significant difference between the two reference cases with 5 and 50
sensors, suggesting that the NN model developed by rudy_prediction_2022, is
amenable to far more sparse sensor distributions than is suggested by those
authors. Additionally, we see no clear distinction between the results of the
model trained using our predicted optimal sensor locations and the uniformly
sampled reference cases. This highlights the limitations of the mutual
information as a practical tool for engineering design. To further highlight
this limitation and exclude the possibility of our conclusions being
influenced by oversampling we train the model using the 5 optimal sensor
locations individually. In other words we train 5 models, each with a single
sensor (the $n^{th}$ optimal OMI prediction) as its input. The same error
metrics for this experiment are plotted in figure 16. Again we see that the
optimal sensor location $s=0.14$ performs no better, and in many cases worse,
than the suboptimal locations – see for example the value of $\alpha^{*}$ at
$\tau=0$ and $7$. These results strongly suggest that optimal sensing based
purely on mutual information does not adequately capture the extreme event
mechanism identified in §3.2 and thus fails as a practical tool for optimal
sensing.

Name | Sensor Placement | Input Observable | Number of Sensors $(N)$
---|---|---|---
$OMI_{N}$ |  $OMI$ algorithm | $P(t)$ | $N$
$uni_{5}$ | $uniform$ | $P(t)$ | $5$
$uni._{50}$ | $uniform$ | $P(t)$ | $50$
Table 1: Summary of neural net models compared in figure 15. Figure 15:
Maximum adjusted area under the precision recall curve, $\alpha^{*}$ and mean
absolute error, $MAE$, for the models summarized in table 1 as a function of
the lead time $\tau$. Figure 16: Maximum adjusted area under the precision
recall curve, $\alpha^{*}$ and mean absolute error, $MAE$, for models trained
using a single sensor as a function of the lead time $\tau$.

##  5 Wavelet preprocessing for extreme event prediction

Despite its robustness to sensor location, the LSTM network considered in §4
is expensive to train, so here we explore a different avenue of model
reduction: preprocessing the data through offline identification of the
extreme event dynamics. In §3.2 we find that the extreme event dynamics are
directly related to the dynamics of a single frequency component. Here we show
that exploiting this observation through the event indicator (15) allows for
the forecasting of the extreme events using very simple network models.

###  5.1 Methods

The extreme event indicator $\bm{\gamma}$ defined in (15) is not only highly
correlated with the bursting events, as it represents an isolated frequency
component, but is also free of noise. This makes it amenable to accurate
numerical differentiation. We therefore define the following transformation

$D_{n}\left(\bm{\gamma}\right):~{}\
\bm{\gamma}\rightarrow\left[\bm{\gamma},\dot{\bm{\gamma}},...,\bm{\gamma}^{(n)}\right],$
(29)

which allows us to track not only the value of $\bm{\gamma}$, but also its
growth rate. This is crucial as we seek to forecast bursting for nonzero lead
times $\tau$, and therefore it is imperative for the model to observe growth
and not just magnitude. The differentiation operation in (29) is essentially a
phase shift of the signal and thus aids the forecasting capabilities of the
model, i.e. the predictions for $\tau>0$. We find that for this flow, a single
derivative ($n=1$) is sufficient, and including a second derivative did not
meaningfully improve results. Note that this differentiation is applied
offline, and thus does not affect the computational cost of training the
model. We propose a preprocessing procedure described in algorithm 2, in which
we replace (7) with

$\bm{\Gamma}(t)\rightarrow q(t+\tau).$ (30)

Here the input data is defined as

$\bm{\Gamma}\equiv
D_{1}\left(\bm{\gamma}\right)=[\bm{\gamma},\dot{\bm{\gamma}}],$ (31)

and $\bm{\gamma}$ is defined in (15). We utilize a fully connected neural
network $f:\mathrm{R}^{n}\rightarrow\mathrm{R}^{1}$ with layers

$\rightarrow FC8\rightarrow FC16\rightarrow FC16\rightarrow FC8\rightarrow
FC1\rightarrow,$ (32)

with the swish activation function (ramachandran_searching_2017) applied
between each layer. A regularization constant of $0.01$ was applied to the
activation layers – no regularization was applied to the kernels. Note that
unlike the LSTM network (26) utilized by rudy_prediction_2022, this network
does not map sequences to sequences, it simply maps values of $\bm{\Gamma}$ at
time $t$, to values of $q$ at time $t+\tau$. To train the model we use both a
standard and output-weighted mean absolute error loss function,

$MAE=\sum_{j}|\hat{q}_{j}-q_{j}|,$ (33)
$MAE_{OW}=\sum_{j}\frac{|\hat{q}_{j}-q_{j}|}{f_{q}(q_{j})}.$ (34)

Here $\hat{q}_{j}$ is the model prediction, $q_{j}$ is the training data, and
$f_{q}(q_{j})$ is the probability density function of the training data
evaluated at $q_{j}$. rudy_output-weighted_2021 found that output-weighted
loss functions significantly improve prediction of outlier events in a variety
of flows including airfoil and Kolmogorov flow. While those authors use the
mean square error, we find that in our case the mean absolute error
consistently performed slightly better. Our model is summarized graphically in
figure 17.

In order to quantify the uncertainty of our model, we perform an ensemble
analysis resulting in a mean prediction $\hat{q}(t)$ and variance
$\sigma_{q}(t)$. We set aside $80\%$ of the data set for training, and for
each iteration of the ensemble we randomly select $75\%$ of that training data
($60\%$ of the total) to use for training. The remaining $20\%$ of the data is
used for testing. All results presented here are computed exclusively using
this test data.

Algorithm 2 Wavelet Preprocessing

1\. Input data: $\mathbf{P}(t)$

2\. Compute spectrogram: $|\hat{\mathbf{P}}(t)|$

3\. Select Maximum: $\bm{\gamma}(t)\equiv\underset{f\neq
f_{v}}{\max}\left(|\hat{\bm{P}}(f,t)|\right)$

4\. Differentiate: $\bm{\Gamma}=[\bm{\gamma},\dot{\bm{\gamma}}]$

5\. Output: $\bm{\Gamma}(t)$

Figure 17: Illustration of wavelet preprocessing algorithm 2. Elements
outlined in red represent the contributions of this work.

###  5.2 Results: Basis Comparison

To illustrate the advantages of the wavelet basis and the output-weighted loss
function we compare the model predictions using the three basis types: $P$,
$\tilde{P}$, and $[\bm{\gamma},\dot{\bm{\gamma}}]$ and the two loss functions
$MAE$ and $MAE_{OW}$. In all cases we use four evenly spread sensor locations
at $s=0.05,0.35,0.65,0.95$. The first two sensor location represent areas
identified in §3.1 and §3.2 as predictive of the drag. The latter two
locations are chosen as to not neglect the underside of the airfoil. In all
cases train an ensemble of 10 models for 200 epochs.

Figures 18 and 19 compare the model prediction of the drag coefficient to the
true value for $\tau=0$ and $\tau=7$. At $\tau=0$, the filtered pressured
model and the wavelet model, shown in blue and red respectively, perform well,
with only the raw pressure model, shown in green, suffering from significant
noise corruption. In fact, in this case the smooth pressure model slightly
outperforms the wavelet basis – best seen by comparing the predicted
probability density functions in figure 20. This is because there exists an
accurate linear mapping from the smooth pressure $\tilde{\mathbf{P}}(t)$ to
the drag $q(t)$, and in taking the wavelet transform of the filtered pressure
some information is lost – limiting the potential accuracy of the wavelet
model. See appendix A.1 for a brief discussion on this linear mapping.
Furthermore, we observe that for zero lead time, the models trained using the
$MAE_{OW}$ loss function perform slightly better than those trained using the
$MAE$ loss – The benefit is most pronounced for the raw pressure data.

The benefits of the wavelet basis and the output weighted loss functions do
not become apparent until considering non-zero lead times – a reflection of
the nonlinearity of the time shift operation $q(t)\rightarrow q(t+\tau)$. In
this case, when using the standard $MAE$ loss, all three models entirely fail
to capture the extreme events. However, using the output-weighted loss the
wavelet basis retains much of the accuracy observed for $\tau=0$, while the
performance of the raw and filtered pressure models deteriorate significantly.
The filtered pressure model still traces the occurrence of the extreme events,
but suffers from significant noise corruption leading to number of false
positive predictions. This phenomenon is even more pronounced for the raw
pressure model, which, as expected, suffers from even greater noise
corruption. The distinction between the three basis types is less evident in
the predicted probability density functions shown in 20. Here we see again
that the models trained using the standard $MAE$ loss fail to capture the
tails of the distribution entirely, however with the $MAE_{OW}$ loss both the
filtered pressure and wavelet models capture the general shape of the
distribution. The wavelet model does however capture the small peaks around
$q=-1.5$ and $q=3$ slightly better than the others.

To quantify the forecasting capabilities of each model we track the number of
extreme events predicted as a function of time. For this purpose we define an
extreme event as a local maximum, whose value is more than 2 standard
deviations greater than the mean. Thus, a time instant $t_{j}$ is considered
to represent an extreme event $t_{EE}$ if it satisfies the following
conditions,

$t_{EE}:~{}\ t_{j}~{}\ s.t.~{}\ \left.\left[\frac{\partial q}{\partial
t}\right|_{t_{j}}=0~{}\ \&~{}\ q(t_{j})>2\sigma_{q}\right].$ (35)

We then define the number of extreme events $N_{EE}(t_{1},t_{2})$ as the
number of extreme events in the interval, $t_{1}$ to $t_{2}$, or more
explicitly,

$N(t_{n},t_{0})=\sum_{j=j_{0}}^{j_{n}}\delta_{t_{j},t_{EE}},$ (36)

where $t_{j}$ is treated as a discrete series and $j_{0}$ and $j_{n}$ are the
indices of $t_{n}$ and $t_{0}$ respectively. In order to avoid over-penalizing
noise, we enforce a minimum separation of the identified extreme events equal
to the characteristic period of the extreme event frequency:
$T_{EE}=1/f_{EE}$. While it relies on two user defined parameters: the extreme
event threshold and the minimum peak separation, this metric provides a useful
quantification of the forecasting capabilities of each model.

This metric is plotted in figure 21 for the filtered pressure and wavelet
models for $\tau=[0,3,7,10]$ – we omit the raw pressure model due to its poor
performance. Since the models trained using the $MAE$ loss and perform so
poorly for $\tau>0$ we plot only the results obtained using the $MAE_{OW}$
loss. Again we see that for $\tau=0$ both models perform similarly, with the
filtered pressure model slightly outperforming the wavelet model. However, as
the lead time $\tau$ increases, the wavelet model retains much of it’s
accuracy, while the filtered pressure model on the other hand dramatically
overestimates the number of extreme events. This is due to the significant
noise in the filtered pressure model, the magnitude of which is often
comparable to the underlying signal.

We also compute the MAE, MSE, $\alpha^{*}$, defined in (28), as well as the
error in the total number of extreme events predicted. These are plotted in
figure 22. Again, both pressure models overestimate the number of extreme
events defined by (35), however for the aggregate error metrics ($MAE$ and
$MSE$) as well as the $\alpha^{*}$ the differences are less pronounced. Both
the filtered pressure and wavelet model significantly outperform the raw
pressure model, but the difference between them is minimal.

The discrepancy in performance between the three basis types can be understood
through the simple nature of the model architecture. The large amplitude of
the high frequency fluctuations in the raw pressure is comparable and
sometimes even larger than the bursting amplitude – see figure 2 – therefore a
one-to-one map is destined to fail. This phenomenon is mitigated by filtering
the vortex shedding frequency out of the pressure data – in this case there is
indeed a linear map for zero lead time. However, for non-zero lead times the
amplitude of the fluctuations at the extreme event frequency are significant
enough to introduce significant ambiguity in a one-to-one map. Conversely, the
wavelet basis is free of noise and fluctuates on a time scale associated with
the mean time between extreme events, thereby greatly improving the
feasibility of such a simple mapping. Moving forward we exclusively use the
wavelet pre-processed input data $\bm{\Gamma}$.

(a)

(b)

Figure 18: Time series of test set predictions for raw pressure (green),
filtered pressure (blue), and wavelet (red) models, compared to true values
(black) for $\tau=0$. Loss function: $MAE$ (a), $MAE_{OW}$ (b).

(a)

(b)

Figure 19: Time series of test set predictions for raw pressure (green),
filtered pressure (blue), and wavelet (red) models, compared to true values
(black) for $\tau=7$. Loss function: $MAE$ (a), $MAE_{OW}$ (b).

(a)

(b)

Figure 20: Probability density function of test set predictions for raw
pressure (green), filtered pressure (blue), and wavelet (red) models, compared
to true values (black) for $\tau=0$ (a) and $\tau=7$(b). Loss function: $MAE$
(left), $MAE_{OW}$ (right).  Figure 21: Predicted number of extreme events by
filtered pressure model (blue) and wavelet model (Red) compared to true number
(black) for from top left $\tau=[0,3,7,10]$.  Figure 22: From top left: $MAE$,
$MSE$, error in number of predicted extreme events, and $\alpha^{*}$. Models:
raw pressure (green), filtered pressure (blue), and wavelet (red).

###  5.3 Results: Optimal Sensing

We now assess how to best exploit the wavelet preprocessing algorithm through
optimal sensor selection. Here we focus exclusively on $\tau=7$. As in §4 we
use algorithm 1, however in this case the acquisition function requires
training the network (and is minimized, not maximized). In particular we
consider the following two acquisition functions which penalize uncertainty in
the model prediction,

$a_{iu}=\frac{1}{T}\int^{T}_{0}\sigma_{q}(t)dt=\frac{1}{N}\sum^{N}_{j=1}\sigma_{q,j},$
(37)

$a_{pw}=\frac{1}{T}\int^{T}_{0}\frac{\sigma_{q}(t)}{p_{q}(q(t))}dt=\frac{1}{N}\sum^{N}_{j=1}\frac{\sigma_{q,j}}{p_{q}(q_{j})}.$
(38)

We refer to these as “integrated uncertainty” (iu) and “probability weighted”
(pw) respectively. To distinguish the effects of the loss and acquisition
functions, we perform three iterations of algorithm 1 with each acquisition
function with both the $MAE$ and $MAE_{OW}$ loss functions. As the acquisition
functions need to be evaluated at every sensor location, we use a slightly
reduced ensemble of 7 models and train over only 70 epochs during the active
search. The resulting optimal sensor locations are summarized in table 2. Once
the optimal sensor locations are found, we retrain an ensemble of 10 models
using those optimal sensors locations for 200 epochs.

Figure 23 compares the output pdf and time series of the mean predictions of
the model trained with each loss function and each acquisition function – each
using their respective optimal sensor locations. As a comparison we also
include the predictions of a reference model with three evenly spaced sensors.
For both the time series and the pdf, we observe that regardless of the
acquisition function the models trained using the standard $MAE$ loss fail,
while the models trained with the output-weighted $MAE_{OW}$ predict the
bursting events relatively accurately. To further compare the models we plot
the time series of the uncertainty bound, $\hat{q}(t)\pm\sigma_{q}$ in figure
24. Consistent with the results of §4, we observe the model to be robust to
specific sensor locations, with little distinction between the three sensor
distributions.

Interestingly, inspection of figures 23(a) and 24 indicates that the output-
weighted acquisition function (38) performs slightly worse than the others –
exhibiting some false positive fluctuations between $t=920$ and $t=960$.
Furthermore, the results using the non output-weighted acquisition function
(37) do not exhibit any meaningful improvement over the reference case. These
findings indicate that the specific locations of the sensors are of secondary
importance when compared to the effects of the output-weighted loss function
and the wavelet preprocessing. The latter of which extracts the bursting
events from the input data a priori. These results indicate that further
emphasizing extreme events through strategies such as output-weighted optimal
sensing is not only unnecessary, but could result in a loss of accuracy during
the quiescent periods.

From a practical point of view, the similarity of these results is
significant. The active regions of the flow, and thus the optimal sensor
locations predicted here and in §4, are likely to vary with Reynolds number
and angle of attack. However, aircraft experience a wide range of flow speeds
and orientations, making the robustness to sensor location a valuable asset.
These results support the possibility of a sparse sensing strategy which is
applicable for a wide range of airfoil designs and is robust to dynamic
changes in angle of attack. A parameter study over Reynolds number and flow
geometry to confirm this hypothesis is beyond the scope of the present work,
but is the topic of ongoing research.

| $IU$ | $PW$
---|---|---
$MAE$ | $[0.31,0.15,0.25]$ | $[0.09,0.84,0.05]$
$MAE_{OW}$ | $[0.21,0.15,0.60]$ | $[0.58,0.76,0.13]$
Table 2: Optimal sensor locations predicted using algorithm 1 with acquisition
functions (37) and (38) with $MAE$ and $MAE_{OW}$ loss functions.

(a)

(b)

Figure 23: Comparison of test set predictions for reference (green), $IU$
(blue), and $PW$ (red) acquisition functions to true values (black). Time
series for $MAE$ (top) and $MAE_{OW}$ (bottom) (a), and probability density
function for $MAE$ (left) and $MAE_{OW}$ (right) (b). Figure 24: Comparison of
test set predictions for reference (green), $IU$ (blue), and $PW$ (red)
acquisition functions to true values (black). Shaded area represents mean
$\pm$ one standard deviation of ensemble prediction.

##  6 Discussion

We have investigated the mechanisms driving the non-periodic bursting
phenomena observed in the two-dimensional flow over a NACA 4412 airfoil at
finite angle of attack. We have conducted a detailed analysis of the
spatiotemporal statistics of the airfoil surface pressure and its connection
to the extreme events observed in the drag force. Through a wavelet analysis
we found that the surface pressure exhibits multi-scale behaviour with three
distinct time scales. In addition to the dominant vortex shedding frequency,
the flow exhibits a slowly varying quiescent time scale and a second energetic
frequency component – at approximately one third the vortex shedding
frequency. We established that the extreme excursions of the drag first
observed by rudy_prediction_2022, correspond to instabilities of this latter
frequency component in the surface pressure.

These findings were corroborated by an analysis of the wavelet transformed
vorticity field. This analysis revealed that during quiescent times the
extreme event manifold evolves independently of the vortex shedding manifold,
however occasionally the extreme event manifold undergoes a transient
instability which links the fortunes of these two generally disparate time
scales. This instability is comprised of two steps, first the extreme event
manifold draws energy from the higher frequency vortex shedding flow, then at
the extreme event frequency, there is an abrupt nonlinear energy transfer from
smaller to larger spatial length scales. Interestingly these findings are
contrary to the far more common situation where linear instabilities transfer
energy from a slowly evolving (or stationary) mean flow to faster time scales
and smaller length scales. Therefore, while we have identified the slow-fast
system at the heart of the bursting events, the exact mechanism by which the
instabilities in the pressure and vorticity are translated to the aerodynamic
forces is not yet clear, and remains the topic of ongoing research. For
example, it is still unclear what causes the global (integrated over the full
domain) magnitude of the extreme event mode to decrease during the extreme
events – see the lower panel of figure 11 – or why the temporal correlation of
the vortex shedding mode exhibits fluctuations resembling the extreme event
frequency – see the upper panel of the same figure. Furthermore, here we have
considered only a single angle of attack, and further study is required to
establish how the orientation of the flow impacts both the active regions of
the airfoil and the extreme event frequency.

From a modeling perspective, we pursued two separate strategies. First, in §4
we investigated the implications of these results for the existing LSTM
architecture developed by rudy_prediction_2022 – which takes raw pressure as
its input. We considered an optimal sensing strategy based purely on the
statistics of the data, and therefore did not require the computationally
costly step of training model. Using the LSTM model this mutual information
based algorithm failed to predict sensor locations which performed better than
a simple uniform sensor distribution. This failure of the purely mutual
information based sensor placement demonstrated the limitations of a purely
statistical offline sampling strategy and highlighted the limitations of
mutual information as a practical tool. Additionally, the model complexity
incurred by the LSTM layers needed to process rapidly varying time series such
as the fluctuating surface pressure remains cumbersome regardless of the
sparsity of the sensor array.

Second, we also developed a preprocessing algorithm (see figure 17) to extract
the time varying magnitude of the extreme event frequency component from the
raw pressure signal. By isolating the wavelet coefficient associate with the
extreme event frequency we eliminate the high frequency fluctuations resulting
in a signal which slowly fluctuates on a time scale associated with the mean
time between extreme events – which are by definition rare. This enables the
(approximate) one-to-one mapping of the wavelet coefficient, which is free of
rapid high amplitude noise, to the drag for lead times $\tau>0$. This then
eliminates the need for a costly LSTM architecture and allows for accurate
prediction using a simply connected feed forward neural network. These results
are consistent with findings of cousins_sapsis; cousins_reduced-order_2016 who
used a spacial wavelet transform wavelet transform to extract unstable spatial
length scales to efficiently predict rogue waves in variety of dynamical
systems including the Majda–McLaughlin–Tabak model and the modified nonlinear
Schrodinger equation.

While this preprocessing drastically reduces the noise in the signal, it can,
and in our case does, eliminate some potentially useful information as well.
As noted in §5.2 and in appendix A.1 for $\tau=0$ there exists an accurate
linear map from the surface pressure to the drag coefficient. By isolating a
single wavelet coefficient, some of the information in the pressure signal is
lost leading to the wavelet model performing slightly worse than the filtered
pressure model for $\tau=0$. However, the predictions from the wavelet model
are far more robust to increasing values of $\tau$. At $\tau=7$ the
predictions of the wavelet model have degraded only slightly, while the raw
and filtered pressure models exhibit significant noise corruption. The higher
the frequency of oscillation, the more nonlinear the transformation
$q(t)\rightarrow q(t+\tau)$ – resulting in the degradation of the filtered and
raw pressure models for $\tau>0$.

This preprocessing alleviates the need for recursive or convolutional network
architectures as used by authors such as hou_machine-learning-based_2019;
rudy_prediction_2022. However, even with this highly extreme event targeted
algorithm we found that training the model using an output-weighted loss
function is necessary for accurate predictions. Most interestingly, we find
that with these training interventions the specific locations of the sensors
is of secondary importance. This is incredibly advantageous as it suggests
that the predictive capabilities of our approach are robust to dynamic changes
in angle of attack or free stream velocity – however this requires further
study. Our findings suggest that improving the prediction of rare events does
not necessarily require more complex models, but can be achieved by
identifying observables which reflect the underlying physical mechanisms and
through tailored training strategies – as also discussed by various authors
including farazmand_variational_2017; sapsis_output-weighted_2020;
rudy_output-weighted_2021; blanchard_output-weighted_2021. We believe the
herein proposed wavelet based analysis is applicable to a wide range of slow-
fast systems and remains the topic of ongoing research.

##  Appendix A Appendix

###  A.1 Linear Mapping

Here we illustrate the linear map from the vector valued filtered pressure
signal $\tilde{\mathbf{P}}(t)$ to the scalar drag coefficient $q(t)$. Let
$\tilde{\mathbf{P}}\in\mathbb{R}^{N\times 100}$ and
$\mathbf{q}\in\mathbb{R}^{N\times 1}$ be their discrete representations – $N$
is the number of data points (time steps). We then seek a linear
representation $\hat{\mathbf{q}}=\tilde{\mathbf{P}}\mathbf{a}$ that minimizes
the $L_{2}$ norm $\|\mathbf{q}-\hat{\mathbf{q}}\|^{2}$. The optimal
coefficient vector is given by

$\mathbf{a}^{*}=\tilde{\mathbf{P}}_{train}^{+}\mathbf{q}_{train},$ (39)

where + denotes the pseudo inverse, and the subscript train refers to the
subset of data used for training. Here we use the first $10\%$ of the data to
fit the regression (39) and the last $20\%$ for testing. Figure 25 compares
the predictions of the linear regression to the truth for $\tau=0$ and
$\tau=7$. For $\tau=0$ the linear prediction is indistinguishable from the
truth, while for $\tau=7$ the linear model completely fails. This is a
reflection of the highly nonlinear nature of the time shift operation
$q(t)\rightarrow q(t+\tau)$.

Figure 25: Prediction of drag coefficient from linear regression model (blue)
compared to true value (black). $\tau=0$ (upper panel), $\tau=7$ (lower
panel).

## Funding Sources

## Acknowledgments

We thank Samuel Rudy and Tanner Harms for their constructive feedback. We also
acknowledge support from the Army Research Office (grant no. W911NF-17-1-0306)
and the Air Force Office of Scientific Research (grant no. MURI
FA9550-21-1-0058)

## References

  * Addison (2016) P. S. Addison.  _The Illustrated Wavelet Transform Handbook: Introductory Theory and Applications in Science, Engineering, Medicine and Finance, Second Edition_.  CRC Press, Boca Raton, 2 edition, Dec. 2016.  ISBN 978-1-315-37255-6.
  * Bayındır (2016) C. Bayındır.  Early detection of rogue waves by the wavelet transforms.  2016\.
  * Blanchard and Sapsis (2021) A. Blanchard and T. Sapsis.  Output-Weighted Optimal Sampling for Bayesian Experimental Design and Uncertainty Quantification.  _SIAM/ASA Journal on Uncertainty Quantification_ , 9(2):564–592, Jan. 2021.  ISSN 2166-2525.
  * Brenner et al. (2019) M. P. Brenner, J. D. Eldredge, and J. B. Freund.  Perspective on machine learning for advancing fluid mechanics.  _Physical Review Fluids_ , 4(10):100501, Oct. 2019\.
  * Brunton et al. (2020) S. L. Brunton, B. R. Noack, and P. Koumoutsakos.  Machine Learning for Fluid Mechanics.  _Annual Review of Fluid Mechanics_ , 52(1):477–508, 2020.
  * Chaloner and Verdinelli (1995) K. Chaloner and I. Verdinelli.  Bayesian experimental design: A review.  _Statistical Science_ , 10(3):273–304, 1995.  doi: 10.1214/ss/1177009939.
  * Cousins and Sapsis (2014) W. Cousins and T. P. Sapsis.  Quantification and prediction of extreme events in a one-dimensional nonlinear dispersive wave model.  _Physica D_ , 280:48–58, 2014.
  * Cousins and Sapsis (2016a) W. Cousins and T. P. Sapsis.  Reduced order precursors of rare events in unidirectional nonlinear water waves.  _Journal of Fluid Mechanics_ , 790:368–388, 2016a.  ISSN 14697645.  doi: 10.1017/jfm.2016.13.
  * Cousins and Sapsis (2016b) W. Cousins and T. P. Sapsis.  Reduced-order precursors of rare events in unidirectional nonlinear water waves.  _Journal of Fluid Mechanics_ , 790:368–388, Mar. 2016b.
  * Dong (2015) S. Dong.  A convective-like energy-stable open boundary condition for simulations of incompressible flows.  _Journal of Computational Physics_ , 302:300–328, Dec. 2015\.  doi: 10.1016/j.jcp.2015.09.017.
  * Farazmand and Sapsis (2017) M. Farazmand and T. P. Sapsis.  A variational approach to probing extreme events in turbulent dynamical systems.  _Science Advances_ , 3(9):e1701533, Sept. 2017\.  doi: 10.1126/sciadv.1701533.
  * Farazmand and Sapsis (2019a) M. Farazmand and T. P. Sapsis.  Closed-loop adaptive control of extreme events in a turbulent flow.  _Physical Review E_ , 100:033110, 2019a.
  * Farazmand and Sapsis (2019b) M. Farazmand and T. P. Sapsis.  Extreme Events: Mechanisms and Prediction.  _Applied Mechanics Reviews_ , 71(5):1–19, 2019b.  ISSN 0003-6900.  doi: 10.1115/1.4042065.
  * Fischer et al. (2008) P. F. Fischer, J. W. Lottes, and S. G. Kerkemeier.  Nek5000 Webpage.  2008\.
  * Fukami et al. (2020) K. Fukami, K. Fukagata, and K. Taira.  Assessment of supervised machine learning methods for fluid flows.  _Theoretical and Computational Fluid Dynamics_ , 34(4):497–519, Aug. 2020.
  * Gomez et al. (2019) D. F. Gomez, F. Lagor, P. B. Kirk, A. Lind, A. R. Jones, and D. A. Paley.  Unsteady DMD-Based Flow Field Estimation From Embedded Pressure Sensors in an Actuated Airfoil.  In _AIAA Scitech 2019 Forum_ , AIAA SciTech Forum. American Institute of Aeronautics and Astronautics, Jan. 2019.
  * Gopalakrishnan Meena et al. (2018) M. Gopalakrishnan Meena, K. Taira, and K. Asai.  Airfoil-Wake Modification with Gurney Flap at Low Reynolds Number.  _AIAA Journal_ , 56(4):1348–1359, 2018.
  * Guth and Sapsis (2019) S. Guth and T. P. Sapsis.  Machine Learning Predictors of Extreme Events Occurring in Complex Dynamical Systems.  _Entropy_ , 21(10):925, Oct. 2019.  doi: 10.3390/e21100925.
  * Hou et al. (2019) W. Hou, D. Darakananda, and J. D. Eldredge.  Machine-Learning-Based Detection of Aerodynamic Disturbances Using Surface Pressure Measurements.  _AIAA Journal_ , 57:5079–5093, Dec. 2019.
  * Kim et al. (2012) D.-H. Kim, J.-W. Chang, and J. Chung.  Low-Reynolds-Number Effect on Aerodynamic Characteristics of a NACA 0012 Airfoil.  _Journal of Aircraft_ , May 2012.
  * Le Provost et al. (2020) M. Le Provost, W. Hou, and J. Eldredge.  Deep learning and data assimilation approaches to sensor reduction in estimation of disturbed separated flows.  In _AIAA Scitech 2020 Forum_ , AIAA SciTech Forum. Jan. 2020.
  * Lissaman (1983) P. B. S. Lissaman.  Low-Reynolds-Number Airfoils.  _Annual Review of Fluid Mechanics_ , 15(1):223–239, 1983.
  * Maulik et al. (2020) R. Maulik, K. Fukami, N. Ramachandra, K. Fukagata, and K. Taira.  Probabilistic neural networks for fluid flow surrogate modeling and data recovery.  _Physical Review Fluids_ , 5(10):104401, Oct. 2020\.  Publisher: American Physical Society.
  * Menon and Mittal (2020) K. Menon and R. Mittal.  Aerodynamic Characteristics of Canonical Airfoils at Low Reynolds Numbers.  _AIAA Journal_ , 58(2):977–980, Feb. 2020.
  * Mohamad and Sapsis (2018) M. A. Mohamad and T. P. Sapsis.  Sequential sampling strategy for extreme event statistics in nonlinear dynamical systems.  _Proceedings of the National Academy of Sciences_ , 115(44):11138–11143, Oct. 2018.
  * Mojahed et al. (2021) A. Mojahed, L. A. Bergman, and A. F. Vakakis.  New inverse wavelet transform method with broad application in dynamics.  _Mechanical Systems and Signal Processing_ , 156:107691, July 2021.  doi: 10.1016/j.ymssp.2021.107691.
  * Qi and Majda (2020) D. Qi and A. J. Majda.  Using machine learning to predict extreme events in complex systems.  _Proceedings of the National Academy of Sciences_ , 117(1):52–59, Jan. 2020.
  * Ramachandran et al. (2017) P. Ramachandran, B. Zoph, and Q. V. Le.  Searching for Activation Functions, Oct. 2017.  arXiv:1710.05941 [cs].
  * Rudy and Sapsis (2021) S. Rudy and T. Sapsis.  Output-weighted and relative entropy loss functions for deep learning precursors of extreme events, Dec. 2021.  arXiv:2112.00825 [cs, math].
  * Rudy and Sapsis (2022) S. H. Rudy and T. P. Sapsis.  Prediction of Intermittent Fluctuations from Surface Pressure Measurements on a Turbulent Airfoil.  _AIAA Journal_ , Mar. 2022.  Publisher: American Institute of Aeronautics and Astronautics.
  * Sapsis (2020a) T. Sapsis.  Output-weighted optimal sampling for Bayesian regression and rare event statistics using few samples.  _Proceedings of the Royal Society A_ , 476:20190834, 2020a.
  * Sapsis (2020b) T. P. Sapsis.  Output-weighted optimal sampling for Bayesian regression and rare event statistics using few samples.  _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 476(2234):20190834, Feb. 2020b.  ISSN 1364-5021, 1471-2946.  doi: 10.1098/rspa.2019.0834.
  * Sapsis (2021) T. P. Sapsis.  Statistics of Extreme Events in Fluid Flows and Waves.  _Annual Review of Fluid Mechanics_ , 53(1):85–111, 2021.
  * Srirangarajan et al. (2013) S. Srirangarajan, M. Allen, A. Preis, M. Iqbal, H. B. Lim, and A. J. Whittle.  Wavelet-based Burst Event Detection and Localization in Water Distribution Systems.  _Journal of Signal Processing Systems_ , 72(1):1–16, jul 2013.  ISSN 1939-8018.  doi: 10.1007/s11265-012-0690-6.
  * Symon et al. (2018) S. Symon, K. Rosenberg, S. T. M. Dawson, and B. J. McKeon.  Non-normality and classification of amplification mechanisms in stability and resolvent analysis.  _Physical Review Fluids_ , 3(5):053902, May 2018\.  ISSN 2469-990X.  doi: 10.1103/PhysRevFluids.3.053902.
  * Wan et al. (2018) Z. Y. Wan, P. Vlachas, P. Koumoutsakos, and T. Sapsis.  Data-assisted reduced-order modeling of extreme events in complex dynamical systems.  _PLOS ONE_ , 13(5):e0197704, May 2018.
  * Wang et al. (2014) S. Wang, Y. Zhou, M. M. Alam, and H. Yang.  Turbulent intensity and Reynolds number effects on an airfoil at low Reynolds numbers.  _Physics of Fluids_ , 26(11):115107, Nov. 2014\.
  * Williamson (1996) C. H. K. Williamson.  Vortex Dynamics in the Cylinder Wake.  _Annual Review of Fluid Mechanics_ , 28(1):477–539, 1996.
  * Yang et al. (2021) Y. Yang, A. Blanchard, T. Sapsis, and P. Perdikaris.  Output-Weighted Sampling for Multi-Armed Bandits with Extreme Payoffs.  _arXiv_ , 2021.