# Active learning of data-assimilation closures using graph neural networks

Michele Quattromini, Michele Alessandro Bucci, Stefania Cherubini, Onofrio Semeraro

► To cite this version:

Michele Quattromini, Michele Alessandro Bucci, Stefania Cherubini, Onofrio Semeraro. Active learning of data-assimilation closures using graph neural networks. Theoretical and Computational Fluid Dynamics, 2025, 39, pp.17. 10.1007/s00162-025-00737-1. hal-05379430

HAL Id: hal-05379430

https://hal.science/hal-05379430v1

Submitted on 24 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Active learning of data-assimilation closures using
Graph Neural Networks

Michele Quattromini<sup>1,2*</sup>, Michele Alessandro Bucci<sup>3</sup>,
Stefania Cherubini<sup>1</sup>, Onofrio Semeraro<sup>2</sup>

<sup>1</sup>*Dipartimento di Meccanica, Matematica e Management, Politecnico di Bari, Via Orabona 4, Bari, 70126, Italy.*

<sup>2</sup>LISN-CNRS, Université Paris-Saclay, Street, Orsay, 91440, France.

<sup>3</sup>Department, Safran Tech, Digital Sciences & Technologies, Street, Magny-Les-Hameaux, 78114, France.

*Corresponding author(s). E-mail(s): michele.quattromini@poliba.it;
Contributing authors: michele-alessandro.bucci@safrangroup.com;
stefania.cherubini@poliba.it; onofrio.semeraro@universite-paris-saclay.fr;

Abstract

The spread of machine learning techniques coupled with the availability of high-quality experimental and numerical data has significantly advanced numerous applications in fluid mechanics. Notable among these are the development of data assimilation and closure models for unsteady and turbulent flows employing neural networks (NN). Despite their widespread use, these methods often suffer from overfitting and typically require extensive datasets, particularly when not incorporating physical constraints. This becomes compelling in the context of numerical simulations, where, given the high computational costs, it is crucial to establish learning procedures that are effective even with a limited dataset. Here, we tackle those limitations by developing NN models capable of generalizing over unseen data in low-data limit by: i) incorporating invariances into the NN model using a Graph Neural Networks (GNNs) architecture; and ii) devising an adaptive strategy for the selection of the data utilized in the learning process. GNNs are particularly well-suited for numerical simulations involving unstructured domain discretization and we demonstrate their use by interfacing them with a Finite Elements (FEM) solver for the supervised learning of Reynolds-averaged Navier-Stokes equations. We consider as a test-case the data-assimilation of meanflows past generic bluff bodies, at different Reynolds numbers $50 \le Re \le 150$, characterized by an unsteady dynamics. We show that the GNN models successfully predict the closure term; remarkably, these performances are achieved using a
very limited dataset selected through an active learning process ensuring the gen-
eralization properties of the RANS closure term. The results suggest that GNN
models trained through active learning procedures are a valid alternative to less
flexible techniques such as convolutional NN.

**Keywords:** Graph Neural Network (GNN), surrogate model, Machine Learning for fluid dynamics, Reynolds-averaged Navier-Stokes (RANS) Equations, Active Learning

# 1 Introduction

Machine learning (ML) applications are spreading across the scientific community
in the most diverse fields, leading to the development of novel techniques and algo-
rithms encompassed within the emerging field of scientific machine learning. This
breakthrough is fostered by the massive augmentation of computer-aided activities,
the constant increase of high performance computational power, the recent develop-
ments in the field of deep learning and the large availability of data [24]. In a very
schematic way, for a regression problem, ML techniques enable to identify mappings
between observables of a system (inputs) and quantities of interest (outputs) we aim
to predict by leveraging data; when these analysed data are governed by deterministic
or statistical laws, in principle, these mappings correspond to approximating models.
Alternative applications can be found in classification problems, clustering, control
in a large variety of scientific fields, including fluid mechanics where the number of
contributions combining ML and standard techniques of analysis has been constantly
growing in the last years [8, 9, 22, 41, 52, 57].

These applications found a natural ground into modelling problems, ranging from
the closure of Reynolds-Average Navier-Stokes equations for turbulent flows, to wall
model Large Eddy Simulations [4, 28, 59] or optimization problems combining adjoint-
based techniques and machine learning [37, 53]. Among the numerous authors that
addressed this problem, [30] applied classification methods for identifying regions of
uncertainties where the closure term of the Reynolds-averaged Navier-Stokes (RANS)
might fail; [49] combined data assimilation with neural networks (NN) modelling of
the Reynolds stress using limited observation. Other approaches leverage baseline
models such as the the Spalart-Allmaras closure [47], physics-informed NN (PINNs)
[20], regression methods [45], decision trees [19], ensemble methods [35] or genetic
programming [55, 58]. For a broader overview, we refer to [19] and [6], where the
different levels of approximation are discussed together with a critical take on the
limitations of the approach. In the specific case of an eddy viscosity closure model, the
recent studies by [31] demonstrated the effectiveness of Tensor-Basis Neural Networks
(TBNNs) in learning a General Eddy Viscosity model type [39]. TBNNs capitalize on
the tensor decomposition approach proposed by Pope to account for invariances and
streamline the number of parameters to be learnt. The inductive bias introduced by
this modelling approach restricts TBNN application to nearly homogeneous flows with
high Reynolds numbers, where local effects predominate [11]. On the other hand, at low
Reynolds numbers, we may consider as alternative approach a data assimilation one,
where the closure model corresponds to the control parameter of an adjoint-based loop [21]; without explicitly introducing a tensor structure or Boussinesq approximations, this method is well suited for non-homogeneous flows at lower Reynolds numbers. Here, we take inspiration from these applications and mainly focus on unsteady flows developing past bluff bodies at low Reynolds numbers $50 \le Re \le 150$; we consider RANS as baseline, although alternative choices can be also adopted, such as Euler equations or linearized Navier Stokes equations in resolvent form, where the parameter to be tuned is the dissipation term [36, 38, 43]. With respect of standard approaches, we focus here on a supervised learning strategy where the closure model is identified by inference from available data. In principle, we could identify universal closure models directly from data, having at disposal an infinite amount of them. In practice, in real cases, we may deal with a limited amount of data or few localized measurements; these limitations can impact on the use of methods such as neural networks (NN), where the expressivity and generic structure make them suitable for a large class of models, but prone to generalisation problems. In this sense, ML models risk to be representative solely of the datasets included in the training process; thus, it becomes compelling given the available data to circumvent these problems by inputting well selected data during the training or providing prior knowledge through modelling [10, 46].

With this in mind, we will test if it is possible to identify generic closure models from data defined on unstructured meshes while only relying on a small amount of data chosen on principled criteria. In order to answer these questions, a first ingredient is the introduction of graph neural networks (GNN) [27]; this architecture is characterized by complex multi-connected nets of nodes that can be naturally adapted to unstructured meshes: the convolution in a GNN is performed by aggregating information from neighbouring nodes, thus overcoming the limitations imposed by the geometry in contrast with convolutional NN (CNN). Moreover, GNNs: i) show remarkable generalization capabilities as compared to standard network models [44]; ii) are differentiable; iii) provide the possibility of directly targeting the learning of the operator via discrete stencils [see 46]. Due to these features, this architecture has recently attracted attention in fluid mechanics. A review is available on the subject authored by [32], while examples are given by the works of [51] or [18], where wall shear stress are modelled for LES simulations based on GNN. Here, we take inspiration from [17], where a GNN-based architecture incorporating permutation and translation invariance is combined with the statistical solver problem; [17] proofed that the architecture – referred to as deep statistical solver – has some universal approximation properties, and it is capable of operator learning.

In this contribution, GNN are combined with numerical simulations performed using finite elements method (FEM). The GNN-FEM interface allows the use of NN predictions in post-processing FEM analysis since it provides a two-way coupling between NN and FEM environments. However – as already mentioned – the volume of data can represent a bottleneck if the available amount is insufficient. Moreover, also quality of data impacts on the prediction properties in pure data-driven modelling [10], in particular when data at hand are not sufficient in representing correctly the distribution of the overall dataset. Inconsistent or unbalanced data distributions may lead the model to being trained on subsets of data that do not adequately represent the
underneath physics of the problem. Thus, as second, fundamental ingredient we adopt active learning, aimed to increase the dissimilarity between data points and ensure that the model distills all discriminant features necessary to perform the required task effectively. At the best of our knowledge, this application is among the first in fluid mechanics where selection of data is performed by applying an active learning criterion. Readers are directed to [42] for a comprehensive overview of active learning methods in deep learning. We follow the work by [13] where scalar products of gradients associated with the update of the model weights of the GNN are considered as similarity metric between samples; the chosen criterion allows to dynamically increase the training dataset by introducing data that promote diversity in the dataset.

The remainder of the paper is organised as follows. We state the fluid mechanics problem in Sec. 2 and Sec. 3, where we discuss the set of equations used for the numerical simulations, the flow cases and the computational setup (Sec. 3.3). The theoretical background is complemented by a discussion of the data-driven modelling based on GNN in Sec. 4, and how this architecture is coupled with numerical simulations resolved on unstructured meshes. The second part of the paper focuses on the results (Sec. 5): GNN models are trained on datasets composed by the mean flows (input) and the Reynolds stress (target); the baseline is provided by the closure problem in cylinder flows at Reynolds numbers $50 \le Re \le 150$ (Sec. 5.2). Generalization properties are assessed using 4 different benchmarks defined in Sec. 5.1, not included in the training sets. This baseline is thus compared with models trained with different strategies of data selection. First, we consider cases where data augmentation is performed using simulations of flows past bluff bodies of random geometry (Sec. 5.3); second, active learning is introduced in Sec. 5.3.3. The paper finalises with conclusions in Sec. 6.

## 2 Governing equations

In this work, we will focus on incompressible two-dimensional (2D) fluid flows developing past bluff bodies in unsteady flow regimes. We consider the Navier-Stokes (NS) equations for incompressible flows. By denoting the spatial Cartesian coordinates as $\mathbf{x} = (x, y)$, the dynamics of the velocity field $\mathbf{u}(\mathbf{x}, t)$ and pressure field $p(\mathbf{x}, t)$ are governed by

$$ \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla) \mathbf{u} = -\nabla p + \frac{1}{Re} \nabla^2 \mathbf{u} \qquad (1a) $$

$$ \nabla \cdot \mathbf{u} = 0. \qquad (1b) $$

Eq.s 1 are made dimensionless using the characteristic length scale $D$ (e.g. the cylinder diameter), the velocity of the incoming uniform flow $U_\infty$ as the reference velocity and $\rho U_\infty^2$ as the reference pressure. Based on these quantities, the Reynolds number is defined as $Re = U_\infty D/\nu$, where $\nu$ is the kinematic viscosity.

As baseline equations for the data assimilation process, we will consider the Reynolds-averaged Navier-Stokes (RANS) equations, the time averaged formulation of the NS Equations [21]. Notice that, although typically used for turbulent flows, the
RANS approach is well suited also for unsteady flow, for instance for the determination of unsteady periodic flow dynamics through weakly nonlinear analyses [33, 48]. Thus, we introduce the Reynolds decomposition

$$ \mathbf{u}(\mathbf{x}, t) = \bar{\mathbf{u}}(\mathbf{x}) + \mathbf{u}'(\mathbf{x}, t), \qquad (2) $$

where the vector field $\mathbf{u} = (u, v)^T$ is rewritten as sum of the time-averaged velocity vector field $\bar{\mathbf{u}} = (\bar{u}, \bar{v})^T$ and the fluctuation velocity vector field $\mathbf{u}' = (u', v')^T$. Formally, any unsteady flow can be described through this decomposition. Plugging the Reynolds decomposition in Eq. 2 into the Navier-Stokes equations (Eq. 1) and time averaging, we get the system of equations

$$ \bar{\mathbf{u}} \cdot \nabla \bar{\mathbf{u}} + \nabla \bar{p} - \frac{1}{Re} \nabla^2 \bar{\mathbf{u}} = \mathbf{f} \qquad (3a) $$

$$ \nabla \cdot \bar{\mathbf{u}} = 0, \qquad (3b) $$

where $\bar{p}$ is the average pressure field. In this framework, the forcing term $\mathbf{f}$, acting as a closure term, is identified as the Reynolds stress tensor. In principle, it can be directly computed when data are available as

$$ \mathbf{f} = -\nabla \cdot (\overline{\mathbf{u}'\mathbf{u}}). \qquad (4) $$

In practice, computing the forcing term in this way requires the solution of DNS or time-resolved experimental measurements, as the fluctuations do not directly depend on the meanflow. This issue is commonly referred to as the closure problem and, for turbulent flows, a number of approximated models based on the Boussinesq hypothesis such as $k-\epsilon$ or $k-\omega$ [56] or more complex models such as explicit algebraic Reynolds stress model [54] or differential Reynolds stress models [12] have been introduced. Here, we will rather focus on a data-assimilation approach in the footprints of [21] and [37], where assumptions such as the Boussinesq one are not used. In this context, we aim at approximating the closure term as a forcing term of an optimization process without relying on an adjoint-based loop. A model based on Graph Neural Networks (GNN) is trained as a surrogate model where the output is the forcing term (Eq. 4) and the input is given by the meanflow $\bar{\mathbf{u}}$.

## 3 Design of experiment and numerical setup

The main goal of the present work is to develop a GNN-based surrogate model emulating the closure term $\mathbf{f}$ given as the target output and the meanflow as input. Thus, in order to train the GNN, a fluid dynamics dataset is built. The training dataset is generated using high-fidelity Direct Numerical Simulations (DNS) performed with an in-house Finite Element Method (FEM) solver implemented in the FEniCS python [2] framework. This solver directly resolves the full Navier-Stokes equations without relying on turbulence models, ensuring accurate representation of the flow dynamics. The computational setup for the simulations, including mesh details and solver validation, is described in Sec. 3.3 and App. A. The DNS results, including drag and lift
<figure><img src="image_2.png" /><figcaption>**Fig. 1:** (a) Stream-wise component of the meanflow $\bar{u}$ and vorticity isolines $\omega = \nabla \times \mathbf{u}$ for the flow past a cylinder at $Re = 150$. (b) For the same case, the stream-wise component of the closure term $\mathbf{f}$ is shown. In both cases, only a portion of the domain is shown.</figcaption></figure>

coefficients for the benchmark cylinder case, were validated against established literature results to ensure accuracy and consistency. In this section, we consider bluff bodies geometries at different Reynolds numbers, such that a comprehensive dataset is obtained to handle a larger variety of flow scenarios. In the following, the reference case – i.e. the cylinder flow – is introduced in Sec. 3.1, while the flows past random geometries are discussed in Sec. 3.2. We briefly report on the numerical simulations performed for this scope in Sec. 3.3 and more deeply in App. A.

## 3.1 Cylinder flow

The unsteady wake developing past a cylinder is a well documented benchmark in fluid dynamics literature, and it is often found as reference case in works of data-assimilation for the assessment of novel techniques; it exhibits steady behaviour until a critical Reynolds number of $Re_c \approx 46.7$, when a supercritical Hopf bifurcation occurs [23, 40]. At higher Reynolds numbers, the baseflow becomes an unstable solution and the unsteady flow develops into a limit cycle known as von Karman street (Fig. 1(a)). This behaviour can be observed until $Re = 150$, for two-dimensional (2D) cases. Beyond this point, alongside periodic vortex formation, irregular velocity fluctuations begin to emerge [3]; note that transition to turbulence occurs when three dimensional cases are considered, which is beyond the scope of the present work. Instead, this work aims at a quantitative analysis of GNN model emulators trained in low-data limits leveraging active learning processes; as such, we will benchmark the proposed algorithm on simpler 2D scenarios exhibiting the limit cycle behaviour, in the range $50 \le Re \le 150$; we perform numerical simulations at different $Re$ numbers in the mentioned range, with $\Delta Re = 10$, in order to collect the necessary data to train the GNN model and compose the dataset. The numerical simulations are resolved in time, while the meanflow $\bar{u}$ and the forcing term $\mathbf{f}$ (Eq. 4) are computed by averaging on-the-fly until convergence. As an example of this data, Fig. 1(a) illustrates the stream-wise component of the meanflow $\bar{u}$ alongside the vorticity isolines $\omega = \nabla \times \mathbf{u}$, while Fig. 1(b) shows the stream-wise component of the closure term $\mathbf{f}$ for $Re = 150$.
<figure><img src="image_3.png" /><figcaption>**Fig. 2:** Flow past a bluff body with a randomly generated shape at $Re = 90$. (a) Geometry of the bluff body with $\alpha = 60^{\circ}$ and unitary diameter of the circumscribed circle around the body. The resulting meanflow and the vorticity isolines $\omega = \nabla \times \mathbf{u}$ are shown in (b), while (c) shows the forcing stress term in a portion of the domain.</figcaption></figure>

## 3.2 Flow around random shapes

As mentioned previously, the dataset can be enriched by including different geometries; thus, alongside with the flow past a cylinder, we include in the dataset flow fields around obstacles of random shapes. In principle, we should determine the critical Reynolds number ($Re_c$) at which the flow around each of the random geometries under analysis develops into unstable baseflows. This would require an extensive campaign of simulations and a comprehensive stability analysis, along the footprints of the recent works by [15, 16]. Therefore, we adopt a more pragmatic approach by retaining the geometries that in the range $50 \le Re \le 150$ develop unsteady flows.

Regarding the definition of the objects geometries, each shape is defined by a set of splines connected at the location of 4 control points, assuring a $C1$ continuity between them. The control points (P1 to P4 in Fig. 2a) are located along two orthogonal axes rotated of an angle $\alpha$ with respect to the Cartesian coordinates system of reference. Their distance from the origin $O(0,0)$ is defined in the range $|s| \in [0.1, 0.6]$, with a discrete step size of $\Delta s = 0.1$. An additional degree of freedom is added for controlling the angle $\alpha \in [0^{\circ}, 90^{\circ}]$, with a discrete step size of $\Delta\alpha = 30^{\circ}$.

Finally, the random shapes generation script prevents the repetition of the same set of degrees of freedom defining the shape such that each geometry appears only once in the final dataset. Note that the centre of the circumscribed circle of each random shape does not align necessarily with the centre of the coordinate system. This is not an accidental feature: in the context of NN training, this strategy corresponds to introducing some positional noise in the training dataset, thus promoting a more robust and prone-to-generalize learning of the underlying data. Such variability in the training data further contributes to reduce the risk of overfitting and enhances the NN ability to generalize its predictions to new, unseen geometries.
<figure><img src="image_4.png" /><figcaption>Fig. 3: Sketch of the computational domain geometry. The diameter of the circumscribed circle of the bluff body, the height and length of the domain are given in non-dimensional units.</figcaption></figure>

### 3.3 Numerical setup

In our numerical setup, we use as a characteristic dimension the diameter $D$ of the circumscribed circle of each of the generated bluff bodies. Based on this dimension, the domain extension is $L_x = 27$ in the stream-wise direction and $L_y = 10$ in the transverse direction; the centre of the system coordinates is placed at $O(0, 0)$, at a distance $\Delta x = 9$ from the inlet in the stream-wise direction and $\Delta y = 5$ from the symmetry boundaries.

The flow evolves from the inlet, where a dimensionless uniform velocity $\mathbf{u} = (1, 0)^T$ is imposed, normalized with respect to the reference velocity $U_\infty$. Following the study of [21], we introduce the following set of boundary conditions

$$ \left\{ \begin{array}{l} u = 1, v = 0 \text{ at the inlet,} \\ u = 0, v = 0 \text{ on the cylinder surface,} \\ \partial_y u = 0, v = 0 \text{ on symmetry boundaries,} \\ \frac{1}{Re} \partial_x u - p = 0, \partial_x v = 0 \text{ at the outlet.} \end{array} \right. \qquad (5) $$

All the numerical simulations are initiated with null flow fields at $t = 0$. The required statistics, i.e. meanflow and second order moments, are computed on-the-fly during the simulations. The final time $T$ of the simulation is determined by a convergence criterion based on the L2-norm of the difference between two subsequent meanflows below a threshold of $10^{-8}$.

Finite element method (FEM) is used for the spatial discretization, based on the FEniCS python library [2]; time marching is performed by second order backward differentiation formula (BFD). The meshes used for the numerical simulations are refined near the obstacle and in the wake region to best capture the flow dynamics. Depending on the shape of the obstacle, our meshes count approximately 13500 nodes for the whole domain; more details on the adopted numerical schemes and convergence of simulations are provided in A.

While this work focuses on 2D cases, the methodology could be extended to 3D flows.
Such an extension would involve addressing increased computational challenges due to the higher dimensionality of the graph representation, which significantly increases the number of nodes and edges. Approaches such as multi-GPU parallelization, subgraph-based training frameworks [5] or scalable GNN frameworks [34] could mitigate these challenges, paving the way for future applications of the proposed method to three-dimensional configurations.

# 4 Methodology

In this section, we aim at providing a concise overview of some of the fundamental features that characterize a GNN. An exhaustive and general description of the main features of this architecture can be found in [27]. Here, we focus on the main process on which a GNN is based: the diffusion of the information extracted from nodes and connections across the whole network using the message passing (MP) algorithm; Sec. 4.1 is mainly dedicated to MP, while Sec. 4.2 reports on the data structure. The training process is described in Sec. 4.3.

## 4.1 Neural Network architecture

Message passing can be divided into three fundamental steps, which we describe in the following as applied to our specific case. A pictorial sketch summarising this process is shown in Fig. 4.

1. **Message creation** – on each node $i$, an embedded state associated with the vector $\mathbf{h}_i$ is created. It is initialized with a zero state and it will embed information as long as the message passing proceeds. The dimension $d_h$ of the vector $\mathbf{h}_i$ is constant across all nodes and it is a fundamental model hyperparameter that defines, along with the number of update layers $k$, the expressivity [26] of the GNN, i.e. the ability of the model to represent functions of a certain complexity. The embedded state does not have a direct physical meaning.

2. **Message propagation** – in this step, the information is propagated. To capture the convective/diffusive dynamics of the underneath system, i.e. reproducing a model for the RANS equations, we transmit the message in both directions: from each node $i$ to its neighbour nodes $j$ and vice versa, from the nodes $j$ to the starting node $i$ (see Fig. 4 for a sketch). Mathematically, for each couple of connected nodes $i$ and $j$, we can write the following relation

$$ \phi_{i,j}^{(k)} = \zeta^{(k)} (\mathbf{h}_{i}^{(k-1)}, \mathbf{a}_{ij}, \mathbf{h}_{j}^{(k-1)}). \qquad (6) $$

In what follows, the subscript will indicate the node index, while the superscript the update index. The vector $\mathbf{h}_i^{(k-1)}$ is the embedded state from the previous layer $k-1$, $\mathbf{a}_{ij}$ denotes the directed connections between the node $i$ and the nodes $j$, while $\zeta^{(k)}$ is a generic differentiable operator such as Multi-Layer Perceptron (MLP), see for instance [24].
3. Message aggregation – finally, on each node $i$, these collected information are aggregated, providing an updated embedded state $\mathbf{h}_i^{(k)}$ for each of them

$$ \mathbf{h}_i^{(k)} = \mathbf{h}_i^{(k-1)} + \alpha \Psi^{(k)} \left( \mathbf{h}_i^{(k-1)}, \mathbf{G}_i, \phi_{i, \rightarrow}^{(k)}, \phi_{i, \leftarrow}^{(k)}, \phi_{i, \leftarrow \rightarrow}^{(k)} \right), \quad (7) $$

where the inputs are the meanflow and the Reynolds numbers, contained in the vector $\mathbf{G}_i = \{\bar{\mathbf{u}}, Re\}$ and are fed at each update; $\phi_{i, \rightarrow}^{(k)} = \frac{1}{N_i} \sum_{j \in N_i} \phi_{i, j}^{(k)}$ is the message that the node $i$ sends to its $N_i$ neighbours, averaged over the number of neighbours; $\phi_{i, \leftarrow}^{(k)} = \frac{1}{N_i} \sum_{j \in N_i} \phi_{i, j}^{(k)}$ is the message that the node $i$ receives from its $N_i$ neighbours, averaged over the number of neighbours; $\phi_{i, \leftarrow \rightarrow}^{(k)} = \phi_{i, i}^{(k)}$ the message that the node $i$ sends to itself to avoid loss of information as the process advances. $\Psi^{(k)}$ is a generic differentiable operator that can be approximated by an MLP also in this case; $\alpha$ is a relaxation coefficient, included in the model hyperparameters, that allows to scale each update of the embedded states with the previous one during the message passing loop. This step marks the exact point where external inputs are introduced into the training scheme. Unlike the initialization phase, where the embedded states $h_i$ were zero initialized without incorporating flow-specific information, this step explicitly combines $\mathbf{G}_i$ with the other quantities defined in Eq. 7, stacking them together at each message passing iteration.

Note that the use of the mean operator in $\phi_{i, \rightarrow}^{(k)}$ and $\phi_{i, \leftarrow}^{(k)}$ allows us to collect, for each of the $i$-th embedded states, the information propagating from a number of neighbours that can differ for each of the $i$-th nodes. Nonetheless, different permutation-invariant functions can be introduced in this step (maximum, sum, concatenation...).

At the end of the message passing process, the information on each node has been handled as described and thus carries with it information from any other point of the graph including their relative distances. In order to allow the message on each node to reach any other node of the mesh, our intuition is that – in principle – the number of updates $k$ should cover the longest geodesic path that can be defined on the mesh [17]; in practice, the number of updates in the present work is fixed and optimized using genetic algorithms (see [1] and B for details). The latest embedded state is projected back to a physical state using a decoder, once again approximated by a differentiable function that in our case is an MLP.

## 4.2 Data structuring

The application of GNN combined with unstructured meshes revolves around the idea of representing mesh data as a graph. In this representation, each node of the mesh is coupled one-to-one to a node in the GNN. In order to achieve this overlap, we structure the data in a tensor form that can be handled as a graph while preserving the adjacency property from the mesh. We build, for each case:

* **Node Feature Matrix A:** This matrix $A \in \mathbb{R}^{n_i \times d_h}$ represents the features of each node in the mesh. Here, $n_i$ is the number of nodes of the mesh, and $d_h$ is
<figure><img src="image_5.png" /><figcaption>**Fig. 4:** Pictorial sketch of the Graph Neural Network message-passing process: given the *i*-th node, **h**<sub>*i*</sub> is the embedded state defined on *i*; *a*<sub>*ij*</sub> are the directed connections between the nodes; **G**<sub>*i*</sub> represents the external inputs injected in each node (see Sec. 4.1 for details).</figcaption></figure>

the dimension of the embedded state in which we project the FEM quantities; each
node is assigned with a unique index.

* **Edge Index Matrix C:** The matrix $C \in \mathbb{R}^{c \times 2}$ defines the edges in the mesh. Here $c$ represents the number of edges of the mesh, and each row contains a pair of indices denoting connected nodes.

* **Edge Attribute Matrix D:** The matrix $D \in \mathbb{R}^{c \times 2}$, contains the computed distances in $x$ and $y$ directions between the connected nodes. This distances include their signs, and $D$ shares the same row indices $c$ as $C$.

The matrices are defined as follows:

$$
\mathbf{A}_{n_i, d_h} = \begin{bmatrix} a_{1,1} & a_{1,2} & \cdots & a_{1, d_h} \\ a_{2,1} & a_{2,2} & \cdots & a_{2, d_h} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n_i, 1} & a_{n_i, 2} & \cdots & a_{n_i, d_h} \end{bmatrix},
$$

$$
\mathbf{C}_{c,2} = \begin{bmatrix} i_1 & j_1 \\ i_2 & j_2 \\ \vdots & \vdots \\ i_c & j_c \end{bmatrix}, \quad \mathbf{D}_{c,2} = \begin{bmatrix} x_{i_1} - x_{j_1} & y_{i_1} - y_{j_1} \\ x_{i_2} - x_{j_2} & y_{i_2} - y_{j_2} \\ \vdots & \vdots \\ x_{i_c} - x_{j_c} & y_{i_c} - y_{j_c} \end{bmatrix}.
$$

Each column of the matrix $\mathbf{A} \in \mathbb{R}^{n_i \times d_h}$ is assigned as a feature vector to a single neuron on the MLPs structure employed in the MP step of the GNN, namely $\zeta$ (see Sec. 4.1). These feature vectors are set to zero at the beginning of the training process
(see Sec. 4.3). The MLP structure is defined by the dimension $d_h$ of the embedded state, rather than by the number of features per node. This ensures that the MLP architecture remains unaltered regardless of the number of nodes in the mesh. Indeed, changing the number of nodes in the mesh only affects the row index $n_i$ in matrix $A \in \mathbb{R}^{n_i \times d_h}$, the number of connections $c$ in matrix $C \in \mathbb{R}^{c \times 2}$, and the distances in $D \in \mathbb{R}^{c \times 2}$ but keeps the MLP structure intact, as it is defined on the dimension $d_h$ of the embedded state. This ensures that the same MLP architecture can be used for different FEM simulations, regardless of varying geometries or node numbers. Moreover, the GNN is not informed with absolute spatial coordinates as inputs, meaning it does not learn any specific mapping between spatial positions and the target output. This ensures that the GNN maintains invariance with respect to the position of the bluff body and mesh refinement, enhancing its generalization capability across different simulation setups. The ability of interacting with unstructured meshes, learning from various geometries, and reconstructing forcing fields on different geometries is rather straightforward when considering GNN architectures as compared with standard techniques such as CNN, thanks to their versatility. Finally, the model is equipped with a FEM-GNN interface that enables the mapping of numerical results from the GNN to a vector field computed or stored in the FEM environment.

## 4.3 Training algorithm

The training framework adopted for the presented GNN is described below, with reference to Fig. 5. The starting point is $A^0$ which represents the matrix composed by all the embedded states defined on the nodes, initialized at zero state. Along with the externally injected quantities ($G$, the meanflow and the Reynolds number $Re$ in our case), it is provided to the message passing process $MP^1$ (see Sec. 4.1). This latter will output an updated version of the embedded state tensor, $A^1$ which will pass through a decoder, $D^1$, a multi-layer perceptron (MLP) trainable function in charge of reconstructing a meaningful physical state ($\hat{f}^1$) from the embedded state which, in our case, is the closure term of the RANS equations. Finally, this physical state is compared with the ground truth that comes from the DNS ($\hat{f}$) via a loss function $\ell = (\hat{f} - \hat{f})^2$. The entire process is repeated in each layer until the final layer $k$ is reached. The last prediction ($\hat{f}^k$) represents the output of the entire algorithm. Following the procedure framework in [17], all the intermediate loss values are considered, in order to robustify the learning process, through a global loss function:

$$L = \sum_{k=1}^{\bar{k}} \gamma^{\bar{k}-k} \left( \frac{1}{n_i} \sum_{i=1}^{n_i} \ell_i \right), \quad (8)$$

in which $n_i$ is the number of nodes, $\bar{k}$ is the number of update layers and $\gamma$ is an hyper-parameter used to control the importance of the contributions associated with each update of the embedded states. The relative importance of the message is weighted by the exponential term $\gamma^{\bar{k}-k}$, with $\gamma < 1$; thus, the latest updates $k$ will be the richest in information and have the highest importance in the process. In B, the list of hyperparameters and their values after optimization are reported.
<figure><img src="image_6.png" /><figcaption>**Fig. 5:** Figure adapted from [17] showing the overall framework of the GNN training procedure. $MP^k$ are the message passing algorithms; $D^k$ are the $k$ decoders trainable MLPs; $A^k$ are the $k$ matrices containing the embedded states from each node; $G$ is the vector containing the input injected in the GNN. Refer to Sec. 4.3 for details.</figcaption></figure>

# 5 Results

The present section is dedicated to the presentation of the results. We introduce a proof-of-concept of the approach in Sec. 5.2; in this first analysis, we train the GNN with a dataset composed of data obtained from simulations of the flow past a cylinder at different Reynolds numbers. The aim is to verify the robustness of the approach with respect to unseen conditions, using the benchmark cases in Sec. 5.1. In Sec. 5.2, we focus on the appropriate choice of the dataset and how the selection of training data can impact the generalization capabilities of the GNN. In Sec. 5.3, we include in the dataset flow fields from simulations of wakes past randomly generated bluff bodies (details of shape generation in Sec. 3.2); the effects of dataset augmentation and dataset expansion are discussed in Sec. 5.3.1 and Sec. 5.3.2, respectively, while data selection by active learning criteria is discussed in Sec. 5.3.3. A quantitative comparative analysis of these training approaches is summarized in Sec. 5.4; in the following, we introduce the 4 reference cases used for all the comparisons discussed in this section.

## 5.1 Test Cases

We introduce 4 benchmarks to assess the performance of the GNN model. The test cases, labelled from Case 1 to Case 4, are designed to evaluate specific aspects of the GNN's prediction. We cover a broad spectrum of scenarios, including variability in Reynolds number and geometry, changes in the position with respect to the inlet and number of the obstacles.

1. **Case 1:** in this scenario, we consider the flow past a cylinder, where the Reynolds number is increased to $Re = 200$. The considered Reynolds number introduces a flow regime characterized by stronger wake dynamics and increased complexity, as it exceeds the training dataset interval, ranging in $50 \le Re \le 150$. Thus, it is used to test the model's capability to extrapolate beyond the training data.
2. **Case 2:** here, we introduce as a test case the flow past a bluff body of random shape not included in the training dataset at $Re = 120$, in order to assess the model generalization capability to unseen shapes.

3. **Case 3:** this case involves data from simulation of a flow past a bluff body of random shape, not present in the training dataset, at $Re = 100$. The obstacle is positioned downstream of the reference position used in the dataset. The main focus in Case 2 and Case 3 remains on generalization to new shapes and configurations without introducing unnecessary variability in flow characteristics.

4. **Case 4:** this test considers a two side-by-side cylinders configuration; two bluff body obstacles are present in the flow, in contrast to the single obstacle used for the training. The Reynolds number is $Re = 90$ and the aim is to test geometries inducing different dynamics, here the one stemming from multiple interacting bluff bodies. At the same time, this lower Reynolds number ensures the case remains well within the training dataset range.

Fig. 6 shows the streamwise component of the meanflow velocity $\bar{u}$, the isolines of the vorticity $\omega = \nabla \times \mathbf{u}$ (left column) and streamwise component of the forcing stress $\mathbf{f}$ (right column) computed using numerical simulations for the 4 benchmarks Case 1-4.

Two metrics will be used to evaluate the performance of the GNN. The first one will be a comparison between the ground truth $\mathbf{f}$ and the GNN prediction for each of the cases, considered as the relative error $\epsilon$, based on the L2-norm, defined as

$$ \epsilon = \frac{\|\mathbf{f} - \hat{\mathbf{f}}\|}{\|\mathbf{f}\|} = \frac{\left[\int_{\Omega}(\mathbf{f} - \hat{\mathbf{f}})^2 d\Omega\right]^{1/2}}{\left(\int_{\Omega}\mathbf{f}^2 d\Omega\right)^{1/2}}. \quad (9) $$

The second metric employed arises from the necessity to assess the accuracy on the meanflow reconstruction $\hat{\mathbf{u}}$ from the GNN prediction with respect to the ground truth meanflow $\mathbf{u}$. It is defined as

$$ \delta = \frac{\|\mathbf{u} - \hat{\mathbf{u}}\|}{\|\mathbf{u}\|} = \frac{\left[\int_{\Omega}(\mathbf{u} - \hat{\mathbf{u}})^2 d\Omega\right]^{1/2}}{\left(\int_{\Omega}\mathbf{u}^2 d\Omega\right)^{1/2}}. \quad (10) $$

In the previous equations $\Omega$ is the computational domain (Sec. 3.3), $\hat{\mathbf{f}}$ is the GNN prediction and $\mathbf{f}$ is the closure term of the RANS equations coming from the DNS. Plugging $\hat{\mathbf{f}}$ in Eq. 3 and solving the inverse problem we can obtain $\hat{\mathbf{u}}$, a reconstruction of the meanflow based on the GNN prediction.

The proposed framework is rather versatile, as its architecture relies on capturing local and global flow relationships through graph representations. With appropriate training data, it can be applied to a variety of flow scenarios beyond those explored in this study. While the focus here is on bluff body configurations, the methodology's reliance on generalizable features makes it well-suited for adaptation to other fluid dynamics problems, provided the underlying data is representative of the target flow conditions.
<figure><img src="image_7.png" /><figcaption>**Fig. 6:** Test cases used as a benchmark of the learning strategies discussed in Sec. 5. From the top raw to the bottom one: Case 1, cylinder flow at $Re = 200$; Case 2, flow past a random-shaped bluff body at $Re = 120$; Case 3, flow past a random-shaped bluff body shifted in the computational domain at $Re = 100$; Case 4, flow past a two side-by-side cylinders configuration at $Re = 90$. For each of the cases, the left column shows the stream-wise component of the meanflow $\bar{u}$, along with the vorticity isolines $\omega = \nabla \times \mathbf{u}$, while the right column show the stream-wise component of the forcing stress $\mathbf{f}$.</figcaption></figure>

## 5.2 Proof-of-concept training: flow past a cylinder flow

Here, we consider the baseline results that in the following we will indicate as proof-of-concept (PC). The training dataset contains data obtained from DNS simulations of the flow past a cylinder; it is composed by 11 pairs of meanflows $\bar{u}$ (input) and Reynolds forcing stress $\mathbf{f}$ (target), in the interval $50 \le Re \le 150$, with a stride of $\Delta Re = 10$. As previously mentioned, all the cases in the dataset exhibit a von Karman street instability in the wake of the bluff body.

Fig. 7 shows the training and validation loss curves for this training approach. The training loss behaviour underlines the model's ability to capture the dynamics of the fluid flows and to learn the patterns in the training dataset. However, a significant
<figure><img src="image_8.png" /><figcaption>**Fig. 7:** Curves for the training and validation loss for the proof-of-concept baseline, trained until 3000 epochs. The training dataset is composed by 11 meanflow-forcing pairs stacked at different Reynolds numbers in the range $50 \le Re \le 150$ with $\Delta Re = 10$. The validation set is formed by the test cases presented in Sec. 5.1.</figcaption></figure>

gap can be observed between the training and validation loss curves. This discrepancy suggests a potential problem of overfitting to the training data. Further evidence supporting this claim comes from experiments with dropout, which reduced the gap between training and validation losses, indicating improved generalization. Similarly, increasing the model's expressivity by using deeper architectures or higher-dimensional embeddings amplified this divergence, with training loss decreasing further while validation loss stabilized or slightly increased—behavior characteristic of overfitting. This problem is not unexpected, considering that the training dataset is composed exclusively of data obtained from simulations past the same geometry. Thus, when the GNN model is tested with unseen shapes or fluid dynamic conditions, the performance drops.

Regarding the test cases, in Case 1 (Fig. 8a), the GNN prediction shows good accuracy in reproducing the fluid structures at higher Reynolds number $Re = 200$, with $\epsilon = 0.1153$. Fig. 8(b) shows Case 2, where the difference between the GNN prediction and the DNS data is primarily visible in the near wake region. This region is crucial for the development of the unstable dynamics as it corresponds to the so-called wavemaker region [23]; for cases as such of non-standard geometrical shape, the flow vortex shedding behaviour can be quite different from the one observed in standard cylinder cases. The $\epsilon$ norm in this case is $\epsilon = 0.2995$. Case 3 in Fig. 8(c) tests the GNN generalization capabilities when the bluff body changes; the main discrepancies are again observed in the near wake of the bluff body. The $\epsilon$ norm in this case is $\epsilon = 0.2084$. Finally, Case 4 in Fig. 8(d) is particularly challenging as it tests generalization capabilities in presence of multiple bluff bodies. The GNN model's predictions capture the major features of the forcing stress of each of the cylinder. However, discrepancies are mainly observed in the area between the two cylinders and extending in the far wake region, where interactions between the Reynolds stress fields of the two cylinders occur. The $\epsilon$ norm in this case is $\epsilon = 0.8704$, significantly higher compared to the previous cases.
<figure><img src="image_9.png" /><figcaption>**Fig. 8:** Stream-wise component comparison of the Reynolds stress tensor. The GNN's model is trained with a dataset composed by 11 cases of cylinder shape bluff bodies, ranging in $50 \le Re \le 150$, $\Delta Re = 10$. (a) Case 1, $Re = 200$, cylinder bluff body shape; (b) Case 2, $Re = 120$, random shape bluff body; (c) Case 3, $Re = 100$, random shape shifted bluff body; (d) Case 4, $Re = 90$, flow past two side-by-side cylinders.</figcaption></figure>

In all the cases discussed so far, the presence of some numerical errors in all the predictions is noted, but these are primarily attributed to statistical noise or inherent errors typical of neural network models. Neural networks indeed, by their nature, include elements of statistical uncertainty due to factors such as the stochastic nature of their training algorithms (e.g., random initialization of weights, batch selection during training), and the approximate nature of the model that represents the underlying physics. Beside this aspect, this preliminary analysis based on the benchmark cases suggests that the GNN model performance could be further improved with a training set including a larger variety of bluff body shapes. Therefore, in what follows, we study the effects of the data augmentation by including in the training dataset random shaped bluff bodies with the aim of enhancing the model's robustness and generalization capabilities.

## 5.3 Data augmentation and active learning: fluid flows past random geometries

In the following, we study the effects of data amount and quality on the training process and generalization properties of the GNN models. This study is structured into distinct steps.

1. **Data augmentation (DA)** (Sec. 5.3.1); this first approach involves dataset augmentation, meaning that we incorporate bluff bodies of random shapes cases into the existing dataset, while maintaining the same dataset size used in the PC case,
<figure><img src="image_10.png" /><figcaption>**Fig. 9:** Curves for the training and validation loss for when DA is performed by introducing 11 cases of flows past random geometries in the dataset. The Reynolds number varies in the range $50 \le Re \le 150$ with $\Delta Re = 10$. The training runs until 3000 epochs. The validation set is formed by the test cases presented in Sec. 5.1.</figcaption></figure>

i.e. 11 cases. This method is aimed at diversifying the range of geometries used during the training of the model, without increasing the data volume.

2. **Dataset expansion** (Sec. 5.3.2); in this second approach, we expand the dataset by including 33 cases of bluff bodies of random shapes. We assess the quality of the used data in terms of sensitivity of the model to specific configurations by performing a k-fold validation.

3. **Active learning (AL) data selection** (Sec. 5.3.3); in this last approach, the training set is built by adding progressively data chosen by similarity criterion. The goal is to develop a surrogate model that can generalize at its best to unseen cases based on a given dataset.

### 5.3.1 Data augmentation with random shapes

In this phase of the study, we employ a stratified random sampling criterion to select shapes among the randomly generated configurations detailed in Sec. 3.2. This latter involves the generation of the entire set of shapes for each Reynolds number in the interval $50 \le Re \le 150$, $\Delta Re = 10$ and randomly select 1 shape for each value of the Reynolds number, for a total of 11 shapes that compose the new training dataset. The chosen geometries are unique. On the selected shapes we perform DNS simulations in order to obtain the meanflow (input) and the forcing stress (target) used to train the GNN.

Fig. 9 shows the training and validation loss curves for this second training approach. The training loss demonstrates a consistent downward trend, highlighting the learning effectiveness of the GNN model. Notably, while there remains a significant gap between the training and validation loss curves, this gap is less pronounced when compared to the initial training approach (Fig. 7). This reduced gap is indicative of diminished overfitting and suggests that the introduction of a wider variety of shapes and flow conditions into the training dataset leads to improved GNN model's generalization capabilities.
<figure><img src="image_11.png" /><figcaption>**Fig. 10:** Stream-wise component comparison of the Reynolds stress tensor. The GNN's model is trained with a dataset composed by 11 cases of random shaped bluff bodies, ranging in the interval $50 \le Re \le 150$, $\Delta Re = 10$. (a) Case 1, $Re = 200$, cylinder bluff body shape; (b) Case 2, $Re = 120$, random shape bluff body; (c) Case 3, $Re = 100$, random shape shifted bluff body; (d) Case 4, $Re = 90$, flow past two side-by-side cylinders.</figcaption></figure>

Regarding the GNN predictions on test cases, Case 1 (Fig. 10a) shows larger discrepancies in the forcing prediction in the neighbourhood of the cylinder compared to the training approach discussed in Sec. 5.2, $\epsilon = 0.5033$. An interpretation for this result is given by the nature of the PC training dataset that solely relies on cylindrical geometries; thus, the resulting GNN model is specialized in the prediction of the cylinder flows. In contrast, the current training strategy involves flows past random shapes: from one hand, this choice broadens the GNN's flexibility to capture diverse flow conditions by enhancing generalization capabilities; on the other hand, it simultaneously decreases the prediction accuracy for cylinder shaped cases.

Case 2 (Fig. 10b) and Case 3 (Fig. 10c) show mainly numerical and statistical noise in the difference fields while the main features of the Reynolds stress tensor are well reproduced in the prediction. Contrary to what is observed in Case 1, this training approach leads to an efficient GNN model in predicting the features of the flows past random shaped bluff bodies. Error norm is comparable for the two test cases, namely $\epsilon = 0.2063$ for Case 2 and $\epsilon = 0.1999$ for Case 3. Finally, Case 4 is characterized by inaccuracies primarily located in the wake region similar to the one already observed in the PC case (for a reference, the reader can compare Fig. 8d and Fig. 10d). However, it's important to note that these errors are quantitatively less significant to those observed in the first training approach, with $\epsilon = 0.6312$.
<figure><img src="image_12.png" /><figcaption>**Fig. 11:** Training and validation loss curves for models trained using 10 *k*-fold up to 3000 epochs. (a) shows the training curves, (b) the corresponding validation curves for each fold.</figcaption></figure>

### 5.3.2 Dataset Expansion

A common technique to improve generalization in neural network models is to enlarge the volume of data the model is trained with. This approach reduces the risk for the NN model to overfit to specific conditions observed in a small dataset and is more likely to capture the underlying phenomena. Therefore, we study the effect of tripling the amount of data in the training dataset. With the same stratified random sampling approach described in Sec. 5.3.1, we select 3 cases for each of the 11 Reynolds numbers chosen in the interval of reference, resulting in a 33 cases training dataset. In order to analyse the sensitivity of the GNN model's generalization capabilities with respect to the training dataset, we employ a *k*-fold validation test. First, we divide the training dataset into *k* groups and train *k* different models. For each of the *k* models, we use *k*−1 groups as training dataset and the *k*-th remaining group as the validation dataset. Although this approach is not feasible in practical applications, it is very informative for assessing preliminarily the impact of the quality of data on the final prediction.

Fig. 11(a) shows the training curves of the 10 models, one for each of the *k* folds, while Fig. 11(b) displays their corresponding validation curves. The training curves are rather close, independently from the chosen dataset, thus suggesting a rather consistent behaviour of the chosen GNN architecture. However, the validation loss varies significantly. This variability can be attributed to the distribution of the training dataset. Inconsistent or unbalanced data distributions may lead the model to being trained on subsets of data that do not adequately represent the overall dataset, thus impacting negatively on the validation results. To address this challenge and identify the most effective dataset, we introduce AL in the following.

### 5.3.3 Active Learning data selection

AL allows to dynamically adjust the dataset, targeting the data that contribute most significantly to the model's generalization capability. We exploit a similarity criterion
from the GNN perspective, following the study by [13]. The key idea is to avoid including in the training dataset cases that do not lead to any significative improvement of the model from the generalization viewpoint; instead, in order to induce diversity in the GNN model, we need to select cases that are as most as different from each other from the GNN perspective. We can achieve this goal by comparing for each of the available data pairs the vector gradient of the cost function with respect to the learnable parameters of the GNN $\theta$ by means of scalar products. As comparison metric, we employ the cosine similarity angular distance

$$ \cos(\beta) = \frac{\mathbf{a} \cdot \mathbf{b}}{\|\mathbf{a}\| \|\mathbf{b}\|}, \qquad (11) $$

where $\mathbf{a} \in \mathbb{R}^m$ and $\mathbf{b} \in \mathbb{R}^m$ are two generic $m$-dimensional vectors and $\beta$ is the angle between them. The basic idea is that each gradient vector steers the neural network's state to a specific direction within the solution space. Therefore, including multiple data points whose gradients point in the same direction might be redundant. Instead, our aim is to select and include in the training dataset as many different directions as possible, in order to explore more extensively the solution space (see also C).

The selection process stops when a predefined similarity threshold is reached, such that only the pairs exhibiting a similarity below this threshold are added in the dataset. Conversely, cases that show a similarity value exceeding the threshold are discarded. However, the similarity threshold does not have a direct, interpretable meaning and needs to be tuned as an input parameter of the process, according to the required performances. In general, a lower similarity threshold means that more cases will be included in the training dataset. While this can enhance models behaviour to diverse scenarios, it may also lead to a more complex and time-consuming training process. On the other hand, a too restrictive similarity threshold might exclude potentially valuable training data. Therefore, an optimal similarity threshold is chosen as a trade-off between the training computational costs and the model generalization performances.

In our study, we explore three different similarity threshold values, $\cos(\beta) \in [0.7, 0.8, 0.9]$. This approach allows to observe the impact of varying levels of data inclusion on the model's training and performance. Results shown in Fig. 13 and Fig. 12 are obtained with a similarity threshold of 0.8, which results in only 6 selected random bluff body cases for the training dataset; the results for the other threshold values are detailed in Sec. 5.4). Fig. 12 reports the training and validation loss curves: the reduced disparity observed between the two curves – compared to the results from the previous learning approaches – underscores a notable decrease in the issue of overfitting. The training is initiated with a single case in the training dataset, specifically the cylinder case at $Re = 120$.

Once the convergence of the vector gradients is reached (see C), a new case in training set is selected; since this process is a seek-and-include algorithm that enlarges the dataset at each step, two different approaches can be applied. The first strategy is inspired by the curriculum learning framework [7], where each selected case is added to the ongoing training. Thus, the model is progressively updated. A second approach, applied in this work, consists in reinitializing the GNN weights every time the dataset
<figure><img src="image_13.png" /><figcaption>**Fig. 12:** Training and validation loss curves are shown when AL with similarity threshold 0.8 is performed until 3000 epochs. The training dataset is formed by 6 random shaped bluff bodies. Their Reynolds number varies in the interval $50 \le Re \le 150$. The validation set is formed by the test cases presented in Sec. 5.1.</figcaption></figure>

is enlarged with a new pair. The rationale behind this choice is based on the dynamics of the solution space, that changes when new data are added in the training set; in this sense, there is no guarantee that the new minima will be "closer" to the previous ones than to the initialization point of the GNN. Nonetheless, we observed that for the analysed flow cases the two training strategies lead to negligible differences in the final results.

Considering the benchmarks for the model assessment, Case 1 leads to $\epsilon = 0.5507$ between the GNN predictions of the Reynolds stress tensor and the DNS ones. The error is $\epsilon = 0.1518$ for the random shape bluff body case (Case 2), $\epsilon = 0.3595$ for the downstream shifted random shape bluff body (Case 3) and $\epsilon = 0.5243$ for the two side-to-side cylinder configuration (Case 4). The results respect the symmetry of the solutions with respect of the $y$ axis and only minor inconsistencies can be observed in the far wake; we note that the errors are concentrated in the region immediately downstream of the bluff body, but we interpret these errors as numerical or stochastic in nature rather than being prediction inaccuracies of the fluid structures.

## 5.4 Quantitative comparison

In conclusion, we summarize the results in Tab. 1. A direct comparison between the first two training approaches (PC and DA) reveals that enriching the training dataset with random geometries enhances the generalization capabilities of the GNN. This behaviour is not unexpected, although an exception is observed when extrapolating at higher Reynolds numbers (Case 1), where we observe that the PC approach outperforms DA. As already discussed, this behaviour can be understood by observing that the GNN model trained following the approach PC is specialized in predicting the flow past a cylinder as it is trained on this specific geometry.

On the other hand, the AL training approach allows to obtain GNN models demonstrating overall superior performance in terms of generalization capabilities as compared to PC and DA approaches, in particular for Case 2 and Case 4, while in
<figure><img src="image_14.png" /><figcaption>**Fig. 13:** Stream-wise component comparison of the Reynolds stress tensor. The GNN's model is trained with a dataset composed by 6 cases of bluff bodies of random shape selected with the AL approach, ranging in the interval $50 \le Re \le 150$, $\Delta Re = 10$. (a) Case 1, flow past a cylinder at $Re = 200$; (b) Case 2, flow past a random shaped bluff body at $Re = 120$; (c) Case 3, random shaped bluff body at $Re = 100$, shifted downstream; (d) Case 4, flow past two side-by-side cylinders at $Re = 90$.</figcaption></figure>

Case 3 performance are essentially comparable. The mean error decreases in all the 4 test cases in the AL approach, indicating an improvement in the global performance of the model. An important aspect regarding the active selection of data for the training is the amount of data. Good performances are achieved with only 6 pairs when a threshold 0.8, although it is also possible to observe a strong variation in the total number of pairs used for the training as a function of the chosen threshold. It is crucial to observe, however, that for all the AL cases the introduction of a selection criterion guarantees that the chosen data points lead to good performance in terms of generalization.

When comparing our work with existing literature, the notable aspect is that our approach achieves comparable accuracy with significantly fewer training cases. For instance, in [14] the authors utilize 2000 cases in their training dataset for steady-state incompressible flow around a cylinder at $Re = 10$. In contrast, we use only 19 cases for the largest dataset used and still manage to generalize to different Reynolds numbers and bluff body positions, with comparable accuracy results. In [29], 500k cases are used for the training dataset although an unsteady flow is predicted using CNN. It is worth noting that the use of a GNN architecture enables to generalize on different geometries and $Re$ numbers, an aspect that is not addressed in [29]. Finally in [50], a GNN is employed for predicting time-averaged steady flow. The dataset consists of 12800 data points while the GNN model has a complexity of over 30M parameters, prohibitive for most practical applications. Our GNN architecture, on the contrary, can
<table><thead><tr><th rowspan="2">Cases</th><th colspan="2">PC</th><th colspan="2">DA</th><th colspan="2">AL</th><th colspan="2"></th></tr><tr><th>ε</th><th>δ</th><th>ε</th><th>δ</th><th>ε</th><th>δ</th><th>ε</th><th>δ</th></tr></thead><tbody><tr><td>Case 1</td><td><b>0.1153</b></td><td><b>0.0118</b></td><td>0.5033</td><td>0.1987</td><td>0.8671</td><td>0.1859</td><td>0.5507</td><td>0.1165</td><td>0.4132</td><td>0.1162</td></tr><tr><td>Case 2</td><td>0.2995</td><td>0.0280</td><td>0.2063</td><td>0.0160</td><td>0.1256</td><td><b>0.0079</b></td><td>0.1518</td><td>0.0166</td><td><b>0.1124</b></td><td>0.0092</td></tr><tr><td>Case 3</td><td>0.2084</td><td>0.0131</td><td><b>0.1999</b></td><td>0.0090</td><td>0.3058</td><td>0.0107</td><td>0.3595</td><td>0.0183</td><td>0.2123</td><td><b>0.0082</b></td></tr><tr><td>Case 4</td><td>0.8707</td><td>0.2683</td><td>0.6312</td><td>0.1701</td><td>0.6751</td><td>0.1930</td><td><b>0.5243</b></td><td><b>0.1565</b></td><td>0.5719</td><td>0.1653</td></tr><tr><td>Training data</td><td colspan="2">11</td><td colspan="2">11</td><td colspan="2">3</td><td colspan="2">6</td><td colspan="2">19</td></tr></tbody></table>

**Table 1:** Comparison on the 4 cases defined as benchmark for the three different training approaches: the proof-of-concept (PC) training (Sec. 5.2), the data augmentation (DA) training (Sec. 5.3.1), and the active learning (AL) strategy (Sec. 5.3.3). For the latter, we consider three similarity threshold values $cos(β) \in [0.7, 0.8, 0.9]$. The chosen metric $ε$ and $δ$ are defined respectively in Eq. 9 and Eq. 10. In the last row the number of pairs used during the training process is reported.

count up to approximately 900k parameters. In conclusion, we believe that the combination of GNN models and active learning makes our method more accessible and practical for broader applications thanks to the parsimony of the data requirements.

## 6 Conclusion

The application of machine learning techniques in fluid mechanics is often characterized by shortcomings such as overfitting, lack of robustness of the prediction with respect to unseen cases and data-hungriness. In this article, we proposed a novel approach that combines a neural architecture based on Graph Neural Networks (GNN), numerical solvers based on Finite Element Method (FEM) and an active learning procedure in order to tackle some of these limitations. Here, we consider a data-assimilation schemes that does not rely on an optimization process and use as baseline equations the Reynolds-averaged Navier-Stokes (RANS) equations. GNN models are trained as a surrogate to predict the forcing/closure term, obtained as an output of the supervised learning, while a given mean flow serves as input. The GNN architecture is particularly suitable in this study due to its adaptability to unstructured meshes and its generalization capability, as compared to other literature approaches. Moreover, this architecture allows frugal training within the low-data limit, as compare to alternative, more expensive in terms of required data, architectures.

A two-fold interface between FEM and GNN environment has been developed to transform a FEM vector field into a numerical tensor that can be handled by a NN structure and vice versa, preserving critical information throughout the process. As test-bed, we focussed on two dimensional, incompressible flows past obstacles at low Reynolds numbers, namely in the range $50 \le Re \le 150$. At these regimes, the presence of obstacles triggers instabilities developing in unsteady flows. We started by studying a cylindrical geometry in the range $50 \le Re \le 150$ as initial benchmark, in order to assess the extent to which a model based on this training dataset can be used also for unseen cases. Not surprisingly, we found good performance at unseen Reynolds number for the cylinder case. On the other hand, when the flow around the bluff bodies of
random geometry is considered, it is observed lack of accuracy in the prediction and overfitting.

In order to tackle these limitations, we explored the impact of the training data
on the generalization capabilities of the GNN in terms of quantity and quality of
data. First, we considered an extended dataset. Our results indicate that the quality
and volume of data notably affect the spectrum of unseen cases on which the model
can generalize to. Particularly, the inclusion of diverse fluid flow conditions into the
training dataset improves the overall generalization capabilities for the vast majority
of cases. Finally, we introduced an active learning data selection criterion based on
the analysis of gradient similarity, with the aim of building a dataset extending the
distribution of the data. At the best of authors knowledge, this is one of the first
applications in the community of fluid mechanics where a systematic selection of the
data is performed addressing the generalization of the NN model prediction to unseen
cases by maximizing the quality of the predictions and at same time minimizing the
amount of data used in the training set. The results clearly indicate the possibility
of improving the performance of the model, also in terms of generalization, while
keeping a rather small amount of data in the training set. The criterion is especially
relevant in the contexts where computational resources for training surrogate models
are limited, and a trade off between accuracy of the predictions and training compu-
tational cost is sought. It is stressed that the datasets used during the training process
are relatively small as compared to other approaches appeared in ML literature, as
the most expensive one consists of less than 20 pairs of snapshots. These datasets
are selected through a criterion that minimizes the number of data points required
and is robust enough to be applied to larger datasets to efficiently reduce them while
maintaining essential information.

We are currently extending our findings by incorporating physical constraint in
the learning loop through the adjoint equations associated with the assimilation loop.
This integration aims to further refine the predictive performance and generalization
capabilities of the GNN leveraging the physics constraints. The final goal of these
data-assimilation schemes is to adapt our training approach to cases where solely
limited or corrupted measurements of the flow are available, such as those based on
sparse probe measurements, noisy or incomplete data. Future work will also aim to
address the observed challenges in high-shear regions, where the inherent smoothing
effect of GNNs can lead to underrepresentation of sharp gradients. Proposed solu-
tions, including dataset augmentation and gradient-aware loss regularization, will be
investigated to enhance the model's accuracy in these critical zones. Finally, a very
promising research route would be the extension of the proposed GNN approach to
3D turbulent flow cases, as well as to more complex frameworks such as Unsteady
Reynolds-averaged Navier-Stokes (URANS) or even Large Eddy Simulation (LES),
where unsteady dynamics are explicitly resolved or modeled.

Acknowledgements

The PhD fellowship of M. Quattromini is supported by the Italian Ministry of Univer-
sity. This study has been partially funded under the National Recovery and Resilience
Plan (NRRP), Mission 4 Component 2 Investment 1.3 - Call for tender No. 1561 of 11.10.2022 Project code PE0000021, Project title Network 4 Energy Sustainable Transition - NEST, funded by the European Union - NextGenerationEU. The support from the Agence Nationale de la Recherche (ANR) through grant ANR-21-REASON is gratefully acknowledged.

## Appendix A Numerical simulations details

Time-resolved numerical simulations were performed using a python code based on the FEniCS [2] library. The inputs and the outputs of the Graph Neural Network model were obtained by averaging these data on-the-fly. More details on the numerical scheme can be found in [60]; the main script is based on the code described in [25], where an extensive discussion is reported. Here we summarize the essential aspects and provides details on the validation. From the numerical viewpoint, the spatial discretization is obtained by introducing the weak formulation based on the finite element method (FEM). In particular, the finite element used is the Taylor-Hood element, with second order elements P2 for velocity and first order elements P1 for pressure. The implemented spatial integration scheme reads as

$$ 
\left\{ 
\begin{array}{l}
\left( \frac{3\mathbf{u}^n - 4\mathbf{u}^{n-1} + \mathbf{u}^{n-2}}{2\Delta t} \right) + (\mathbf{u}^{n-1} \cdot \nabla)\mathbf{u}^n + (\mathbf{u}^n \cdot \nabla)\mathbf{u}^{n-1} \\
\quad - (\mathbf{u}^{n-1} \cdot \nabla)\mathbf{u}^{n-1} - \frac{1}{Re}\Delta\mathbf{u}^n + \nabla p^n = \mathbf{0} \\
\nabla \cdot \mathbf{u}^n = 0,
\end{array}
\right. 
\quad (A1) 
$$

where $\mathbf{u} = (u, v)^T$ represents the velocity vector, $p$ the pressure and $Re$ the Reynolds number. The time marching is performed by second order backward differentiation formula (BFD): the $n$ apex indicates the values of a quantity at the current time, with $n-1$ the values at the previous time step and $n-2$ its values two time steps before. The time step - denoted as $\Delta t$ - is chosen for granting the Courant-Friedrichs-Lewy condition, $CFL \le 0.5$ for each of the random generated shapes. The convective term is treated using the Newton and Picard methods, such that a linear system of equations is solved in each step of the temporal iteration. The basic mesh is based on the one used for the cylinder flow, our reference case; the mesh is refined in the most sensitive regions around the obstacle and in the wake past the cylinder [23]: these are the areas where high refinement is required for increasing the accuracy of the computations in terms of correct frequencies and growth rate of the instabilities. Beside the well known dynamics, the rich literature allows a detailed comparison for validation purposes. Thus, numerical solver and meshes were validated using the cylinder flow, in particular by comparing the drag coefficient $C_d$ and the lift coefficient $C_l$. Fig. A1 shows the evolution in time of $C_d$ and $C_l$, while the table Tab. A1 shows the numerical comparison of (time-averaged) $C_d$ and $C_l$ (amplitude) with respect to values in literature. The first coefficient compares remarkably well with the values in literature, while a slightly higher $C_l$ is found: this can be related to the chosen numerical box that is slightly small along the $y$ direction, thus producing some numerical blockage effects; nonetheless, for the nature of the present work, the results are satisfactory. Finally, the chosen
<figure><img src="image_15.png" /><figcaption>Fig. A1: Drag coefficient $C_d$ and lift coefficient $C_l$ of flow past a cylinder.</figcaption></figure>

<table><caption>Table A1: Comparison of $C_d$ average and $C_l$ amplitude</caption><thead><tr><td>Cylinder</td><td>C<sub>d</sub> average</td><td>C<sub>l</sub> amplitude</td></tr></thead><tbody><tr><td>Etienne and Pelletier</td><td>1.36</td><td>0.67</td></tr><tr><td>Li et al.</td><td>1.34</td><td>0.69</td></tr><tr><td>Liu et al.</td><td>1.31</td><td>0.69</td></tr><tr><td>Present result</td><td>1.32</td><td>0.73</td></tr></tbody></table>

bluff bodies are characterized by mechanisms resembling the phenomena observed in the cylinder flow, thus characterized by comparable time and space scales; the corresponding meshes are characterized by a refinement similar to the one adopted for the reference case.

## Appendix B Hyperparameters optimization

Numerous parameters define the structure of a neural network, usually denoted with the term hyperparameters. We can distinguish between model hyperparameters and process hyperparameters.

* A model hyperparameter defines the capacity of the neural network, i.e. the ability of the model to represent functions of high complexity. Thus, the capacity is directly related to the possibility of approximating a large variety of nonlinear functions.

* A process hyperparameter defines the training phase. Tuning these hyperparameters deeply modify the duration of the training, its computational costs and the way the weights are adjusted while the model evolves.
These terms are defined a-priori, before training the model, thus they need to be manually tuned as they can't be adjusted or learnt during the training process. In general, this does not pose problems in deep learning when the expressivity of the NN (e.g. the number of neurons, the number of layers, etc.) is enough to represent the complexity of the problem under investigation. On the other hand, since the GNN is trained to fulfil a specific task and because the computational cost of the GNN inference has to be affordable compared to the most sophisticated turbulence models available today, we tried to keep the GNN as parsimonious as possible.

To this end, the hyperparameters defining the architecture require optimization. Standard gradient based optimizers cannot be employed when dealing with integer numbers, like the number of neurons or layers. For this purpose gradient-free algorithms can be used. There exists dedicated libraries that can automate the tuning process through all the possible sets of hyperparameters by trying and appropriately pruning the unpromising sets of them. In this work we apply the library Optuna [1], an open-source package that combines efficiently searching and pruning algorithms. By exploring the complex solution hyperspace, a number of combinations of hyperparameters is found, among which the one outperforming the others in terms of monitored validation metrics is given by the following set

1. Embedded dimension, 35

2. Number of GNN layers, $k = 40$

3. Update relaxation weight, $\alpha = 6 \times 10^{-1}$

4. Loss function weight, $\gamma = 0.1$

5. Learning rate, $LR = 3 \times 10^{-3}$, as maximal/starting value.

## Appendix C Similarity criteria algorithm details

The similarity criteria algorithm designed to compare different data from the neural network perspective is based on the analysis of the vector gradients of a metric function (Eq. C3) with respect to the $\theta$ parameters of the neural network. In particular, the similarity comparison between two generic $m$-dimensional vectors $\mathbf{a} \in \mathbb{R}^m$ and $\mathbf{b} \in \mathbb{R}^m$ is computed using the cosine similarity, defined as

$$ \cos(\beta) = \frac{\mathbf{a} \cdot \mathbf{b}}{\|\mathbf{a}\| \|\mathbf{b}\|}, \qquad (C2) $$

where $\beta$ is the angle between the two vectors. In this context the metric we use is the Mean Absolute Error (MAE), a piecewise linear function defined as

$$ MAE = \sum_{i=1}^{n_i} |x_i - y_i|, \qquad (C3) $$

in which $x_i$ is the NN prediction on the node $i$, $y_i$ the ground truth and $n_i$ the number of nodes.

The choice of MAE is crucial for our analysis. In scenarios where the Mean Squared Error (MSE) is used, we empirically observe marked oscillations in the direction of
<figure><img src="image_16.png" /><figcaption>Fig. C2: Visualization of gradient auto-similarity convergence over multiple training epochs using MAE for 9 distinct cases within the training data set.</figcaption></figure>

the auto-similarity of the vector gradient throughout successive training epochs. Conversely, when employing MAE, the auto-similarity of the vector gradient tends to approach unity, suggesting a stable direction in the solution space, for the data under analysis. A visual representation of the consistent convergence of the vector gradient auto-similarity as function of the epochs is shown in Fig. C2, for a training process involving 9 cases in the training dataset. Training begins with a specific initial dataset; when every data point in the training dataset achieves an auto-similarity convergence exceeding a predefined threshold of 0.99, the training is stopped and we can assume that the vector gradient's direction for each instance in the training dataset has stabilized.

The following step is to assess the similarity between cases within the training dataset and those outside it. This aims at identifying the most diverse cases among those not included in the training set, which will then be added to the training dataset to enhance diversity. Firstly, the GNN runs for additional 10 epochs to obtain the vector gradients of each out-of-training dataset instance. Then, a similarity matrix is computed by cross calculating the similarity between each in-training instance and each out-of-training instance. The case that shows the lowest similarity score is also the most diverse one and enables to promote diversification in the training dataset based on available data. Note that the values are normalized using a z-score value

$$z = \frac{S - \mu}{\sigma}, \qquad (C4)$$

where S represents the similarity score for a specific data point, μ is the mean of all similarity scores, and σ their standard deviation. In summary, the approach outlined here serves as a robust method for comparing and evaluating the similarities in data
behaviour leveraging the neural network model, thereby enhancing the efficacy of the training process.

## References

[1] Akiba T, Sano S, Yanase T, et al (2019) Optuna: A next-generation hyper-parameter optimization framework. In: Proceedings of the 25rd Association for Computing Machinery International Conference on Knowledge Discovery and Data Mining

[2] Alnæs M, Blechta J, Hake J, et al (2015) The FEniCS project version 1.5. Archive of Numerical Software 3(100)

[3] Anatol R (1958) Naca report 1191

[4] Bae HJ, Koumoutsakos P (2022) Scientific multi-agent reinforcement learning for wall-models of turbulent flows. Nature Communications 13(1):1443

[5] Bai J, Ren Y, Zhang J (2021) Ripple walk training: A subgraph-based training framework for large and deep graph neural network. arXiv preprint arXiv:200207206 [cs.LG]

[6] Beck A, Kurz M (2021) A perspective on machine learning methods in turbulence modeling. GAMM-Mitteilungen 44(1):e202100002

[7] Bengio Y, Louradour J, Collobert R, et al (2009) Curriculum learning. In: Proceedings of the 26th Annual International Conference on Machine Learning, pp 41–48

[8] Brunton SL, Hemati MS, Taira K (2020) Special issue on machine learning and data-driven methods in fluid dynamics. Theoretical and Computational Fluid Dynamics 34(4):333–337

[9] Brunton SL, Noack BR, Koumoutsakos P (2020) Machine learning for fluid mechanics. Annual Review of Fluid Mechanics 52:477–508

[10] Bucci MA, Semeraro O, Allauzen A, et al (2023) Curriculum learning for data-driven modeling of dynamical systems. The European Physical Journal E 46(3):12

[11] Cai J, Angeli PE, Martinez JM, et al (2024) Revisiting tensor basis neural network for reynolds stress modeling: Application to plane channel and square duct flows. Computers & Fluids p 106246

[12] Cécora RD, Radespiel R, Eisfeld B, et al (2015) Differential reynolds-stress modeling for aeronautics. AIAA Journal 53(3):739–755

[13] Charpiat G, Girard N, Felardos L, et al (2019) Input similarity from the neural network perspective. In: Wallach H, Larochelle H, Beygelzimer A, et al
(eds) Advances in Neural Information Processing Systems, vol 32. Curran Associates, Inc., URL https://proceedings.neurips.cc/paper_files/paper/2019/file/c61f571dbd2fb949d3fe5ae1608dd48b-Paper.pdf

[14] Chen J, Hachem E, Viquerat J (2021) Graph neural networks for laminar flow prediction around random two-dimensional shapes. Physics of Fluids 33(12):123607

[15] Chiarini A, Quadrio M, Auteri F (2021) Linear stability of the steady flow past rectangular cylinders. Journal of Fluid Mechanics 929. https://doi.org/10.1017/jfm.2021.819

[16] Chiarini A, Quadrio M, Auteri F (2022) A new scaling for the flow instability past symmetric bluff bodies. Journal of Fluid Mechanics 936:R2. https://doi.org/10.1017/jfm.2022.99

[17] Donon B, Liu Z, Liu W, et al (2020) Deep statistical solvers. Advances in Neural Information Processing Systems 33:7910-7921

[18] Dupuy D, Odier N, Lapeyre C, et al (2023) Modeling the wall shear stress in large-eddy simulation using graph neural networks. Data-Centric Engineering 4:e7

[19] Duraisamy K, Iaccarino G, Xiao H (2019) Turbulence modeling in the age of data. Annual Review of Fluid Mechanics 51:357-377

[20] Eivazi H, Tahani M, Schlatter P, et al (2022) Physics-informed neural networks for solving reynolds-averaged navier-stokes equations. Physics of Fluids 34(7):075117

[21] Foures DPG, Dovetta N, Sipp D, et al (2014) A data-assimilation method for Reynolds-averaged Navier-Stokes-driven mean flow reconstruction. Journal of Fluid Mechanics 759:404-431

[22] Garnier P, Viquerat J, Rabault J, et al (2021) A review on deep reinforcement learning for fluid mechanics. Computers & Fluids 225:104973

[23] Giannetti F, Luchini P (2007) Structural sensitivity of the first instability of the cylinder wake. Journal of Fluid Mechanics 581:167-197

[24] Goodfellow I, Bengio Y, Courville A (2016) Deep learning. Massachusetts Institute of Technology

[25] Guégan T (2022) Contrôle d'écoulements turbulents par apprentissage automatisque. Université de Poitiers

[26] Gühring I, Raslan M, Kutyniok G (2020) Expressivity of deep neural networks. arXiv preprint arXiv:200704759 34

[27] Hamilton WL (2020) Graph representation learning. Synthesis Lectures on Artificial Intelligence and Machine Learning 14(3):1-159
[28] Lapeyre CJ, Misdariis A, Cazard N, et al (2019) Training convolutional neural networks to estimate turbulent sub-grid scale reaction rates. Combustion and Flame 203:255-264

[29] Lee S, You D (2019) Data-driven prediction of unsteady flow over a circular cylinder using deep learning. Journal of Fluid Mechanics 879:217-254

[30] Ling J, Templeton J (2015) Evaluation of machine learning algorithms for prediction of regions of high Reynolds-averaged Navier Stokes uncertainty. Physics of Fluids 27(8):085103

[31] Ling J, Kurzawski A, Templeton J (2016) Reynolds averaged turbulence modelling using deep neural networks with embedded invariance. Journal of Fluid Mechanics 807:155-166

[32] Lino M, Fotiadis S, Bharath AA, et al (2023) Current and emerging deep-learning methods for the simulation of fluid dynamics. Proceedings of the Royal Society A 479(2275):20230058

[33] Mantić Lugo V, Arratia C, Gallaire Fmc (2014) Self-consistent mean flow description of the nonlinear saturation of the vortex shedding in the cylinder wake. Phys Rev Lett 113:084501

[34] Martinkus K, Lucchi A, Perraudin N (2021) Scalable graph networks for particle simulations. URL https://arxiv.org/abs/2010.06948, 2010.06948

[35] McConkey R, Yee E, Lien FS (2022) On the generalizability of machine-learning-assisted anisotropy mappings for predictive turbulence modelling. International Journal of Computational Fluid Dynamics 36(7):555-577

[36] Morra P, Semeraro O, Henningson DS, et al (2019) On the relevance of reynolds stresses in resolvent analyses of turbulent wall-bounded flows. Journal of Fluid Mechanics 867:969-984

[37] Patel Y, Mons V, Marquet O, et al (2024) Turbulence model augmented physics-informed neural networks for mean-flow reconstruction. Physical Review Fluids 9(3):034605

[38] Pickering E, Rigas G, Schmidt OT, et al (2021) Optimal eddy viscosity for resolvent-based models of coherent structures in turbulent jets. Journal of Fluid Mechanics 917:A29

[39] Pope SB (1975) A more general effective-viscosity hypothesis. Journal of Fluid Mechanics 72(2):331-340

[40] Provansal M, Mathis C, Boyer L (1987) Benard-von Karman instability: transient and forced regimes. Journal of Fluid Mechanics 182:1-22. https://doi.org/
[41] Raissi M, Yazdani A, Karniadakis GE (2020) Hidden fluid mechanics: Learning velocity and pressure fields from flow visualizations. Science 367(6481):1026–1030

[42] Ren P, Xiao Y, Chang X, et al (2021) A survey of deep active learning. Association for Computing Machinery computing surveys 54(9):1–40

[43] von Saldern JG, Reumschüssel JM, Kaiser TL, et al (2023) Self-consistent closure modeling for linearized mean field methods. In: AIAA AVIATION 2023 Forum, p 4351

[44] Sanchez-Gonzalez A, Heess N, Springenberg JT, et al (2018) Graph networks as learnable physics engines for inference and control. In: International Conference on Machine Learning, pp 4470–4479

[45] Schmelzer M, Dwight RP, Cinnella P (2020) Discovery of algebraic Reynolds-stress models using sparse symbolic regression. Flow, Turbulence and Combustion 104:579–603

[46] Shukla K, Xu M, Trask N, et al (2022) Scalable algorithms for physics-informed neural and graph networks. Data-Centric Engineering 3

[47] Singh AP, Duraisamy K (2016) Using field inversion to quantify functional errors in turbulence closures. Physics of Fluids 28(4):045110

[48] SIPP D, LEBEDEV A (2007) Global stability of base and mean flows: a general approach and its applications to cylinder and open cavity flows. Journal of Fluid Mechanics 593:333–358. https://doi.org/10.1017/S0022112007008907

[49] Ströfer CA, Xiao H (2021) End-to-end differentiable learning of turbulence models from indirect observations. Theoretical and Applied Mechanics Letters 11(4):100280

[50] Thuerey N, Weißenow K, Prantl L, et al (2020) Deep learning methods for reynolds-averaged navier–stokes simulations of airfoil flows. AIAA Journal 58(1):25–36

[51] Toshev AP, Galletti G, Brandstetter J, et al (2023) Learning lagrangian fluid mechanics with e (3)-equivariant graph neural networks. In: International Conference on Geometric Science of Information, Springer, pp 332–341

[52] Vinuesa R, Brunton SL (2022) Enhancing computational fluid dynamics with machine learning. Nature Computational Science 2(6):358–366

[53] Volpiani PS, Meyer M, Franceschini L, et al (2021) Machine learning-augmented turbulence modeling for RANS simulations of massively separated flows. Physical Review Fluids 6(6):064607
[54] Wallin S, Johansson AV (2000) An explicit algebraic reynolds stress model for incompressible and compressible turbulent flows. Journal of Fluid Mechanics 403:89-132

[55] Weatheritt J, Sandberg R (2016) A novel evolutionary algorithm applied to algebraic modifications of the RANS stress-strain relationship. Journal of Computational Physics 325:22-37

[56] Wilcox DC, et al (1998) Turbulence modeling for CFD, vol 2. DCW industries La Canada, CA

[57] Yang S, Vinuesa R, Kang N (2024) Enhancing graph u-nets for mesh-agnostic spatio-temporal flow prediction. 2406.03789

[58] Zhao Y, Akolekar HD, Weatheritt J, et al (2020) Rans turbulence model development using CFD-driven machine learning. Journal of Computational Physics 411:109413

[59] Zhou Z, He G, Yang X (2021) Wall model based on neural networks for les of turbulent flows over periodic hills. Physical Review Fluids 6(5):054610

[60] Zienkiewicz O, Taylor R, Nithiarasu P (eds) (2014) The Finite Element Method for Fluid Dynamics, seventh Edition edn. Butterworth-Heinemann, Oxford, https://doi.org/10.1016/B978-1-85617-635-4.00016-9, URL https://www.sciencedirect.com/science/article/pii/B9781856176354000169