# On the Conflict of Robustness and Learning in Collaborative Machine Learning

Mathilde Raynal
SPRING Lab, EPFL  Carmela Troncoso
SPRING Lab, EPFL

###### Abstract

Collaborative Machine Learning (CML) allows participants to jointly train a
machine learning model while keeping their training data private. In scenarios
where privacy is a strong requirement, such as health-related applications,
safety is also a primary concern. This means that privacy-preserving CML
processes must produce models that output correct and reliable decisions _even
in the presence of potentially untrusted participants_. In response to this
issue, researchers propose to use robust aggregators that rely on metrics
which help filter out malicious contributions that could compromise the
training process. In this work, we formalize the landscape of robust
aggregators in the literature. Our formalization allows us to show that
existing robust aggregators cannot fulfill their goal: either they use
distance-based metrics that cannot accurately identify targeted malicious
updates; or propose methods whose success is in direct conflict with the
ability of CML participants to learn from others and therefore cannot
eliminate the risk of manipulation without preventing learning.

##  1 Introduction

Collaborative Machine Learning (CML) allows users to jointly train a machine
learning model. The main argument to favor CML instead of centralized learning
techniques is _privacy_ : CML enables training of models without the users’
data leaving their device. Users share learning proxies, e.g., gradients or
model weights [fl], which are typically referred to as updates. CML
architectures differ in the means to orchestrate updates amongst participants,
with the main two approaches being server-aided (through Federated Learning
[fl]), and peer-to-peer (through Decentralized Learning [dl]).

The privacy of these CML designs has been studied in depth [privacyDL]. In
this paper, we focus on _robustness_ –i.e., the reliability and correctness of
the model output (e.g., prediction or classification result) – a property
crucial in many fields where the privacy of CML is typically considered an
advantage. For example, in health applications [health1, health2, health3,
health4, health5, health6, health7, helath8, health9, health10], where
incorrect diagnosis can threaten the life of patient; or in autonomous driving
[fldriving, FLdeepdriving, feddrive, driving4, driving5, driving6, driving7]
where wrong predictions threaten the safety of passengers, bystanders, and the
environment. Even in fields where robustness is not essential for safety,
robustness may be needed for economic reasons, such as avoiding manipulations
that change content visibility in advertisement applications [ads], or change
credit risk prediction [credit].

Lack of robustness can have two causes. First, from non-malicious failures,
such as users having low-quality data or an unreliable network connection that
prevents them from participating in the network. Second, it can be the result
of attacks performed by malicious participants whose goal is to poison the
final model and influence its behavior in some undesirable way.

In this paper, we evaluate the robustness of CML in the presence of an
adversary. Following the literature, we assume this adversary can send
arbitrary model updates to other participants during training e.g., [rfa,
Section 3.2], [krum, Section 2] [scclipping, Section 1], [rofl, Section 2.1],
[trimmedmedian, Section 1], [baffle, Section 3.C]. In this scenario, naive CML
implementations, such as averaging all received updates [fl], do not provide
robustness: the aggregation function processes its inputs in a predictable
manner, and hence a malicious user can craft updates to arbitrarily influence
the output of the aggregation [open].

Robust training techniques from the centralized setting such as examining data
quality [clean, demon, fail, roni] or testing for adversarial examples
[advexp], do not apply to CML. This is because in CML users do not share raw
and they are not guaranteed to share a common distribution. To address this
limitation, researchers have developed a number of CML-tailored sophisticated
functions to improve robustness [rfa, trimmedmedian, tm2, scclipping,
ogclipping, krum, bridge, bylan, order, cao, cert, romoa, autofed, mixing,
byrdie, xie, ubar, localpoison, fltrust, svmdl, vote, rofl]. These functions
enable participants to evaluate the model updates they receive, and reject
those that are deemed detrimental to the final model, i.e., participants learn
as if there were only honest users in the network. These functions have been
referred to as robust aggregators [scclipping], byzantine resilient defenses
[byrdie], scoring mechanism [krum], or even poisoning detectors [local].

While the efficiency of state-of-the-art robust aggregators against a
selection of attacks has been demonstrated, so far there is yet no evaluation
that discusses the effectiveness of these algorithms against an strategic
adversary that would adapt to the defense nor of their impact on learning.
This follows the typical trend in systems’ evaluation in which researchers
iterate attacks and defenses to arrive to the state-of-the-art. This was the
case in the search for robust algorithms in the centralized machine-learning
literature [carlini], and this trend is repeating in the collaborative
setting. There is a growing amount of attacks against state-of-the-art robust
aggregators [local, scclipping, alie, imp, history, limits] and new robust
aggregators protecting against state-of-the-art attacks [scclipping, local,
tm2, ogclipping, history]. None of these works provide a formal proof of
robustness guarantees against a strategic adversary. This lack of formalism is
even more problematic in light of the recent findings which demonstrate how
standard experimental setups can misrepresent robustness [pitfalls].

In this paper, we avoid this cat-and-mouse process by systematizing the study
of CML robust aggregators, and directly characterizing the limits of
collaborative robust learning in presence of malicious strategic adversaries.
This characterization reveals a fundamental limit to robustness as robust
aggregators currently define it: obtaining robustness in CML is in direct
conflict with learning. In other words, allowing participants to learn always
incur a risk of honest users being manipulated. Such risk might be acceptable
in some applications, like in advertising; or inadmissible due to potential
harms to users, like in health or self-driving applications.

Concretely, we make the following contributions:

  * •

We introduce a taxonomy for robustness in CML that systematizes the different
kinds of robust aggregation techniques in the literature. We find that most
robust aggregators belong to one of two categories: those that evaluate the
robustness of updates based on their distance, e.g., L2 distance to a
reference point; and those that evaluate the robustness of updates based on
their behavior, e.g., error rate.

  * •

We show that distance-based robust aggregators that allow for a level of
learning must also allow for adversarial manipulations of the same order of
magnitude. This implies that _distance is not an adequate proxy for evaluating
the maliciousness of updates_. We empirically demonstrate our findings in two
realistic applications of CML: learning in healthcare and autonomous driving.

  * •

We demonstrate that when using behavior-based robust aggregators, the quality
of the evaluation that users make is inversely proportional to their ability
to learn from collaboration. Thus, the only users that can make good decisions
and accurately reject malicious updates (and therefore learn in a robust
manner) are those whose prior knowledge implies that they do not benefit much
from participating in collaborative learning. Conversely, users who have the
highest learning potential cannot produce accurate robustness evaluations, and
remain vulnerable to manipulation despite their use of robust aggregation. We
empirically validate our results on two example behaviour-based aggregators.

##  2 Systematization of Robustness Metrics

A CML system consists of $n$ users who hold private local data. User $i$ holds
data $S_{i}$ formed by tuples of samples $x$ and labels $y$ sampled from the
domain $\mathcal{D}=(\mathcal{X},\mathcal{Y})$. Users are connected according
to a graph $\mathcal{G}$, in which nodes are users and edges represent
connections between users. We call this graph _topology_. A common topology in
CML is a star, used in Federated Learning [fl], in which all nodes communicate
with each other via a central entity.

In CML, users initialize their models with common public parameters $m^{0}$
and iterate over the following three steps until a stop condition is met: (1)
Local training: At the beginning of epoch $t$, users apply gradient descent on
their local model parameters $m^{t}$ using their local data $S_{i}$, and
compute an update $\theta_{i}^{t}$. (2) Communication: Users share their model
update $\theta_{i}^{t}$ with their neighbors— or the server in the case of
federated learning. (3) Aggregation: Users (or the server in federated
learning) aggregate the updates, e.g., compute their average or their median.
Users obtain an updated model $m^{t+1}$ using the result of their local
aggregation. In federated learning, the server updates the global model and
sends it to the users.

The primary goal of the collaboration is to enable users to obtain a model
which performs better on the learning task than a model trained solely on
their local data. In many cases, the application also requires this model to
be robust, i.e., the participants must not be able to influence the model in
any undesired way. While undesired is subjective and potentially relative to
the application, in all cases the goal is to prevent manipulations that result
on erroneous model outputs.

In CML, the main tool to ensure robustness of the final model is the use of
robust aggregators during the training process. Robust aggregators aim to
filter out detrimental updates before performing the aggregation, such that
only honest users contribute to training. Robust aggregators rely on an
evaluation metric to distinguish benign updates from malicious ones. We use
$R(\theta;\mathcal{D})$ to refer to the evaluation of the update $\theta$ on
the inputs from $\mathcal{D}$, using the metric $R$. The output of the
evaluation is not necessarily restricted to $\\{benign,malicious\\}$ but can
be a number that participants use in conjunction with a threshold value
$\delta$, e.g., $R(\theta;\mathcal{D})<\delta$, to identify benign updates.
Using a threshold makes the evaluation more flexible, e.g., to adapt the
robustness evaluation to the learning task or application, or to the learning
stage (e.g., use a stricter evaluation at the final steps of the training
process).

A robust aggregator is a concrete instantiation of the metric $R$ and a means
to evaluate the exchanged updates, i.e., each proposition of a metric $R$ is a
proposition of a definition of robustness. The metric $R$ must have two
properties. First, $R$ must capture the difference between benign and
malicious updates such that the robust aggregator can effectively reject
malicious updates: if an update does not influence the model in any way that
leads to erroneous results, then $R(\theta;\mathcal{D})$ must be 0 (or below
$\delta$, when using a threshold). Second, participants must be able to
correctly compute $R$, even when they have access to a limited dataset and
limited knowledge on the learning task.

###  2.1 Existing Approaches to Robustness

The metric $R$ should, ideally, be informative about the maximal error of the
model $m^{t+1}$ resulting of the use of the update $\theta$ on inputs from the
domain $\mathcal{D}$, i.e.,
$R(\theta,\mathcal{D})\doteq\underset{(x,y)\sim\mathcal{D}}{max}[m^{t+1}(x)-y]$.
Yet, computing such metric in the collaborative setting is impossible because
it requires computing the output of the model on all samples from
$\mathcal{D}$. In addition to being a potentially unbounded computational
task, individual users do not know the whole domain of inputs (by definition
of the problem).

As a consequence, metrics $R(\theta,\mathcal{D})$ proposed in the literature
resort to the evaluation of proxies that can be easily computed. These proxies
include (but are not limited to): average accuracy of the update on a fixed
test set [xie], invariance of prediction under input perturbation
[invariance], smoothness of decision regions [advex] when the update is a
model itself, or more generally the low norm of weights [l2] or L2-distance to
the average of previously received updates [trilemma].

In general, the evaluation of an update $\theta$ is based on a comparison
between this update and a reference that we denote as $\theta^{*}$. The
reference can be a gradient or a model, typically of the same format as
$\theta$. In the literature, we find that the reference $\theta^{*}$ is either
fixed before training given prior knowledge of the participants, or computed
dynamically per epoch given information available at the time of evaluation.
In the latter case, the evaluator can use local knowledge (e.g., a local
update), or knowledge obtained from the other participants’ in the network
(e.g., the average of all updates, assuming that the majority of users is
honest). The closer the update being evaluated ($\theta$) is to this reference
($\theta^{*}$), according to $R$, the more likely $\theta$ is to be benign.
When $\theta=\theta^{*}$, $R(\theta,\mathcal{D})=0$.

In the CML literature, we find two classes of metrics used to assess
‘closeness’, which we formalize in the next two sections. Systematizing robust
aggregators through classes of metrics allows us to reason about the
properties of a class instead of individual robust aggregators. Such
systematization also enables to evaluate the robustness of new aggregators
without the need of an analysis from scratch.

####  2.1.1 Distance-based Evaluation

The most common class of robust aggregators are those that evaluate a
characteristic of the update itself, e.g., the update weights [rfa,
trimmedmedian, tm2, scclipping, ogclipping, krum, bridge, bylan, order, cao,
cert, romoa, autofed, mixing, byrdie, rofl, fltrust, ubar, sentinel].
Approaches in this class measure the closeness between the update $\theta$ and
a reference $\theta^{*}$ using a distance metric $\Delta$ —i.e. a metric that
has reflexivity, non-negativity, symmetry, and respects the triangle
inequality. We refer to this class of robust aggregators as distance-based. We
formalize them in Definition 2.1 and give examples of possible $\theta^{*}$
and $\Delta$ in Examples LABEL:ex:krum, LABEL:ex:rofl, and LABEL:ex:scc.

######  Definition 2.1 (Distance-based robustness).

We say a robust aggregator is distance-based if its evaluation function can be
expressed as:

$R(\theta,\mathcal{D})\doteq\Delta(\theta,\theta^{*})$

for a distance metric $\Delta$ and a reference update $\theta^{*}$.


# Deep learning-based predictive modelling of transonic flow over an aerofoil

Li-Wei Chen, Nils Thuerey
Technical University of Munich, D-85748 Garching
<EMAIL_ADDRESS>
(22 March 2024)

###### Abstract

Effectively predicting transonic unsteady flow over an aerofoil poses inherent
challenges. In this study, we harness the power of deep neural network (DNN)
models using the attention U-Net architecture. Through efficient training of
these models, we achieve the capability to capture the complexities of
transonic and unsteady flow dynamics at high resolution, even when faced with
previously unseen conditions. We demonstrate that by leveraging the
differentiability inherent in neural network representations, our approach
provides a framework for assessing fundamental physical properties via global
instability analysis. This integration bridges deep neural network models and
traditional modal analysis, offering valuable insights into transonic flow
dynamics and enhancing the interpretability of neural network models in
flowfield diagnostics.

##  1 Introduction

Self-sustained shock wave motions on aerofoils are closely linked with the
occurrence of buffeting, causing pronounced shock-induced fluctuations.
Several comprehensive reviews have shed light on this topic, including works
by [1], [2] and [3]. To elucidate the underlying mechanism of self-sustained
shock wave motion, [4] proposed a phenomenological model that considers a
closed feedback loop between the mean position of the shock wave and the
trailing edge, which serves as the acoustic source. This feedback model has
gained widespread recognition and has been extensively employed to scrutinize
transonic flows over aerofoils [5, 6, 7, 8].

More quantitatively, [9] highlighted the significance of global instability in
the onset of the transonic buffet. Building upon this, [10] further explored
the impact of various disturbances and perturbations on buffet mode stability.
Their investigation uncovered medium-frequency unsteadiness resembling Kelvin-
Helmholtz-type instability, in addition to the well-documented low-frequency
shock unsteadiness. Furthermore, [11] utilized resolvent analysis to connect
the buffet mechanism to the shock foot, shedding light on low-Reynolds number
compressible aerodynamics and the influence of perturbations in the vicinity
of the trailing edge. The existence of both low and high frequencies in the
context has been supported by the research conducted by [12] and [13]. They
investigated unsteady transonic flows around a two-dimensional aerofoil,
revealing four transition stages, including steady flow, spontaneous
unsteadiness, a critical Mach number range, and two instability processes,
with a focus on the influence of compressibility effects on these transitions.
To explain the coexistence of the unsteadiness, they utilized Landau’s model
as an interpretive framework.

In addition to the aforementioned advancement of theoretical models, data-
driven modeling methods demonstrate significant potential for modeling and
analyzing transonic flows over aerofoils, encompassing techniques such as
proper orthogonal decomposition [8, 14], dynamic mode decomposition [15, 16],
and data-driven resolvent analysis [17, 18]. Those modal approaches facilitate
the discernment of physically critical features within the complexity of high-
dimensional and nonlinear phenomena. However, their limitations are also well-
documented, as they can struggle with accurately predicting the behavior of
highly nonlinear systems [19].

Recently, the application of deep learning models for predicting the unsteady
evolution of complex fluid dynamics has garnered significant attention and
opened up new possibilities for performing data-driven predictive modeling of
fluid dynamics. Some researchers have harnessed the potential of deep learning
in conjunction with numerical solvers to effectively predict the non-
stationary characteristics of turbulent flows [20, 21, 22, 23, 24]. In
contrast, others have pursued end-to-end supervised training and achieved
remarkable outcomes [25, 26, 27]. Nonetheless, a distinct gap persists in the
availability of neural network-based predictive models for unsteady transonic
flows, especially in situations where interactions between shock waves and
vortical flows play a significant role [28, 29]. As previously mentioned, the
transonic flow over an aerofoil may exhibit two distinct frequency components,
depending on the free stream Mach number [12]. The presence of these dual
frequencies and their interaction mechanism present a significant challenge
for neural network prediction models. Therefore, a specialized approach is
required to capture and predict the complex interaction between these high and
low-frequency phenomena.

Our main contribution is to showcase that a deep learning-based model can
capture the intricate dynamics of high-resolution unsteady transonic flows.
Remarkably, the trained model can predict the evolution of the entire
flowfield surrounding an aerofoil, in scenarios with unseen free-stream Mach
numbers and even with an initial time-mean flowfields. Once trained, its
adaptability allows for handling diverse tasks with multiple objectives.
Furthermore, we show that the resulting model really learned the essence of
the underlying physical system by scrutinizing global instability. Compared to
others based on a Koopman framework [30] or a dynamical mode decomposition
method [31], we show for the first time that a global instability analysis can
be derived for neural networks with a linearized operator-based approach. In
contrast to traditional approaches, the differentiable nature of neural
network representations provides a natural avenue for evaluating fundamental
physical properties of learned dynamics with global instability analysis. This
fusion bridges the gap between deep neural network models and traditional
modal analysis, offering fresh insights into transonic flow dynamics. Our
innovative approach not only offers new perspectives in flowfield diagnostics
but also lays the groundwork for improving the interpretability of neural
network models.

This paper is organized as follows. The physical system and data generation
are briefly presented in §2. The neural network architecture and training
procedure will be described in §3. The detailed experiments and results are
then given in §4 and concluding remarks in §7.

##  2 Physical system and data generation

###  2.1 Physical system

We consider the transonic flow over a NACA0012 aerofoil at a low Reynolds
number, i.e. $15000$. The governing equations of the physical system are
compressible Navier-Stokes equations, which are non-dimensionalized with the
freestream variables (i.e., the density $\rho_{\infty}$, speed of sound
$a_{\infty}$, and the chord length of the aerofoil $c$). The conservative form
about density $\rho$, momentum $\rho u_{\mathrm{i}}$ and total energy $\rho E$
can be expressed in tensor notation as

$\displaystyle\frac{\partial\rho}{\partial t}+\frac{\partial\rho
u_{\mathrm{i}}}{\partial x_{\mathrm{i}}}$ $\displaystyle=0$
$\displaystyle\frac{\partial\rho u_{\mathrm{i}}}{\partial
t}+\frac{\partial\rho u_{\mathrm{i}}u_{\mathrm{j}}}{\partial x_{\mathrm{j}}}$
$\displaystyle=-\frac{\partial p}{\partial
x_{\mathrm{i}}}+\frac{\partial\tau_{x_{\mathrm{i}}x_{\mathrm{j}}}}{\partial
x_{\mathrm{j}}}$ $\displaystyle\frac{\partial\rho E}{\partial
t}+\frac{\partial(\rho E+p)u_{\mathrm{i}}}{\partial x_{\mathrm{i}}}$
$\displaystyle=\frac{\partial(-q_{\mathrm{i}}+u_{\mathrm{j}}\tau_{x_{\mathrm{i}}x_{\mathrm{j}}})}{\partial
x_{\mathrm{i}}}$

where the shear stress and heat flux terms are defined as

$\tau_{x_{\mathrm{i}}x_{\mathrm{j}}}=\mu\frac{M_{\infty}}{Re_{\infty}}[(\frac{\partial
u_{\mathrm{i}}}{\partial x_{\mathrm{j}}}+\frac{\partial
u_{\mathrm{j}}}{\partial x_{\mathrm{i}}})-\frac{2}{3}\frac{\partial
u_{\mathrm{k}}}{\partial x_{\mathrm{k}}}\delta_{\mathrm{i}\mathrm{j}}]$

and

$q_{x_{\mathrm{i}}}=-\frac{\mu}{Pr}\frac{M_{\infty}}{Re_{\infty}(\gamma-1)}\frac{\partial\Theta}{\partial
x_{\mathrm{i}}}.$

Here, Reynolds number is defined as
$Re_{\infty}=\rho_{\infty}\sqrt{u_{\infty}^{2}+v_{\infty}^{2}}c/\mu_{\infty}$;
$\gamma$ is the ratio of specific heats, 1.4 for air; the laminar viscosity
$\mu$ is obtained by Sutherland’s law (the function of temperature); laminar
Prandtl number is constant, i.e. $Pr=0.72$. The relation between pressure $p$
and total energy $E$ is given by

$p=(\gamma-1)\bigl{[}\rho E-\frac{1}{2}\rho
u_{\mathrm{i}}u_{\mathrm{i}}\bigr{]}.$

Note also that from the equation of state for a perfect gas, we have $p=\rho
a^{2}/\gamma$ and temperature $\Theta=a^{2}$. No turbulence model is employed
as we perform 2D high-resolution quasi-direct numerical simulations.

The finite-volume method numerically solves the equations using the open-
source code CFL3D [32, 33]. As depicted in figure 1(a), the mesh resolution
$1024\times 128$ is kept the same for cases, i.e. 128 grid cells in the wall-
normal direction, 320 grid cells in the wake, and 384 grid cells around the
aerofoil surface. The convective terms are discretized with the third-order
upwind scheme and viscous terms with a second-order central difference. An
inflow/outflow boundary based on one-dimensional Riemann invariant is imposed
at about 50 _c_ away from the aerofoil in the (_x_ , _y_) plane. The grid
stretching is employed to provide higher resolution near the surface and in
the wake region, and the minimal wall-normal grid spacing is $6\times 10^{-4}$
to ensure $y^{+}_{n}<1.0$. A no-slip adiabatic wall boundary condition is
applied on the aerofoil surface. The non-dimensional time step is
$0.008c/{U_{\infty}}$.

(a)

(b)

Figure 1: (a) Computational mesh with $1024\times 128$ grid cells; (b) The
contour of the y-component velocity at $M_{\infty}=0.825$ on the fine grids on
the left column v.s. filtered one on coarse grids ($256\times 64$) on the
right column. The bottom shows the fields with the corresponding meshes.
Figure 2: Typical flowfields in the training set.

(a)

(b)

Figure 3: (a) Root mean square values of lift coefficients and (b)
Statistically mean drag coefficients; shaded regions represent $\pm 3$
standard deviations.

###  2.2 Data generation

In the transonic regime, airflow behavior becomes more complex due to the
formation of shock waves and supersonic and subsonic flow areas on the
aerofoil surfaces. To cover all possible flow regimes of interest, the samples
in the training dataset are generated at 18 different Mach numbers (i.e.,
$M_{\infty}$=0.725, 0.75, 0.76, 0.775, 0.778, 0.78, 0.79, 0.8, 0.81, 0.82,
0.825, 0.835, 0.85, 0.86, 0.87, 0.88, 0.89, and 0.9).

As shown in figure 2, at $M_{\infty}=0.75$, the primary source of unsteadiness
in the flowfield stems from high-frequency vortex shedding. When the free-
stream Mach number approaches 0.8, a series of compression waves coalesce to
form a strong shock wave, and the flow structures are dominated by alternately
moving shock waves along the upper and lower sides of the aerofoil (e.g.
$M=0.825$ and $M=0.85$ in figure 2). In this regime, the flowfield exhibits a
combination of high-frequency vortex shedding and low-frequency unsteadiness
attributed to shock-wave dynamics, leading to significant fluctuations in lift
and drag coefficients, as shown in the figure 3. When the Mach number exceeds
$0.89$, the intense shock waves stabilize on both surfaces, and unsteadiness
is predominantly driven by the shedding of vortices in the wake.

To evaluate the learned model, we choose five representative conditions for
the test set: $M_{\infty}=0.71$, $M_{\infty}=0.755$, and $M_{\infty}=0.772$
are in the vortex shedding regime; $M_{\infty}=0.84$ case is in the shock
oscillation regime; $M_{\infty}=0.905$ is the case when the shockwave is
stable. The lift and drag coefficients of the test cases are marked as red
symbols in figure 3.

Following the completion of the transient phase, snapshots are stored at
intervals of every four simulation steps (i.e.,
$dt_{sampling}=0.032c/U_{\infty}$ and $f_{Nyquist}=15.625U_{\infty}/c$). As
the flowfield variables are defined on grid cells, we downsample them by
4$\times$ and 2$\times$ in the circumferential and wall-normal directions with
spatial filtering. Specifically, the density and pressure, as well as the grid
cell center coordinates, are subject to cell-area-weighted averaging, while
the velocities undergo mass-weighted filtering, known as Favre filtering,
preserving the grid’s topological structure and maintaining a one-to-one block
interface in the wake region. The typical filtered field variables are shown
in figure1(b). In each case (or sequence), there are 400 to 1000 snapshots,
which capture at least two periods of low-frequency motion in the flowfield.

Through the coordinate transformation, variables in the physical space can be
uniquely mapped into the computational space and back [34]. Hence, we define a
unit cube of curvilinear coordinates as the canonical space in which neural
networks operate [35]. As our training focuses on a single aerofoil and
computational mesh, the geometric information, in the form of a transformation
matrix and airfoil profile, is not included as part of the neural network’s
input.

##  3 Predictive modelling with deep learning

###  3.1 Learning task

We frame the learning task as the discovery of a function $f_{\theta}$ (with
weights $\theta$) that can map the input state variables $\mathbf{q}_{t}$,
given the condition $M_{\infty}$, to an approximation of the subsequent state
$\mathbf{q}_{t+\Delta t}$,

$\mathbf{q}_{t+\Delta t}=f_{\theta}(\mathbf{q}_{t},M_{\infty}).$ (1)

Here, $\mathbf{q}$ represents the four flowfield variables, including density,
velocities in the $x$ and $y$ directions, and pressure. The learned model is
applied for a sequence of $m$ steps, generating a trajectory of states
$(\mathbf{q}_{t_{0}},\dots,\mathbf{q}_{t_{m}})$, and training is aimed at
effectively capturing the dynamic evolution of a dataset of ground truth
trajectories.$(\mathbf{\tilde{q}}_{t_{0}}\dots\mathbf{\tilde{q}}_{t_{m}})$.

###  3.2 Neural Network Architectures

As architecture for the neural networks we employ an attention-based U-Net
[36]. It consists of three encoder blocks, each progressively capturing
features from the input image through convolution with Group Normalization and
ReLU activation, and downsampling. Following the encoder blocks, there’s a
bottleneck block for information compression with a higher number of channels.
Subsequently, the architecture includes three decoder blocks, which use skip
connections to integrate features from both the bottleneck and corresponding
encoder blocks during the upscaling process. These decoder blocks gradually
reduce the number of channels, culminating in a $1\times 1$ convolutional
output layer that generates pixel-wise predictions. In contrast to the
original implementation, our current study places a specific emphasis on the
c-type mesh topology. As a result, we’ve incorporated custom paddings for all
convolution layers and bilinear upsampling operations, achieved by exchanging
the 1-1 block interface values within the mesh interface. This approach allows
us to preserve the c-type mesh topological structure throughout our feature
maps.

We train networks of varying sizes by adjusting the number of features, as
indicated in table 1. The details are provided in the appendix figure A1.

Model name | # of Parameters | # of Features in encoder & bottleneck
---|---|---
0.51M | 511332 | [16, 32, 64] [128]
2.04M | 2037956 | [32, 64, 128] [256]
8.14M | 8137092 | [64, 128, 256] [512]
Table 1: Network architectures and parameter counts

###  3.3 Rollout training

Prediction operates on the same network parameter (i.e. $\theta$) from
equation 1. The total loss over the unrolled trajectory consisting of $s$
steps is

$\mathcal{L}=\sum_{i=1}^{s}\mathcal{L}(\tilde{\mathbf{q}}_{i+s\Delta
t},f_{\theta}^{s}(\mathbf{q}_{i}))=\sum_{i=1}^{s}\mathcal{L}^{s}.$ (2)

Here, $\mathcal{L}^{s}$ represents the loss evaluated at a specific step $s$.
As the model is differentiable, we can calculate the full optimization
gradients by propagating gradients through predictions. The gradients are thus
evaluated as

$\frac{\partial\mathcal{L}^{s}}{\partial\theta}=\sum_{B=1}^{s}\bigg{[}\frac{\partial\mathcal{L}^{s}}{\partial
f_{\theta}^{s}}\bigg{(}\prod_{b=1}^{s-B}\frac{\partial
f_{\theta}^{s-b+1}}{\partial f_{\theta}^{s-b}}\bigg{)}\frac{\partial
f_{\theta}^{B}}{\partial\theta}\bigg{]}.$ (3)

The model is trained with Adam and a mini-batch size of 4. The group
normalization technique employs four groups within the convolutional layers.
To improve generalization on unseen data and enhance prediction robustness, we
introduce 0.5% Gaussian noise into the input data [37, 29]. Following [21],
the training curriculum gradually increases the number of unrolled steps,
commencing with $s=1$ and subsequently progressing to $s=4$, and then $s=9$,
followed by $s=14$. In each, a learning rate of $6\times 10^{-4}$ is applied
for the initial 250k iterations, after which it exponentially decays to
$6\times 10^{-5}$ for the rest 250k. Thus, a total of 2 million iterations are
used in the whole curriculum. Our observations indicate that the marginal
benefits diminish significantly beyond 14 roll-out steps in the current
context. Hence, we choose the model trained with $s=14$ for our analyses.

##  4 Results and discussion

###  4.1 Inference from a previously unseen condition

In the following, we discuss the overall performance of the learned models and
the predicted aerodynamic characteristics. We investigate test cases with
previously unseen Mach numbers: $M_{\infty}=0.71$, $M_{\infty}=0.755$,
$M_{\infty}=0.772$, $M_{\infty}=0.84$ and $M_{\infty}=0.905$, corresponding to
three shock-free cases, one shock-oscillation case and one stabilized shock-
wave case, respectively. Given the Mach number range of the training data,
cases of $M_{\infty}=0.71$ and $M_{\infty}=0.905$ will be considered as
“extrapolation tests”, the other three will be treated as “interpolation
tests”.

####  4.1.1 Overall evaluation

(a) Loss history over 999 prediction steps.

(b) Loss history over 2499 prediction steps.

Figure 4: Long-term $L_{2}$ errors measured using five large-network-size models (8.14M), each initialized with distinct random seeds. Case | Model name | $N=99$ | $N=999$ | $N=2499$
---|---|---|---|---
| 0.51M | 0.00146$\pm$0.000497 | 0.0138$\pm$0.0102 | /
$M=0.755$ | 2.04M | 0.000956$\pm$0.000272 | 0.00937$\pm$0.00390 | /
| 8.14M |  0.000813$\pm$0.000304 |  0.00695$\pm$0.00517 | /
| 0.51M | 0.00144$\pm$0.000462 | 0.0117$\pm$0.00794 | /
$M=0.772$ | 2.04M | 0.000743$\pm$ 0.000168 | 0.00451$\pm$0.00394 | /
| 8.14M | 0.000511$\pm 6.91\times 10^{-5}$ |  0.00305$\pm$0.00117 | /
| 0.51M |  0.00545$\pm$0.000536 |  0.0167$\pm$0.00465 |  0.0325$\pm$0.00978
$M=0.84$ | 2.04M | 0.00578$\pm$0.000478 | 0.0229$\pm$0.00663 | 0.0470$\pm$0.0134
| 8.14M | 0.00629 $\pm$0.000668 | 0.0304$\pm$0.00386 | 0.0609$\pm$0.00619
| 0.51M | 0.00369$\pm$0.000774 |  0.0186$\pm$0.00686 | /
$M=0.71$ | 2.04M |  0.00355$\pm$ 0.000435 | 0.0191$\pm$0.00586 | /
| 8.14M | 0.00371$\pm$0.000686 | 0.0221$\pm$0.00818 | /
| 0.51M | 0.0412$\pm$0.00564 | 0.0575$\pm$0.0202 | /
$M=0.905$ | 2.04M | 0.0389$\pm$0.00239 | 0.0432$\pm$0.00568 | /
| 8.14M |  0.0374$\pm$0.00255 |  0.0411$\pm$0.00445 | /
Table 2: Accumulated $L_{2}$ errors for $N$ prediction steps.

Due to the different random seeds used for initializing the training of neural
networks, the resulting models exhibit certain variations in predictions. To
address this, we trained five models and calculated the average prediction
error at a given prediction step, as shown in figure 4. The shade areas
represent the $\pm 0.5$ standard deviation. Figure 4(a) gives the loss history
over 999 prediction steps. At free-stream Mach numbers of 0.71, 0.755, and
0.772, since there are no complex shock or compression wave motions, the flow
structure is primarily characterized by periodic vortex shedding.
Consequently, although the error accumulates with the increase in prediction
steps, the overall errors remain small. Notably, the case at Mach 0.71 is an
extrapolation test, which explains its larger error compared to the other two.
We observe a notable difference in error behavior at Mach number 0.905
compared to the other three. Here, the error becomes significant early on but
does not display a clear growth trend. This phenomenon can be attributed to
the considerable deviation of the flow regime from the majority of training
samples: a stable shock wave system with moderate interaction between shock
waves and vortices. While phase deviations are present, their cumulative
effect over time remains relatively modest.

At Mach number 0.84, the moving shock wave over the aerofoil surface exhibits
significant unsteadiness and distortions as it interacts with the boundary
layers. As a result, the errors observed are substantially higher than those
in other scenarios shown in figure 4(a). We are particularly intrigued by this
case, and thus we have conducted extended prediction steps (2499 steps,
equivalent to 9996 simulation steps). As shown in figure 4(b), despite
displaying an upward trend, the prediction remains stable.

As shown in table 2, we calculate the mean $L_{2}$ errors and standard
deviations for the first $N$ prediction steps, where $N=99$, $N=999$ and
$N=2499$ (only for $M_{\infty}=0.84$), respectively. The bold texts show the
highest error or standard deviations. The error values and standard deviations
further illustrate that in complex physical problems such as these, it is
insufficient to rely solely on the interpolation and extrapolation
expectations tested by the simple categorization of data labels (i.e. Mach
numbers). One must also consider whether the statistical distribution of the
test data is similar to that of the training dataset. For instance, at
$M_{\infty}=0.84$, although categorized as an interpolation case, the
structure and evolution of the flow field are sufficiently intricate that they
have deviated from the scope of the samples within the training dataset,
resulting in a pronounced error. Upon comparing the performance of three
different sizes of neural networks, it is evident that the model “8.14M ”
exhibits the most optimal performance overall. While the model “8.14M ” may
not consistently outperform others in individual tests, this observation also
implies that the marginal benefit of further increasing the neural network
size is limited for the current task. Subsequent analyses will delve deeper
into the results obtained with the model “8.14M”.

Now, let’s evaluate the accuracy and errors in flowfields predicted by the
neural network model. We will mainly focus on two representative cases: the
shock-free case at $M_{\infty}=0.755$ and the shock wave oscillation case at
$M_{\infty}=0.84$. Figure 5 compares the flow field at $M_{\infty}=0.755$
after 999 prediction steps with CFD data, showing density, two velocity
components, and pressure from top to bottom, respectively. Although there are
visual differences, one can still observe $O(0.01)$ deviations in the wake
region. Figure 6 illustrates the flow field at $M_{\infty}=0.84$ after 999
prediction steps, compared with CFD results. The velocity and pressure contour
plots demonstrate a noticeable phase shift in the motion of shockwaves along
the lower side of the aerofoil. Despite minor errors, particularly near shock
waves and areas of strong compression, as well as in the wake, the overall
prediction of the flow field structure is remarkably accurate. For an extended
inference test, we also demonstrate that the flow field remains stable even
after 2499 steps (see appendix figure A2).

####  4.1.2 Aerodynamic characteristics

We discover that although the transient flow field experiences deviations due
to the accumulation of errors and phase shifts, the learned model still
effectively captures the statistical quantities and frequency components of
the flowfield. To assess the aerodynamic characteristics, we compare the
predicted results with CFD data in terms of the time-averaged pressure
coefficient and its fluctuations. Additionally, we analyze the FFT spectrum of
the pressure signal measured at a point close to the aerofoil trailing edge,
specifically at $x/c=0.972$, $y/c=-0.00452$.

Figure 7 (a) shows the distribution of the pressure coefficient averaged over
time on the surface of the airfoil, where the shading indicates the range of
$0.5\sigma$ ($\sigma$, standard deviation) above and below the average value,
showing the uncertainty in the prediction of learned models using different
initialization weights. First, let us consider the interpolation tests at Mach
numbers 0.755, 0.772, and 0.84. Figures 7 (a) and 7 (b) show that the
predictions are in excellent agreement with the CFD data. The shade areas
indicating the uncertainty caused by random seeds are not pronounced. In
figure 7(b), the higher values close to $x/c=1.0$ are caused by a series of
compression waves that develop in the region near the trailing edge. Note that
at Mach number 0.84, the pressure fluctuation shows a peak at around
$x/c=0.45$, corresponding to the mean location of the strong shock.
Furthermore, the observation that the pressure fluctuation is higher than in
the other two scenarios corresponds well with the complex wave systems and
strong unsteady processes inherent in the flow field, as depicted in Figure 6.
The peak value around $x/c=0.4$ is attributed to the formation of a strong
shock wave.

Let’s further look at the extrapolation tests. As illustrated in Figure 7(c),
the time-averaged pressure distributions for Mach numbers of 0.71 and 0.905
exhibit a satisfactory correlation with the CFD data. Figure 7(d) shows the
pressure fluctuations, revealing that both cases exhibit weaker unsteadiness
compared to others. At Mach 0.71, vortex shedding dominates within the shock-
free flow regime, while at Mach 0.905, stabilized shock waves prevail. In both
cases, the neural network prediction exhibits a high degree of accuracy,
closely aligning with the CFD data.

Figure 8 shows the FFT spectral of the pressure signals at the probe point
($x/c=0.972$, $y/c=-0.00452$). Each is evaluated using five models, and the
shad shows $\pm 0.5\sigma$. The FFT spectral profiles reveal prominent
frequencies occurring at approximately $fc/U_{\infty}=2.30$ for Mach number
0.755, $fc/U_{\infty}=2.20$ for Mach number 0.772, and $fc/U_{\infty}=1.77$
for Mach number 0.84. These frequencies are indicative of vortex shedding in
the wake. Harmonics of these dominant frequencies manifest in the high-
frequency spectrum. Notably, at Mach number 0.84, a shock wave is observed,
oscillating alternately on the upper and lower aerofoil surfaces at a lower
frequency. This behavior aligns with the occurrence of the peak value
$fc/U_{\infty}=0.0733$ and its harmonics in the FFT profile.

Figure 9 shows the FFT spectral of pressure signals for the extrapolation
results. The test at Mach number 0.71 appears to be quite successful in terms
of FFT spectrum predictions. The primary frequency and its harmonics closely
resemble those observed in other cases without shock waves. The flow regime at
Mach number 0.905, as mentioned previously, exhibits strong and stable shock
waves, and the primary frequency values in the FFT spectrum remain very close
to those in the CFD data. In conclusion, we can see that the neural network
effectively replicates the frequency characteristics of the physical process.

Figure 5: The flowfield from the inference at $M_{\infty}=0.755$ after 999
prediction steps. Figure 6: The flowfield from the inference at
$M_{\infty}=0.84$ after 999 prediction steps.

(a)

(b)

(c)

(d)

Figure 7: Pressure coefficients and fluctuations: (a-b) represent results from
interpolation tests, while (c-d) depict outcomes from extrapolation tests. The
inference is initialized with an unseen instantaneous flowfield. Figure 8:
Spectrum of the pressure signals using FFT. From top to bottom:
$M_{\infty}=0.755$, $M_{\infty}=0.772$ and $M_{\infty}=0.84$. Here, “NN”
represents the neural network, while “CFD” denotes the ground truth, i.e. CFD
data. Figure 9: Spectrum of the pressure signals using FFT. From top to
bottom: $M_{\infty}=0.71$ and $M_{\infty}=0.905$. Here, “NN” represents the
neural network, while “CFD” denotes the ground truth, i.e. CFD data.

###  4.2 Inference from a time-averaged flowfield

Moreover, we undertake a more demanding test for the trained model –
predicting the evolution from the time-averaged flowfield. While this scenario
serves as an exceptionally tough and rigorous ’out-of-distribution’ test,
seldom encountered in other fields, its significance is pronounced in the
realm of fluid dynamics.

We conduct inference using the learned model based on the time-averaged flow
field for 5000 prediction steps, equivalent to 20000 simulation steps in the
CFD simulations. Figure 10 illustrates three interpolation scenarios at Mach
numbers 0.755, 0.772, and 0.84. In these scenarios, the learned model predicts
the complete evolution of the flow field, showcasing the gradual development
of flow instability from the initial time-averaged flow field to a fully
developed, periodic unsteady flow field. Moreover, Figure 11 portrays two
extrapolation scenarios corresponding to Mach numbers 0.71 and 0.905. These
extrapolation cases pose significant challenges to the current methodology.
Specifically, the model struggles to capture the unsteady evolution of the
flow in the case of 0.71. In the case of 0.905, while the level of
unsteadiness is minimal, the fully developed flow field structure remains
acceptable, providing a promising outcome.

(a) $M_{\infty}=0.755$

(b) $M_{\infty}=0.772$

(c) $M_{\infty}=0.84$

Figure 10: The predicted evolution of the flow initialized using time-mean
data at $M_{\infty}=0.755$, $M_{\infty}=0.772$ and $M_{\infty}=0.84$
(interpolation tests). From left to right, the Mach number contour maps are
shown after 1, 6, 11, 51, 101, and 301 prediction steps.

(a) $M_{\infty}=0.71$

(b) $M_{\infty}=0.905$

Figure 11: The predicted evolution of the flow was initialized using time-mean
data at $M_{\infty}=0.71$ and $M_{\infty}=0.905$ (extrapolation tests). From
left to right, the Mach number contour maps are shown after 1, 6, 11, 51, 101,
and 301 prediction steps.

Having removed the transient stage data, we calculate the statistics,
including the time-averaged pressure coefficient and fluctuations. Similar to
the previous analysis, the FFT spectra of the pressure signals are measured at
$x/c=0.972$, $y/c=-0.00452$, close to the aerofoil trailing edge.

Figure 12 presents the distribution of average surface pressure coefficients
and pressure fluctuation, with $\pm 0.5\sigma$ indicating the uncertainty
intervals obtained from five models. From figure 12(a), it is evident that the
distribution of average pressure coefficients closely aligns with the CFD
data, demonstrating a close and highly robust agreement (narrow uncertainty
regions), although the pressure fluctuation distribution exhibits larger areas
of uncertainty as shown in figure 12(b).

The inference from the time-mean at extrapolation scenarios presents greater
complexity. However, there are notable successes. As seen in Figure 12(c), the
distribution of average pressure coefficients exhibits robust performance for
both Mach numbers. Challenges arise at Mach number 0.71, as it appears that
the pressure fluctuations on the wall have not been adequately stimulated, as
shown in Figure 12(d).

As depicted in figure 13, the learned model demonstrates remarkable success in
accurately capturing the peaks and harmonics in the FFT spectra at Mach
numbers of 0.755, 0.772, and 0.84 when predicting from the time-averaged flow
field. Although there is a slightly larger uncertainty region compared to the
preceding tests, overall performance remains commendable. In the extrapolation
experiments illustrated in figure 14, while the FFT spectra distribution at
Mach number 0.905 may exhibit some deviations from expectations, the model’s
performance remains relatively strong. Notably, at Mach number 0.71, where the
absence of shock waves leads to high-frequency vortex shedding dominating low-
frequency motion, the model initially appears to replicate this pattern
accurately, akin to the previous tests. However, as observed earlier, the
pressure fluctuations on the surface are not fully incited, resulting in an
overall lower distribution of FFT signal compared to the CFD data.

In summary, we would like to re-emphasize that this test is exceptionally
challenging, as it deviates significantly from conventional evaluations within
the field of deep learning. Unlike typical tests that involve providing a
transient flow field as an initial condition, this test begins with a time-
averaged field and relies solely on a neural network’s learned ”flow evolution
patterns” to predict subsequent developments and instabilities. While the
model’s performance appears suboptimal in extrapolation scenarios, it
demonstrates excellent performance in interpolation cases. In these instances,
the model accurately infers complete unsteady evolutions and successfully
reproduces statistical quantities such as average pressure distribution and
pressure fluctuation values.

(a)

(b)

(c)

(d)

Figure 12: Pressure coefficients and fluctuations: (a-b) represent results
from interpolation tests, while (c-d) depict outcomes from extrapolation
tests. The inference is initialized with the time-mean data. Figure 13:
Spectrum of the pressure signals using FFT. From top to bottom:
$M_{\infty}=0.755$, $M_{\infty}=0.772$ and $M_{\infty}=0.84$. The inference is
initialized using the time-mean data. Here, “NN” represents the neural
network, while “CFD” denotes the ground truth, i.e. CFD data. Figure 14:
Spectrum of the pressure signals using FFT. From top to bottom:
$M_{\infty}=0.71$ and $M_{\infty}=0.905$. The inference is initialized using
the time-mean data. Here, “NN” represents the neural network, while “CFD”
denotes the ground truth, i.e. CFD data.

##  5 Training time and performance

The performance of trained DNN models is one of the central factors motivating
their use, in particular, in aerodynamic shape optimizations. We evaluate our
models using a regular workstation with 12 cores, i.e. Intel® Xeon® CPU
E5-1650 v4 @ 3.60GHz, with an NVidia GeForce GTX 1080 Ti GPU. Due to the
difference in implementations, we compare solvers in terms of wall clock time,
which is a pure computation time and is averaged over multiple runs without
taking the start-up and initialization overheads into account. A typical two-
dimensional laminar simulation with CFL3D solver with the grid resolution
$1024\times 128$ requires 4000 time steps with a wall clock time of 21 hours
using 8 CPU processes. Similarly, we assess the elapsed wall time required for
1000 predictions on a single batch of data, which equates to 4000 CFD steps,
using neural network models of sizes 0.5M, 2.0M, and 8.1M. The respective wall
times are approximately 7.01 seconds, 7.20 seconds, and 7.15 seconds. It is
noteworthy that the runtime per solution can be notably reduced when
evaluating multiple batch solutions simultaneously. Overall, The neural
network inference has achieved considerable acceleration compared to CFD
simulations.

##  6 Learned-model-based global instability analysis

###  6.1 General description

In the previous sections, it was noted that the learned model not only
accurately forecasts the evolution of intricate flow fields under unseen
conditions but also performs exceptionally well in more stringent tests, where
the development of unsteady flow states stems from a time-averaged flow field.
Encouraged by these test outcomes, we are motivated to explore a fundamental
question: Can a model be trained such that it uncovers the underlying
mechanisms of a dynamical system? If confirmed, the trained model de facto
becomes a differentiable operator, encapsulating the dynamic characteristics
present in the training dataset. To explore this question, this section
utilizes the trained model as an operator for classical global stability
analysis.

The choice of the base flow is crucial to any comprehensive analysis of global
stability. The base flow serves as the reference point for linearizing the
dynamics. However, the flow field is predominantly governed by instabilities
for laminar and scale-resolving simulations, lacking a true stable flow
pattern suitable for linearization. Nonetheless, as the vortex shedding and
self-sustained shockwave (when buffet occurs) are periodic, we hence assume
that the flow has statistically stationary properties over a long-time
average. Note that this long-time-averaged flowfield is neither a solution of
the governing equations nor a stationary point of the learned recursive
operator, The following global linear stability analysis makes sense in
identifying the frequency content of large coherent structures rather than
addressing the stability of the base flow [19].

###  6.2 Discrete and continuous dynamical systems

Because our neural network model is a discrete dynamical system (as in
equation 1), and stability analysis is typically based on continuous dynamical
systems, we first need to establish a connection between discrete and
continuous dynamical systems. This will enable us to then solve the
corresponding eigenvalue problem.

Firstly, we revisit the neural network equation 1 and transform it into 4,
introducing the notation $\mathcal{N}$ under the assumption that $\theta$ and
$M_{\infty}$ are predetermined:

$f_{\theta}(\mathbf{q}_{t},M_{\infty}):=\mathcal{N}(\mathbf{q}_{t})$

To clarify the discussion, we assume the time step $\Delta t$ is uniform and
denote $\mathbf{q}_{t}$ as $\mathbf{q}_{n}$ and $\mathbf{q}_{t+\Delta t}$ as
$\mathbf{q}_{n+1}$. Then we have

$\displaystyle\mathbf{q}_{n+1}$ $\displaystyle=\mathcal{N}(\mathbf{q}_{n})$
(4) $\displaystyle\bar{\mathbf{q}}+\mathbf{q}^{\prime}_{n+1}$
$\displaystyle=\mathcal{N}(\bar{\mathbf{q}}+\mathbf{q}^{\prime}_{n})$ (5)
$\displaystyle\bar{\mathbf{q}}+\mathbf{q}^{\prime}_{n+1}$
$\displaystyle=\mathcal{N}(\bar{\mathbf{q}})+\frac{\partial\mathcal{N}}{\partial\mathbf{q}}\Bigg{|}_{\bar{\mathbf{q}}}\mathbf{q}^{\prime}_{n}$
(6)

Here, $\bar{\mathbf{q}}$ represents the time-averaged flow variables,
calculated by

$\displaystyle\bar{\mathbf{q}}$
$\displaystyle=\sum_{i=0}^{\infty}\mathbf{q}_{i+1}$
$\displaystyle=\sum_{i=0}^{\infty}{\mathcal{N}(\mathbf{q}_{i})}$
$\displaystyle=\overline{\mathcal{N}(\mathbf{q}_{i})}=\mathcal{N}(\bar{\mathbf{q}})-\mathbf{h}$

where $\mathbf{h}$ represents the forcing term due to the nonlinearity of the
operator. Hence, we have

$\displaystyle\mathbf{q}^{\prime}_{n+1}$
$\displaystyle=\frac{\partial\mathcal{N}}{\partial\mathbf{q}}\Bigg{|}_{\bar{\mathbf{q}}}\mathbf{q}^{\prime}_{n}+\mathbf{f}=\mathbf{L}(\bar{\mathbf{q}})\
\mathbf{q}^{\prime}_{n}+\mathbf{h}$ (7)
$\displaystyle\mathbf{L}(\bar{\mathbf{q}})$
$\displaystyle=\mathbf{V_{L}}\mathbf{\Sigma}\mathbf{V_{L}}^{-1}$ (8)

where $\mathbf{L}$ represents the Jacobian matrix of the operator
$\mathcal{N}$, $\mathbf{V_{L}}$ is the eigenvector of $\mathbf{L}$, and
$\mathbf{\Sigma}$ the eigenvalue.

On the other hand, looking at the continuous dynamical system (i.e. temporally
continuous but discretized in space) is expressed as

$\displaystyle\frac{d\mathbf{q}}{dt}=\mathcal{S}(\mathbf{q})$ (9)

where $\mathcal{S}$ represents the nonlinear operator, such as the Navier-
Stokes equation. Assuming $\mathbf{q}=\bar{\mathbf{q}}+\mathbf{q}^{\prime}$,
the linearized ordinary differential equation in a continuous form can be
written as

$\displaystyle\frac{d\mathbf{q}^{\prime}}{dt}$
$\displaystyle=\frac{\partial\mathcal{S}}{\partial\mathbf{q}}\Bigg{|}_{\bar{\mathbf{q}}}\mathbf{q}^{\prime}+\mathcal{S}(\bar{\mathbf{q}})$
(10)
$\displaystyle=\mathbf{A}\mathbf{q}^{\prime}+\mathcal{S}(\bar{\mathbf{q}})$
(11)

The general solution of the nonhomogeneous equation is the sum of the
homogeneous solution and the particular solution, and the discretized form can
be expressed as

$\displaystyle\mathbf{q}^{\prime}_{n+1}$
$\displaystyle=\mathrm{exp}({\mathbf{A}\Delta
t})\mathbf{q}^{\prime}_{n}+\mathbf{f}$ (12)
$\displaystyle=\mathbf{V_{A}}\mathrm{exp}(\mathbf{\Lambda}\Delta
t)\mathbf{V_{A}}^{-1}\mathbf{q}^{\prime}_{n}+\mathbf{f}$ (13)
$\displaystyle\mathbf{A}$
$\displaystyle=\mathbf{V_{A}}\mathbf{\Lambda}\mathbf{V_{A}}^{-1}$ (14)

Here, $\mathbf{A}$ represents the Jacobian matrix of the continuous dynamical
system, $\mathbf{V_{A}}$ is the eigenvector of $\mathbf{A}$, and
$\mathbf{\Lambda}$ the eigenvalue. The similarity of the two discrete
dynamical system recurrence formulas, equations 7 and 12, means that If a
learned model trained on high-fidelity data achieves sufficient accuracy, it
can provide a high degree of validity, potentially leading to
$\mathbf{L}=\mathbf{A}$ in the limit. This equivalence is contingent upon the
data-driven model accurately capturing the dynamics of the system.

Hence, from equations 8 and 14, we have

$\displaystyle\lambda_{i}=\mathrm{log}(\sigma_{i})/\Delta t$ (15)

where $\lambda_{i}$ are the eigenvalues of $\mathbf{A}$ and $\sigma_{i}$ the
eigenvalues of $\mathbf{L}$. Under the assumption of a highly accurate data-
driven model, the eigenvectors $\mathbf{V_{L}}$ can be regarded as well-
estimated counterparts to $\mathbf{V_{A}}$. This implies that, as long as we
can train a sufficiently accurate neural network, we can analyze the dynamics
of the original continuous system by only examining the dynamical
characteristics of the trained model. This is possible without relying on the
governing equations of a continuous dynamical system such as equation 12.
Instead, insights about the original dynamical system can be gained by
analyzing $\mathbf{L}$.

The eigenvalues of $\mathbf{A}$ are associated with the system’s natural
frequencies. The real parts of the eigenvalues indicate the exponential growth
or decay rates, while the imaginary parts correspond to oscillatory behavior.
Nonetheless, learned models usually come with uncertainty, which is typically
caused by random initialization of parameters and mini-batch shuffling. To
obtain a reliable Jacobian matrix $\mathbf{A}$, we employ “model ensembling”
by combining five trained models which are initialized with different random
seeds, and calculate the ensemble average $\mathbf{A}$. The ensemble approach
captures a broader range of multiple solutions, and reduces sensitivity to
specific weight initializations, resulting in a more stable and representative
average gradient. This approach has been widely used in applications such as
learned model-based control and inverse design tasks [38, 27].

###  6.3 Eigenvalue spectrum and dominant modes

Firstly, we look at the eigenvalue spectrum plot at $M_{\infty}=0.772$. In
figure 6.3(a), the eigenvalue spectrum exhibits a large cluster of points in
the region with negative real parts, and several eigenvalues extend from the
stable half plane to the unstable half plane, indicating that the system is in
a critical stability state. In such cases, the system might undergo a
transition from stability to instability or exhibit oscillatory behavior.

Temporal frequencies of the modes, and the relative errors between the modes
and the frequencies identified in the FFT (see figure 8) are given by table 3
The dominant peak in the FFT corresponds to unstable modes $m_{0}$ to $m_{8}$,
with a relative frequency error of above 12%, suggesting that significant
nonlinear effects play a role in selecting the frequency.

Mode | Eigenvalue $\frac{\lambda}{2\pi}$ | FFT frequency | $\Delta\mathrm{Im}\\{\frac{\lambda}{2\pi}\\}$
---|---|---|---
$m_{0}$ | $1.1703308+2.81216j$ | 2.20 | 27.82%
$m_{2}$ | $0.9641248+2.6544704j$ | 2.20 | 20.66%
$m_{4}$ | $0.6119804+2.647502j$ | 2.20 | 20.3%
$m_{8}$ | $0.034572423+2.483544j$ | 2.20 | 12.89%
Table 3: Selected frequencies for $M_{\infty}=0.772$ Mode | Eigenvalue $\frac{\lambda}{2\pi}$ | FFT frequency | $\Delta\mathrm{Im}\\{\frac{\lambda}{2\pi}\\}$
---|---|---|---
$m_{12}$ | $2.8078394+0.27758664j$ |  |
$m_{70}$ | $-0.06118219+\mathbf{0.09932874}j$ | 0.077 | 28.57%
$m_{72}$ | $-0.06684557+0.3259675j$ |  |
$m_{74}$ | $-0.07541487+\mathbf{1.9026724}j$ | 1.77 | 7.34%
Table 4: Selected frequencies for $M_{\infty}=0.84$

In figure 6.3(b), the pressure modal patterns of modes $m_{0}$ to $m_{8}$
exhibit a clear alternation between positive and negative values, indicating
the presence of pressure wave propagation. Particularly noteworthy is the
discernible propagation from the airfoil trailing edge towards the far-field
direction. Notably, within the wake region downstream, the pattern arrangement
suggests oscillations in pressure associated with vortex shedding. This
interpretation finds support in the stream-wise velocity modes (i.e. the real
part of $\hat{u}$) depicted in figure 6.3(c), where it becomes apparent that
the dominant mode driving instability within the flow system corresponds to
vortex shedding.

Upon examining the case of Mach number 0.84, it is observed that the
eigenvalue spectrum exhibits multiple branches extending towards the unstable
upper half-plane. We particularly focus on modes 12, 70, 72, and 74. As
depicted in figure 6.3(a), the branches corresponding to modes 12 and 74
represent high-frequency signals (i.e., the imaginary parts of the
eigenvalues). The associated pressure modal distributions in figure 6.3(b) and
velocity modal distributions in figure 6.3(c) resemble those observed at Mach
number 0.772, indicating the occurrence of vortex shedding resulting from
Kelvin-Helmholtz instability. Furthermore, the alternating arrangement of
positive and negative pressure modal patterns (see figure 6.3(b)) reveals the
direction of pressure wave propagation from the airfoil trailing edge towards
the far field.

Unlike Mach number 0.772, the case at Mach number 0.84 exhibits additional
low-frequency unstable branches. To investigate this, we examined modes 70 and
72 at the axis $\mathbf{Re}{\lambda}=0$. The frequencies corresponding to
these two modes are 0.077 and 1.77, respectively, likely associated with low-
frequency peaks observed in FFT analysis of the flow field. As analyzed
earlier, these frequencies correspond to the low-frequency motion of shock
waves following buffet onset. The pressure modal distributions in figure
6.3(b) and velocity modal distributions in figure 6.3(c) for modes 70 and 72
reveal stripe structures parallel to the flow direction. In the trailing edge
region, these modes correspond to large-scale spatial structures, consistent
with the low-frequency characteristics of shock wave motion.

In summary, by treating the trained neural networks as differentiable
operators, we conducted a global stability analysis to explore the complex
dynamics embedded within the training dataset. This method diverges from
traditional approaches, as it relies solely on training data rather than
continuous system formulations. Our findings highlight the versatility of
neural networks in capturing system behavior and demonstrate the applicability
of global stability analysis in interpreting the learned dynamics of neural
networks. This approach provides valuable insights into the underlying
mechanisms of complex systems and contributes to the growing intersection of
machine learning and classical stability analysis methodologies.

#

(a) Eigenvalue spectrum

(b) Pressure modes, real part.

(c) Stream-wise modes, real part

Figure 15: $M_{\infty}=0.772$

#

(a) Eigenvalue spectrum

(b) Pressure modes, real part.

(c) Stream-wise modes, real part

Figure 16: $M_{\infty}=0.84$

##  7 Concluding remarks

In this study, we propose a novel approach for data-driven prediction of
transonic flow fields. We trained attention U-Net models using roll-out
training and noise addition techniques. The learned models can predict the
evolution of complex flows under conditions that were not presented during
training, including initial instantaneous flow fields at unseen Mach numbers,
and even with time-averaged flow fields. We demonstrate that the neural
network model not only qualitatively predicts the subsequent evolution of the
flow field but also quantitatively forecasts the aerodynamic characteristics
of the airfoil, including average pressure distribution, fluctuating pressure
distribution, and the frequency spectrum.

We suggest that a successful learned model should be able to replicate the
characteristics of dynamical systems. We have found, through classic Jacobian
matrix-based global stability analysis, that the learned model can obtain
unstable eigenvalues and eigenvectors, and the mode shapes and frequencies
were found to match very well with the FFT analysis of the pressure signals.
In the shock-free case at Mach number 0.772, the eigenvalue spectrum plot
shows a high-frequency unstable branch, corresponding to the Kelvin-Helmholtz
mode. In the buffet-onset case at Mach number 0.84, the eigenvalue spectrum
exhibits multiple unstable branches, associated with the high-frequency
Kelvin-Helmholtz mode and the low-frequency shock wave motion.

In general, our work contributes to the broader goal of enhancing the
interpretability of neural network models, as it offers a systematic
methodology for analyzing and interpreting the behavior of such models in the
context of complex dynamical systems. The integration of deep learning and
global instability analysis stands as a promising methodology for addressing
the complexities of flow control and shape optimization. The neural-network-
based resolvent analysis will be an interesting direction.

## Code and data availability

To ensure reproducibility, the full code, and training as well as test
datasets for this work will be made available upon publication.

## Acknowledgements

This work was supported by the European Research Council Consolidator Grant
_SpaTe_ (Grant agreement ID: 863850).

## Declaration of interests

The authors report no conflict of interest.

## Appendix A

Figure A1 shows an overview of the attention U-Net architecture. It comprises
encoder and decoder pathways: the encoder progressively extracts hierarchical
features from input images, while attention gates within the architecture
selectively emphasize relevant features and suppress irrelevant ones. The
decoder pathway then upsamples the features to recover spatial information,
with attention mechanisms refining the reconstructed features. During the
forward pass, the input flowfield undergoes convolutional operations and
attention-based refinement to produce outputs, allowing the model to focus on
crucial input regions, thus enhancing its performance in the prediction. To
explore the influence of neural network size on accuracy, multiple training
runs are conducted using three distinct network sizes (see table 2). This is
accomplished by adjusting the number of features in each layer.

Figure A2 illustrates the flow field at $M_{\infty}=0.84$ after 2499
prediction steps, compared with CFD results. The velocity and pressure contour
plots reveal minor errors near shock waves and areas of strong compression, as
well as in the wake. However, the overall prediction of the flow field
structure is remarkably accurate, with the flow field remaining remarkably
stable even after 2499 steps.

(a) The attention U-Net.

(b) The basic schematic of attention gate.

Figure A1: An overview of the attention U-Net architecture with input and
output specifications. It contains 8.14 million weights. Figure A2: The
flowfield from the inference at $M_{\infty}=0.84$ after 2499 prediction steps.

## References

  * [1] G Y Nieuwland and B M Spee.  Transonic airfoils: Recent developments in theory, experiment, and design.  Annual Review of Fluid Mechanics, 5(1):119–150, 1973.
  * [2] H Tijdeman and R Seebass.  Transonic flow past oscillating airfoils.  Annual Review of Fluid Mechanics, 12(1):181–222, 1980.
  * [3] B.H.K. Lee.  Self-sustained shock oscillations on airfoils at transonic speeds.  Progress in Aerospace Sciences, 37(2):147–196, 2001.
  * [4] B. H. K. Lee.  Oscillatory shock motion caused by transonic shock boundary-layer interaction.  AIAA Journal, 28(5):942–944, 1990.
  * [5] Herbert Olivier, Thomas Reichel, and Mitja Zechner.  Airfoil flow visualization and pressure measurements in high-reynolds-number transonic flow.  AIAA Journal, 41(8):1405–1412, 2003.
  * [6] Q. Xiao, H. M. Tsai, and F. Liu.  A numerical study of transonic buffet on a supercritical airfoil.  AIAA Paper 2004-1056, 2004.
  * [7] Sébastien Deck.  Numerical simulation of transonic buffet over a supercritical airfoil.  AIAA Journal, 43(7):1556–1566, 2005.
  * [8] Li-Wei Chen, Chang-Yue Xu, and Xi-Yun Lu.  Numerical investigation of the compressible flow past an aerofoil.  J. Fluid Mech., 643:97–126, 2010.
  * [9] J. D. Crouch, A. Garbaruk, D. Magidov, and A. Travin.  Origin of transonic buffet on aerofoils.  Journal of Fluid Mechanics, 628:357–369, 2009.
  * [10] Fulvio Sartor, Clément Mettot, and Denis Sipp.  Stability, receptivity, and sensitivity analyses of buffeting transonic flow over a profile.  AIAA Journal, 53(7):1980–1993, 2015.
  * [11] Yoimi Kojima, Chi-An Yeh, Kunihiko Taira, and Masaharu Kameda.  Resolvent analysis on the origin of two-dimensional transonic buffet.  Journal of Fluid Mechanics, 885:R1, 2020.
  * [12] A. Bouhadji and M. Braza.  Organised modes and shock–vortex interaction in unsteady viscous transonic flows around an aerofoil: Part i: Mach number effect.  Computers & Fluids, 32(9):1233–1260, 2003.
  * [13] A Bouhadji and M Braza.  Organised modes and shock–vortex interaction in unsteady viscous transonic flows around an aerofoil: Part ii: Reynolds number effect.  Computers & Fluids, 32(9):1261–1281, 2003.
  * [14] Pradeep Moise, Markus Zauner, and Neil D. Sandham.  Large-eddy simulations and modal reconstruction of laminar transonic buffet.  Journal of Fluid Mechanics, 944:A16, 2022.
  * [15] Jiaqing Kou, Weiwei Zhang, and Chuanqiang Gao.  Modal analysis of transonic buffet based on POD and DMD techniques (in chinese).  Acta Aeronautica et Astronautica Sinica, 37(9):2679–2689, 2016.
  * [16] Andre Weiner and Richard Semaan.  Robust dynamic mode decomposition methodology for an airfoil undergoing transonic shock buffet.  AIAA Journal, 61(10):4456–4467, 2023.
  * [17] Yuta Iwatani, Hiroyuki Asada, and Soshi Kawai.  POD, DMD, and resolvent analysis of transonic airfoil buffet.  2022\.
  * [18] Yuta Iwatani, Hiroyuki Asada, Chi-An Yeh, Kunihiko Taira, and Soshi Kawai.  Identifying the self-sustaining mechanisms of transonic airfoil buffet with resolvent analysis.  AIAA Journal, 61(6):2400–2411, 2023.
  * [19] Kunihiko Taira, Steven L. Brunton, Scott T. M. Dawson, Clarence W. Rowley, Tim Colonius, Beverley J. McKeon, Oliver T. Schmidt, Stanislav Gordeyev, Vassilios Theofilis, and Lawrence S. Ukeiley.  Modal analysis of fluid flows: An overview.  AIAA Journal, 55(12):4013–4041, 2017.
  * [20] R. Maulik, O. San, A. Rasheed, and P. Vedula.  Subgrid modelling for two-dimensional turbulence using neural networks.  Journal of Fluid Mechanics, 858:122–144, 2019.
  * [21] K. Um, R. Brand, P. Holl, R. Fei, and N. Thuerey.  Solver-in-the-loop: Learning from differentiable physics to interact with iterative PDE-solvers.  Advances in Neural Information Processing Systems, 2020.
  * [22] J. F. MacArt, J. Sirignano, and J. B. Freund.  Embedded training of neural-network subgrid-scale turbulence models.  Physical Review Fluids, 6(5):050502, 2021.
  * [23] D. Kochkov, J. A. Smith, A. Alieva, Q. Wang, M. P. Brenner, and S. Hoyer.  Machine learning–accelerated computational fluid dynamics.  Proceedings of the National Academy of Sciences, 118(21):e2101784118, 2021.
  * [24] Björn List, Li-Wei Chen, and Nils Thuerey.  Learned turbulence modelling with differentiable fluid solvers: physics-based loss functions and optimisation horizons.  Journal of Fluid Mechanics, 949:A25, 2022.
  * [25] Rui Wang, Karthik Kashinath, Mustafa Mustafa, Adrian Albert, and Rose Yu.  Towards physics-informed deep learning for turbulent flow prediction, 2020.
  * [26] Kimberly Stachenfeld, Drummond B. Fielding, Dmitrii Kochkov, Miles Cranmer, Tobias Pfaff, Jonathan Godwin, Can Cui, Shirley Ho, Peter Battaglia, and Alvaro Sanchez-Gonzalez.  Learned coarse models for efficient turbulence simulation, 2022.
  * [27] Kelsey Allen, Tatiana Lopez-Guevara, Kimberly L Stachenfeld, Alvaro Sanchez Gonzalez, Peter Battaglia, Jessica B Hamrick, and Tobias Pfaff.  Inverse design for fluid-structure interactions using graph network simulators.  In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh, editors, Advances in Neural Information Processing Systems, volume 35, pages 13759–13774. Curran Associates, Inc., 2022.
  * [28] Ziyang Liu, Renkun Han, Miao Zhang, Yang Zhang, Hao Zhou, Gang Wang, and Gang Chen.  An enhanced hybrid deep neural network reduced-order model for transonic buffet flow prediction.  Aerospace Science and Technology, 126:107636, 2022.
  * [29] Georg Kohl, Li-Wei Chen, and Nils Thuerey.  Turbulent flow simulation using autoregressive conditional diffusion models.  arXiv:2309.01745, 2023.
  * [30] Alexandre Mauroy, Aivar Sootla, and Igor Mezić.  Koopman Framework for Global Stability Analysis, pages 35–58.  Springer International Publishing, Cham, 2020.
  * [31] Benjamin Herrmann, Peter J. Baddoo, Richard Semaan, Steven L. Brunton, and Beverley J. McKeon.  Data-driven resolvent analysis.  Journal of Fluid Mechanics, 918:A10, 2021.
  * [32] Christopher Lockwood Rumsey, Robert T Biedron, and James Lee Thomas.  CFL3D: its history and some recent applications.  NASA TM–112861, 1997.
  * [33] Christopher L Rumsey.  CFL3D contribution to the AIAA supersonic shock boundary layer interaction workshop.  NASA TM–2010-216858, 2010.
  * [34] W. J. Gordon and C. A. Hall.  Construction of curvilinear co-ordinate systems and applications to mesh generation.  International Journal for Numerical Methods in Engineering, 7(4):461–477, 1973.
  * [35] Li-Wei Chen and Nils Thuerey.  Towards high-accuracy deep learning inference of compressible flows over aerofoils.  Computers & Fluids, 250:105707, 2023.
  * [36] Ozan Oktay, Jo Schlemper, Loic Le Folgoc, Matthew Lee, Mattias Heinrich, Kazunari Misawa, Kensaku Mori, Steven McDonagh, Nils Y Hammerla, Bernhard Kainz, et al.  Attention u-net: Learning where to look for the pancreas.  arXiv:1804.03999, 2018.
  * [37] Christopher M. Bishop.  Pattern Recognition and Machine Learning.  Springer, 2006.
  * [38] Kurtland Chua, Roberto Calandra, Rowan McAllister, and Sergey Levine.  Deep reinforcement learning in a handful of trials using probabilistic dynamics models, 2018.