Learning Distributions of Complex Fluid Simulations with Diffusion Graph Networks
Authors: Mario Lino Valencia
Tobias Pfaff
Nils Thuerey
Venue: ICLR 2025 Oral (2025)
Keywords: Graph Neural Networks
Diffusion Models
Physics Simulations
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)
TL;DR: We propose an efficient graph-based latent diffusion model, which allows us to directly sample unsteady flow states from their equilibrium distribution given a mesh discretisation of the system and its physical parameters.

Abstract: Physical systems with complex unsteady dynamics, such as fluid flows, are often poorly represented by a single mean solution. For many practical applications, it is crucial to access the full distribution of possible states, from which relevant statistics (e.g., RMS and two-point correlations) can be derived. Here, we propose a graph-based latent diffusion model that enables direct sampling of states from their equilibrium distribution, given a mesh discretization of the system and its physical parameters. This allows for the efficient computation of flow statistics without running long and expensive numerical simulations. The graph-based structure enables operations on unstructured meshes, which is critical for representing complex geometries with spatially localized high gradients, while latent-space diffusion modeling with a multi-scale GNN allows for efficient learning and inference of entire distributions of solutions. A key finding of our work is that the proposed networks can accurately learn full distributions even when trained on incomplete data from relatively short simulations. We apply this method to a range of fluid dynamics tasks, such as predicting pressure distributions on 3D wing models in turbulent flow, demonstrating both accuracy and computational efficiency in challenging scenarios. The ability to directly sample accurate solutions, and capturing their diversity from short ground-truth simulations, is highly promising for complex scientific modeling tasks.

Published as a conference paper at ICLR 2025
LEARNING DISTRIBUTIONS OF COMPLEX FLUID SIMULATIONS WITH DIFFUSION GRAPH NETWORKS
Mario Lino' Tobias Pfaff
Nils Thuerey'
Technical University of Munich
2 Google DeepMind
ABSTRACT
Physical systems with complex unsteady dynamics, such fluid flows, are of- ten poorly represented by a single mean solution. For many practical applica- tions, crucial to access the full distribution of possible states, from which relevant statistics (e.g., RMS and two-point correlations) can be derived. Here, we propose a graph-based latent diffusion (or alternatively, flow-matching) model that enables direct sampling of states from their equilibrium distribution, given a mesh discretization of the system and its physical parameters. This allows for the efficient computation of flow statistics without running long and expensive numerical simulations. The graph-basc based structure enables operations on unstruc- tured meshes, which is critical for representing complex geometries with spatially localized high gradients, while latent- space diffusion modeling with a multi-scale GNN allows for efficient learning and inference of entire distributions of solutions. A key finding is that the proposed networks can accurately learn full distributions even when trained on incomplete data from relatively short simulations. We apply this method to a range of fluid dynamics tasks, such as predicting pressure dis- tributions on 3D wing models turbulent flow, demonstrating both accuracy and computational efficiency in challenging scenarios. The ability to directly sample accurate solutions, and capturing their diversity from short ground -truh simula- tions, is highly promising for complex scientific modeling tasks.
INTRODUCTION
Numerically solving partial differential equations is essential in many scientific and engineering fields (Karniadakis & Sherwin, 2013), with fluid dynamics being of particular interest across a wide range of disciplines (Verma et al., 2018; Kochkov et al., 2021). Nevertheless, the high computational cost often limits the practical application of these methods and the exploration of large parameter spaces. While advances in deep learning have produced promising surrogate models, these typically focus on predicting mean flows (Lino et al., 2023) or trajectories (Kim et al., 2019; Stachenfeld et al., 2021). For long trajectories, such methods often suffer from instability issues (Lippe et al., 2024). Besides, in many applications, the primary interest is not in individual trajectories or mean flows but in the statistical properties and probability distributions of unsteady flow fields in statistical cquilibrium (Pope. 2000). For instance, RMS fluctuations must be accounted for in the design of airfoils or when deciding the relative placement of aerodynamic components (Caros et al., 2022; Jané-Ippel et al., 2023), and turbulence research inherently revolves around flow field statistics.
While specialized techniques exist for certain applications such as unsteady Reynolds-averaged simulations, which approximate the mean flow (Alfonsi, 2009; Pope, 2004) obtaining accurate distributions and statistics in the general case still requires extensive, Often, long simulations are needed, and in turbulent flows, very small time-steps are necessary. Ad- ditionally, a significant portion of compute time spent during the warm-up phase, before reaching statistical equilibrium (Moin & Kim, 1982), and in complex systems, such weather simulations, reliable results require running ensembles of perturbed simulations (Maher et al., 2021).
connemmnnninng simulations.
To address this, we leverage diffusion model that directly learns the distribution of equilibrium flow states, enabling efficient sampling from the flow state distribution and allowing for the in- expensive computation of any desired distributional metrie. Unlike prevalent grid-based diffusion models (Song et al., 2021; Lienen et al., 2024), our approach employs a graph-based representation, operating directly meshes. This allows for the accurate representation of complex geometries

---

Published as a conference paper at ICLR 2025
and the adaptive allocation of resolution (Pfaff al., 2021) crucial advantage for complex 3D scenarios, which would otherwise require computationally expensive operations on a dense 3D grid.
However, integrating graph neural networks (GNNs) with diffusion models for large systems presents challenges, particularly due to limitations in the efficiency of message passing for prop- agating features across the graph, which can hinder global denoising. To address this, introduce an efficient multi-scale GNN architecture to model the denoising transitions and propose working on a compressed latent mesh. Operating in this latent space not only reduces inference time but also mitigates the introduction of undesired high-frequency noise in the solutions. Additionally, we find that the underlying distributions can be accurately learned from short simulations, even when these simulations lack sufficient diversity to fully represent the flow statistics. We demonstrate that our approach can accurately capture full distributions where other methods, such as Gaussian mixture models or variational autoencoders (VAEs), suffer from noise and mode collapse. In addition, it exhibits strong generalization and can be significantly more efficient than numerical simulations.
RELATED WORK
Numerous deep learning algorithms have been proposed to address scientific problems (Guest et al., 2018; Jumper et al., 2021; Pathak al., 2022), with fluid flow problems receiving considerable atten- tion within this field (Morton et al., 2018; Bar-Sinai et al., 2019; Lino et al., 2023). At the same time, probabilistic models, such as Gaussian mixture models and generative adversarial networks (GANs) (Goodfellow et al., 2014), have enabled the modeling of probability distributions over plausible states of physical systems based on collections of system snapshots (Maulik al., 2020; Drygala et al., 2022; Kim & Lee, 2020). Recently, denoising diffusion probabilistic models (DDPMs) have emerged as powerful class of probabilistic models, initially applied to image generation (Ho al., 2020; Nichol & Dhariwal, 2021), where GANs have struggled to capture the full diversity of the data distribution (Dhariwal & Nichol, 2021). These models have found significant utility in physical modeling well, with applications such as flow field super-resolution and reconstruction (Shu et al., 2023; Li et al., 2023), uncertainty estimation in under-resolved simulations (Liu & Thuerey, 2024), and enhancing the long-tern stability of autoregressive rollouts (Lippe et al., 2024; Rühling Cachay et al., 2024; Kohl et al., 2024). More closely related to our work, Gao et al. (2024) and Lienen et al. (2024) leveraged DDPMs model the probability distribution of fully-developed flows.
However, most of these studies relied on CNNs to parameterize the denoising process, which gen- erally restricts system representation to Cartesian grids. GNNs, on the other hand, offer greater flexibility and allow for the discretization of systems using unstructured meshes (Pfaff et al., 2021; Lino al., 2022). While GNNs have shown efficacy in modeling steady and unsteady simulations involving complex geometries and uneven spatial gradients (Belbute-Peres et al., 2020; Lam et al., 2023), their use combination with DDPMs remains relatively unexplored.
Recent efforts have applied DDPMs to graph structures for tasks such as molecule synthesis (Xu et al., 2022; Hoogeboom et al., 2022: Trippe et al., 2023; Vignac et al., 2023), protein generation (Wu et al., 2024), and inverse protein folding (Yi et al., 2024). More aligned with our work, Wen et al. (2023) combined DDPMs and GNNs to forecast the probability distribution of future system states. Their approach was applied to traffic and air quality prediction, where the system's 's evolution was conditioned on past states and the underlying graph structure. However, their work applied data compression only along the temporal dimension, neglecting along with previous work spatial compression. As a result, these approaches were restricted to small graphs. In contrast, our work targets systems with thousands of nodes in both 2D and 3D. This required modeling the denoising process by multi-scale GNNs. These typically follow a U-Net like architecture (Gao & Ji, 2019), with specialized pooling and unpooling layers (Lino et al., 2022). Additionally, by shifting the generative problem to a compressed latent space, we let the model focus on more meaningful features, while any residual noise in the generated samples is reduced when transforming them back to physical space. unlike previous GNN-based DDPMs. We also address promising, yet overlooked, topic: extrapolating the full distribution of system states from incomplete simulations.
LEARNING TO SAMPLE FROM DISTRIBUTIONS OF SIMULATIONS
We propose Diffusion Graph Networks (DGNs), a model that learns the probability distribution of dynamical states of physical systems, defined by their discretization mesh and their physical
2

---

Published as a conference paper at ICLR 2025
(b)
Training data vs. learned Std.
transient:
statisticca equilibrium stage
Training data Std.
PDF
Ground-truth Std.
Y()
2,00
R = 0.875
PDF system
Z-t0
7:00)
LDGN Std.
Learn full distrisution of 2, states from trajectories
R = 0.980
Figure I: (a) We learn the probability distribution of the systems' converged states provided only short trajectory of length d << 'I' per system. (b) A preview of our turbulent wing experiment. The distribution learned by our LDGN model accurately captures the variance of all states (bottom right), despite seeing only an incomplete distribution for each wing during training (top right).
parameters, by applying the DDPM framework to the mesh nodes. Additionally, we introduce a second model variant, the Latent DGN (LDGN), which operates in pre-trained semantic latent space rather than directly in the physical space.'
We represent the system's 's geometry using mesh with nodes Var and edges Eas. where each node i is located at æi. The system's 's state at time t, (t)), is defined by F continuous fields sampled at the mesh nodes: (t) := (t) E RF" | E V&a h. where (t(t) y(z, t). Simulators evolve the system through an infinite sequence of states, y
(ta). (t))). (t..).... (t.)).... }, starting from an initial state Y(4)). We assume that after an initial transient phase, the system reaches a statistical equilibrium (see Figure la). In this stage, statistical measures of Y. computed over sufficiently long time intervals, time-invariant, even if the dynamics display oscillatory or chaotic behavior (Wilcox, 1998). The states in the equilibrium stage, Z c y. depend only on the system's geometry and physical parameters, and not on its initial state.
In many engineering applications, such as aerodynamics and structural vibrations, the primary focus is not on each individual state along the trajectory. but rather the statistics that characterize the system's 's dynamics (e.g., mean, RMS, two-point correlations). However, simulating a trajectory of converged states 2 long enough to accurately capture these statistics be very computationally expensive, especially for real-world problems involving 3D chaotic systems. To address this, we aim directly sample converged states Ztt) Subsequently, we can analyze the system's 's dynamics by drawing multiple samples.
without simulating the initial transient phase.
(2,, Z2, ... ZNJ. our goal is learn probabilistic model of 3 that enables sampling of converged state Z(t) E Z, conditioned on the system's 's mesh, boundary conditions, and physical parameters. Crucially, this model must capture the underlying probability distributions even when trained on trajectories that are too short to fully characterize their individual statistics. Although this is an ill-posed problem, we hypothesize that, given sufficient training trajectories, it is possible to uncover their statistical correlations and shared patterns, enabling interpolation across the condition space.
Given a dataset of short trajectories from N systems,
3.1 DIFFUSION GRAPH NETWORKS
We use the DDPM framework (Ho et al., 2020; Nichol & Dhariwal, 2021) to generate states ((t) by denoising sample " R | drawn from an isotropic Gaussian distribution. The system's conditional information is encoded in a directed graph G := (V,E), where V = VM and the mesh edges Eι are represented as bi-directional graph edges E. Node attributes Ve [Df | E V} and edge attributes E. {ei (i.j) E E} encode the conditional features, including the relative positions between adjacent node, x; T,. Domain specifcc details on the node and edge encodings can be found in Appendix C and Table 4.
In the diffusion (or forward) process, node features from Z1 E V to Z E R X are vI 11--1 3,,1)
generated by sequentially adding Gaussian noise: Z--1 NiZZ;
"Code available http://giitab.commmmm.mphhggggg
3

---

Published as a conference paper at ICLR 2025
where B, (0, and Z Ztt). Any Z' can be sampled directly via:
Z' Vā, 2 + vI a.c.
with a, := -,, a, = IT.= 0₂ and e (), (Ho al., 2020). The denoising process re-
moves noise through learned Gaussian transitions: po 20 zn) 12)) =N(Z">3 : the mean and variance are parameterized as (Nichol & Dhariwal, 2021):
where
Mg 1-
Z
Br v1 ar E8
Eg exp (v% log 3r + (1- (1) log log B,).
(2)
Var
with 3 == (1 ā 11/-1/(11- 0, Here, Eg E R predicts the noise e in equation (1), and Va E R VixF interpolates between the two bounds of the process' entropy, 3, and B
DGNs predict Eg and Vgg using message message-passing based GNN (Battaglia al., 2018). This takes Z" as input, and it is conditioned on graph G, its node and edge features, and the diffusion step
1 1.- DGN, Z",, G.9 Ve, Ee,r).
We train the DGN using the loss function in equation (14). The full denoising process requires R evaluations of the DGN to transition from 2 R to zo, though more efficient sampling techniques exist (Nichol & Dhariwal, 2021; Song al., 2021).
GNN architecture (Sanchez- Gonzalez et al., 2020; Pfaff et al., 2021). In addition the node and edge encoders, our encoder includes a diffusion-step encoder, which generates a vector Temb E 1 Fun that embeds the diffusion step r. The node encoder processes the conditional node features ef, alongside Temb- Specifically, the diffusion-step encoder and the node encoder operate as follows:
The DGN follows the widely used picoddrrprroesssr-deeee
Temb & o LINEAR SINEMB(r),
+ LINEAR ([ő o LINEAR(U;) (Temb)), Vi V, (4)
where o denotes the activation function and SINEMB is the sinusoidal embedding function (Vaswani al., 2017). The edge encoder applies linear layer the conditional edge features eij. The encoded node and edge features are IR= dimeessoonl vectors (Femb = 4 x Fi). We condition each message-passing layer
by projecting Temb Fi-dimensional space and adding the result to the node features before each of these layers - ie., 0; + 0; (ANEAR(Femb)) Details on message passing can be found in Appendix B.I.
Previous work on graph-based based diffusion models has used sequential message passing to propagate node features across the graph (Hoogeboom et al., 2022; Wen al., 2023). However, this ap- proach fails for large-scale fluid flows, such as the ones studied in this paper, as denoising of global features becomes bottlenecked by the reach of message passing (Section 5). To address this, we adopt a multi-scale GNN for the processor, applying message passing on G and multiple coarsened versions of it in U-Net fashion (Lino et al., 2022; Cao al., 2023). This design leverages the U-Net's 's effectiveness in removing both high- and low-frequency noise (Si et al., 2023). To obtain each lower-resolution graph from its higher-resolution counterpart, we use Guillard's coarsening algorithm (Guillard, 1993), originally developed for fast mesh coarsening in CFD applications. As in the conventional U-Net, pooling and unpooling operations, now based on message passing, are used to transition between higher- and lower-resolution graphs. These operations detailed in Appendix B.3.
3.2 LATENT DIFFUSION GRAPH NETWORKS
Inspired by latent diffusion models for CNNs (Rombach et al., 2022), we extend the DGN frame- work to operate in a lower-dimensional graph-based representation that is perceptually equivalent to 3. This space is defined as the latent space of a Variational Graph Auto-Encoder (VGAE) trained to reconstruct (t). We refer to DGN trained this latent space Latent DGN (LDGN).
In this configuration, the VGAE captures small vortices), while the LDGN focuses on modeling mid- to large-scale patterns (e.g., the wake and vortex street). By decoupling these two tasks, we aim to simplify the generative learning pro- cess, allowing the LDGN to concentrate more meaningful latent representations that less
information (e.g.. spatial gradients and
4

---

Published as a conference paper at ICLR 2025
(a) VGAE
(b) LDGN inference
edge
Ve. E
Condition encader
conditoon
V,,
odgn
V E
Condition enoodar
Ocnaaiso latents
C
Phys coll
Latent noise
magnihidis at the rodees
LDGN
c
Node encuder
Iers
desociar
().1
magnitudees
Z'
F-R
P denniiinng alocs
Z'
Figure 2: (a) Our VGAE consists of a condition encoder, a (node) encoder, and (node) decoder. The multi-scale latent features from the condition encoder serve as conditioning inputs to both the encoder and the decoder. (b) During LDGN inference, Gaussian noise is sampled in the VGAE latent space and, after multiple denoising steps conditioned on the low-resolution outputs from the VGAE's condition encoder, transformed into the physical space by the VGAE's decoder.
sensitive to small-scale fluctuations. Additionally, during inference, the VGAE's decoder helps move residual noise from the samples generated by the LDGN. This approach significantly reduces sampling costs since the LDGN operates smaller graph rather than directly on G.
For our VGAE, we propose an encoder-decoder architecture with an additional condition encoder to handle conditioning imputs (Figure 2a). The condition encoder processes V. and Ec, encoding these into latent node features V and edge features E across I. graphs (g == (V 1 L)), where G g and the size of the graphs decreases progressively, i.e., | 12 [V4). This transformation begins by linearly projecting V. and E. to a Fas fedd-mmmmiional space and applying two message-passing layers to yield V. and E. Then, 1 encoding blocks are applied sequentially:
V E++1
+ MP MP o GRAPHPOOL (V- EE,, for | = ,,
(5)
where MP denotes message-passing layer and GRAPHPOOL denotes graph-pooling layer (see the diagram on Figure 9a).
The encoder produces two Fz FL-dimensional vectors for each node E V the mean μ4 and standard deviation o that parametrize a Gaussian distribution over the latent space. It takes as input a state (t), which linearly projected to a dim-iiinniio vector space and then passed through L 1 sequential down-sampling blocks (message passing + graph pooling), each conditioned on the outputs of the condition encoder (Figure 9b):
GRAPHPOOL MP MP (V LINEAR (V:). LINEAR (E¹)). for 1,2, L (6)
and a bottleneck block:
V < M MP (V LINEAR (V.!) LINEAR (E%))
The output features are passed through node-wise MLP that returns μ4; and a; for each node i E Vz. The latent variables are then computed Ç: BATCHNORM(J4; + orEe), where e; (O I). Finally, the decoder mirrors the encoder, employing a symmetric architecture (replacing graph pooling by graph unpooling layers) to upsample the latent features back the original graph G (Figure 9c). Its blocks are also conditioned on the outputs of the condition encoder. The message passing and the graph pooling and unpooling layers in the VGAE are the same as in the (L)DGN.
The VGAE is trained to reconstruct states Z(t) with KL-penalty towards a standard normal distribution on the learned latent space. Once trained, the LDGN can trained following the ap- proach in Section 3.1. However, the objective is now to learn the distribution of the latent states Ç, defined on the coarse graph g. conditioned on the outputs V and E' from the condition en- coder. As illustrated in Figure 2b, during inference, the condition encoder generates the conditioning features V and E: (for | 1, 2...., L). and after the LDGN completes its denoising steps, the decoder transforms the generated Ço back into the physical feature-space defined on G.
Unlike in conventional VGAEs (Kipf & Welling, 2016), the condition encoder is necessary because, inference time, need encoding of V. and E, graph -². where the LDGN operates. This
5

---

Published as a conference paper at ICLR 2025
(a) Pressure distribution
(b) Flow fields around an Ellipse
(c) Pressure distribution on wing
ellipse
plx,y.
y.x.y.))
v(x,y, 0
pls.y.
p6.. 2,
Figure 3: (L)DGNs can generate diverse states of dynamic mesh-based simulations given the mesh geometry and simulation parameters. We demonstrate this by learning: (a) the pressure ellipse in 2D laminar flow (ElLIPsE task), (b) the velocity and pressure fields resulting from that flow (ELLipsEFLow task), and (c) the pressure wing in 3D turbulent flow (WING task).
encoding cannot be directly generated by the encoder, as it also requires (t) as input, which is unavailable during inference. An alternative approach would to define the conditions directly in the coarse representation of the system provided by g but this representation lacks fine-grained details, leading to sub-optimal results (see Appendix D.3 for details).
EXPERIMENTAL DOMAINS
We consider three experimental domains: (i) pressure on the wall of an cllipse in 2D quasi-periodic laminar flow (ELLIPSE task), (ii) velocity and pressure fields around that ellipse (ELLiPsEFLOW task), and (iii) pressure on a wing in 3D turbulent flow (WING task). The ELLIPSE and ELLIPSEFLLOW data were adapted from Lino et al. (2022), while the dataset for the WING task was generated using Detached Eddy Simulation (DES) with OpenFOAM's s PISO solver.
The ELLIPSEFLOW task involves canonical fluid dynamics problem: predicting the velocity and pressure fields around an elliptical cylinder (Figure 3b). While this task benefits from spatial refine- ment near the surface of the ellipse, the ELLIPSE task fully leverages the graph-based representation by focusing solely on the surface of the immersed object (Figure 3a). The WING experiments target wings in 3D turbulent flow, characterized by intricate vortices that spontaneously form and dissipate on the wing surface (Figure 3c). This task is particularly challenging due to the -ig-ndimensional, chaotic nature of turbulence and its inherent multi-scale interactions across wide range of scales. The geometry of the wings varies in terms of relative thickness, taper ratio, sweep angle, and twist angle. These simulations are computationally expensive, and using GNNs allows us to concentrate computational effort on the wing's 's surface, avoiding the need for costly volumetric fields. A regular grid around the wing would require over 10y cells, in contrast to approximately 7,000 nodes for the surface mesh representation. The surface pressure be used to determine both the aerodynamic performance of the wing and its structural requirements. Fast access to the probabilistic distribution of these quantities would be highly valuable for aerodynamic modeling tasks.
RESULTS
Our main findings indicate that DGNs and LDGNs can generate high-quality fields and accurately reproduce the distribution of converged states, even when trained on incomplete distributions, Both models outperformed the baselines vanilla GNN, Bayesian GNN, Gaussian regression GNN, a Gaussian mixture GNN, and VGAE terms of sample accuracy and distributional aceuracy. LDGNs showed improvement over DGNs, particularly in distributional accuracy and suppressing undesired high-frequency noise.
Unless otherwise specified, all models were trained on 10 consecutive states per system for the ELLIPSE and ELLIPSEFLOW domains, and 250 consecutive states (shortly after the data- getta-genraaing simulator reached statistical equilibrium) for the WING domain. For the ELLIPSE and ELLIPSEFLOw domains, this covers approximately 26-48% % the time points required to capture a full vortex- shedding cycle, while for the WING domain, it represents about 10 % of the time points needed to

---

Published as a conference paper at ICLR 2025
(a) LDGN samples and matching G.T.
(b) G.T. and learned PDFs for 10-time-step trainining trajectories
Ground Truth
Vanil GNN
Bayesian GNN
GR-GNN
R 0.9993
6-3-
logae Ayguuah.
1048
IigoR 3.0501
WA 0.0275
3.7105
14 0.9872
0.2-
0.0
0.2
-0.4
-0.2
o0
07
0.4
046
16
-644
02.52
0.8
0.8
R 0.9958
GM-GNN
VGAE
DGN
LDGN
0.4
677
Wy-e 0.0534
0.33-
0.0222
6777
0.0142
0.22-
W5- 0.0138
14
-0.2
-0.4
6.12
0.22
112
0.2
0.
02
C6
Co
36
03
G.T. top wil G.T. bottom
Fred. Fred bottom
4.52
-.52
452
0.2
Figure 4: For a system on dataset INLIPSE-INDSTT, (a) samples from the LDGN and ground truth, and (b) probability density function from the DGN, LDGN, baseline models, and ground truth. The DGN and LDGN show the best distributional accuracy.
(a) Ground truth and predicted fields
(b) Ground truth and predicted distribution
Ground truth
Vanilla GNN
GM-GNN Are 5.593
Ground truth
GM-GNN
VGAF Af 0.9618
DGN
LDGN Aie 3-94
R 9999
DGN
LDGV
Figure 5: Samples from the DGN, LDGN, baseline models, and ground truth on (a) dataset BLLIPS&FLOWINNDISS and (b) dataset WING-TEST. The DGNs and LDGNs achieve the highest sample accuracy, with the DGNs showing good accuracy but retaining some nigh-frequency noise.
achieve statistically stationary variance. To assess the performance, we compared to ground-truth trajectories spanning three complete periods and time interval long enough to achieve stationary variance, respectively. To accelerate sampling, we adopted the strategy from Song & Ermon (2020), which considers only subset of the denoising steps. We used 50 steps, evenly distributed.
Sample accuracy Visually, the fields generated by our DGNs and LDGN models are faithful to the ground truth, as can be seen in Figures 4a and 5. In the ELLIPSE and ELLIPSEFLOW domains, it is possible to compare each sampled state with the states from a simulated trajectory and find its equivalent state (the one with the highest correlation). since these trajectories are smooth and periodic. For these domains, we report the coefficient of determination (R²) between the generated samples and their ground-truth equivalents as measure of the models' sample accuracy. In the ELLIPSE task, the LDGN achieved the highest accuracy, closely followed by the DGN (Table 5), and, in the ELLIPSEFLOW task, the DGN was the best performing model, followed by the LDGN (Table 6).
Although the DGN and LDGN models' accuracy is comparable, high-frequency noise is evident in the samples generated by the DGN model (Figure 5). To quantify high- fighffrequeccy noise, we computed the graph Fourier transform of the predicted fields and assessed the energy associated with the 10 highest graph eigenvalues for all systems in the FLL.PSEFFLOWIINDIIS and WING-TEST datasets. We found that the DGNs consistently over-estimated high-frequency content in generated samples, while LDGN samples remained close to the ground truth, with on average 3x reduced high- frequency eтrоr (See Appendix D.4 for detailed numbers). This discrepancy is due to the fact that DGNs required to handle high-frequency details directly, whereas LDGNs leverage VGAE to transform potentially noisy latent features back to physical space. The VGAE, being robust to noisy latent features, effectively mitigates such noise.
7

---

Published as a conference paper at ICLR 2025
Distributional accuracy High sample accuracy does not necessarily imply that model is learn- ing the true distribution. In fact, these properties often conflict. For instance, in VGAEs, the KL-divergence penalty allows control over whether to prioritize sample quality mode cover- age. To evaluate how well models capture the probability distribution of system states, we use the Wasserstein-2 distance. This metric be computed in two ways: (i) by treating the distribution at each node independently and averaging the result across all nodes, or (ii) by considering the joint distribution across all nodes in the graph. We denote these metrics as Wirde and W respec- tively. The node-level measure (W2 provides insights into how accurately the model estimates point-wise statistics, such as the mean and standard deviation at each node. However, it does not penalize inaccurate spatial correlations, whereas the graph-wise measure (W?Tmhh)) does. To ensure stable results when computing these metrics, the target distribution for the ELLIPSE and ELLIPSE- FLOW tasks is represented by 60 to 100 consecutive states, depending on the period length, while the predicted distribution is represented by 200 samples. For the WING task, the target distribution is represented by 2,500 consecutive states, and the predicted one by 3,000 samples.
In our experiments, the training trajectories are in- tentionally short. For instance, in the ELLIPSE-TRAIN and ELLIPSEFLOW--T OON-FIII datasets, the trajectories are too short to cover one full oscillation period, meaning they do not explicitly provide full statistical informa- tion about the systems. In the WING- TRAIN dataset, while the trajectories are long enough to capture the mean flow, they fall short of capturing the standard deviation, spatial correlations, or higher-order statis- tics. Despite these challenges, the DGN, and espe- cially the LDGN, are capable of accurately learning the complete probability distributions of the training trajectories and accurately generating new distribution for both in- and ditt-i--ittiitutionn physical settings. This capability is highlighted in Figure 7a, which shows the probability density function (PDF) for a training trajectory of the ELIJIPEE-TRAIS dataset (top right) and the PDF learned by the LDGN for the same system (bottom right). Similarly, Figure lb demonstrates how the LDGN accurately captures the standard deviation across all possible states (bottom right), even though this information is not provided in the training dataset. This excel- lent distributional accuracy was likewise observed for test data, Figures 4b and 6. As demonstrated in Figure 4b, this is rare property - (L)DGN is the only method which could faithfully represent the ground-truth distribution, while other methods show either collapse or wrong weighting.
Inferred Std. on test case
Ground truth
GM-GNN
R 0.602
DGN
LDGN
rC, 0.897
R = 0.958
Figure 6: The LDGN's s standard deviation is the closest to the ground-truth.
The Wasserstein-2 distances for all datasets are summarized Tables 7 to The LDGN consis- tently outperforms other models, followed by the DGN. Experiments in Appendix D.5 suggest that the ability of the diffusion- based models to learn the full distribution from short trajectories likely stems from their capacity to recognize common physical patterns across different systems, leverag- ing this knowledge to enhance its understanding of each individual trajectory. Unlike models that rely on sampling from compressed latent space or make assumptions about probability distribu- tions for parameters or outputs, their gradual denoising approach minimizes overfitting based on only the training data. This allows them extrapolate the training trajectories and generalize to new physical settings more effectively. The LDGN variant outperforms the DGN thanks to learning the distributions more meaningful latent space that is less sensitive to small-scale fluctuations.
Generalization We tested generalization -uist--diitributtion Reynolds numbers, ellipse thick- ness and ellipses rotated by 10-degrees in the ELLIPSE and ELLIPSEFLOW tasks. The DGN and LDGN models generally exhibit good generalization, shown by their sample and distributional accuracy in the out-of-distribution datasets, which is comparable to their performance on the in-distribution datasets (Tables 8). However, we observed that both types of accuracy degrade as the variance of the ground-trutt trajectories increases, even within the in- an-distribution datasets, as seen in Fig- ure 12. This is understandable, as higher diversity of states makes more difficult for the model to identify common features among systems that are close in parameter space, thereby affecting sam- ple quality. Additionally, propoabbayllity states in high-variance trajectories are rarely encountered during training, and can differ significantly from those in lower-variance trajectories, such as states of separated flow near the ellipse's 's leading edge. As a result, these states may be overlooked by the models, leading to reduced distributional accuracy.
8

---

Published as a conference paper at ICLR 2025
Influence of the number of states per training trajectory A key strength of the DGN and LDGN models is their ability to learn all possible states of a system from incomplete data. In the ELLIPSE task, we used 10 consecutive time-steps per training trajectory, while the full periods of these tra- jectories range from 21 to 38 time-steps. We found that the distributional accuracy of the DGN and LDGN models is only slightly affected by further reducing the length of the training trajectories, and, interestingly, this accuracy remains relatively close to the accuracy of training runs with com- plete periods. This behavior is illustrated in Figure 7b, which shows the W and W metrics for the DGN and LDGN models, as well as for four baseline models, on the ELL.PSEE-NNDST dataset. In this figure, the darker gray region indicates training trajectories that are always longer than the oscillation period, while in the lighter gray region at least one trajectory exists that is long enough. To the left of this area, it can be observed that the Wrisk of most baselines drops sharply due to the missing information. In contrast, the DGN and LDGN models are barely affected. It is also worth noting that the Gaussian mixture model (GM-GNN) achieves the lowest Whide among all models when trained on long trajectories. However, this baseline does not model spatial correlations, which explains its consistently high Wr values.
Computational efficiency Direct comparisons between very different implementations are inher- ently difficult. Nonetheless, to provide rough estimates of the performance characteristics for our WING experimental domain, the ground-trute tound-truth simulator, running on CPU threads, required 2,989 minutes to simulate the initial transient phase plus 2,500 equilibrium states sufficient to obtain a converged variance. In contrast, the LDGN model took only 49 minutes on 8 CPU threads and 2.43 minutes on a single GPU to generate 3,000 samples. If we consider the generation of a single converged state (for use as an initial condition in another simulator, for example), the speedup is four orders of magnitude on the CPU and five orders of magnitude on the GPU. Thanks to its latent space, the LDGN model is not only more accurate, but also 8x faster than the DGN model, while requiring only about 55% more training time. Performance details are available in Appendix D. This significant efficiency advantage suggests that these models could be particularly valuable in scenarios where computational costs otherwise prohibitive.
Multi-scale vs. single-scale DGN A key distinction of our GNN architecture is the use of a multi-scale GNN, while all previous works applied message passing directly the input graph (Hoogeboom et al., 2022; Wu et al., 2024; Yi et al., 2024; Wen et al., 2023). DDPMs are known to benefit from hierarchical and global denoising transitions, which enable the model capture and process information at multiple spatial resolutions (Dhariwal & Nichol, 2021; Si et al., 2023). While the previously used single-scale GNNs may suffice for smaller systems, we observed that they failed on larger systems, such as those in the ELLIPSEFLOW and WING domains, due to the limited reach of their receptive field. For fair comparison to heterogeneous previous work, we have considered a single-scale version of our DGN (with the same parametrization and hyperparameters) In the FLIIIPSE task (approximately 70 nodes per graph), our proposed four-scale DGN significantly outperformed its single-scale counterpart in both sample and distributional accuracy, as shown in
Training data PDF vs. Learned PDF
(b) Distributional
Training states
Training data Mytton 1.2856
(craph-wi wize aist. errar)
wyww (node dist.
full trajectory
0.08
C.8
3.6
LDGN
0.04
0.4
620
wython. in 0.2244
435
10F
101
10
102
107
10
0.1
# stales in treining trajectories les
states in trarning torjec
Vanilla GNN
GM-GNN VGAE
GR
LOGN
Figure 7: (a) The LDGN accurately learns the full probability distribution of the ELLIPSE system states, despite being trained on short trajectories (10 consecutive snapshots) that do not cover all possible states. (b) The distributional error of the DGN and LDGN models is not significantly affected by short training trajectories. They consistently exhibit the lowest Wagaaph distance.

---

Published as a conference paper at ICLR 2025
Table 11 of the appendix. In the ELLIPSEFLOW task (~2.3k nodes per graph) and the WING task (~6.8k nodes), the single-scale version did not produce feasible solutions (see Figure 15). We observed that, in the ELLIPSEFLoW task, performance improves as the number of scales is increased with the five-scale model yielding the best results, as also shown in Table 11. This model is also the most computationally efficient, outperforming the two-scale version by 25%
Comparison to probabilistic baselines We compare to a Bayesian GNN, Gaussian mixture model (GM-GNN), a Gaussian regression model (GR-GNN) and a VAE on the ELLIPSE task; to a GM-GNN and VAE on the ELLIPSEFLOW task; and to GM-GNN on the WING task.
Bayesian GNN: Although the Bayesian GNN provided reasonable predictions for both in- and out- of-distribution datasets, and achieved higher sample and distributional accuracy than other baselines, its performance was still substantially inferior to that of the DGN and LDGN models (see Tables and 7). Additionally, the probability distribution of the outputs incorrectly resembles a Gaussian distribution (Figure 5). Another significant drawback is its high training cost (approximately 8 times slower than the DGN training), making impractical for larger physical systems.
Gaussian Mixture Model: The GM-GNN represents learned distribution as a combination of mul- tiple Gaussian distributions. For each node it provides the mean, variance, and weight characterizing each Gaussian distribution. Once the trained GNN is evaluated for a given conditional input, sam- pling is efficient. However, since is performed independently for each node, the resulting samples are spatially discontinuous (Figure 5) and fail to represent the underlying distribution (Figure 4b). This leads to low sample accuracy (Tables 5 and 6) and high w:graph (Tables 7 to 9). Interestingly, in the ELLIPSE task, the GM-GNN outperforms the DGN and LDGN models in terms of Wyode, metric that does not penalize discontinuous samples, when trained on complete trajectories (grey region in Figure 7b). However, when trained on shorter trajectories, its performance clearly deteri- orates (white region in Figure 7b). This suggests that the GM-GNN is less capable at learning the shared physical features across different systems. In comparison, the DGN and LDGN models much better suited for learning from partial distributions, which is especially relevant for 3D chaotic systems due to the high cost of data generation. We also considered the GR-GNN special case of the GM-GNN. It represents the learned distribution as a single Gaussian. However, it generally exhibited a worse performance than the GM-GNN.
Variational Autoencoder: Finally, a VGAE was also considered as a baseline model. This model mirrors the architecture of the VGAE proposed for the LDGN (Section 3.2). but with the number of message- passinn and graph-pooling layers matching the totals used in the (L)DGN models. This baseline achieved a sample accuracy closer that of the DGN and LDGN models (Tables 5 and 6). In the ELIIPSE task, it successfully learned the distribution of states when trained on short trajectories (Figures 4b and 7b). However, in the more complex EI..IPSEFI.OW task, it suffered from mode collapse despite careful selection of the latent space size and the KL penalty (Table 8). We believe that training on short trajectories caused the decoder to become overly dependent on the input conditions rather than on the latent features. In contrast, the VGAE in the LDGN models did not encounter this issue, as their latent space is less compressed.
CONCLUSIONS
We have introduced physical- and latent-space baph-bbsed diffusion models that enable direct sam- pling of physical states in large, geometrically complex dynamical systems from their equilibrium distribution. This approach allows for the efficient computation of flow statistics without the need for long and expensive numerical simulations. Our latent-space model produces high-quality solutions without undesired high-frequency noise in a fraction of the time compared to the physical-space model. Notably, it accurately learns full distributions even when trained on incomplete data from short simulations, marking an important step towards more efficient probabilistic modeling. How- ever, our method has limitations that suggest interesting directions for future work. For example, the samples produced by our method are not temporally correlated, and it would be valuable to ex- tend the model with temporal conditioning. Additionally, incorporating physical constraints (e.g., through guidance) could further improve the quality of the generated solutions (Huang et al., 2024; Lino et al., 2024). Another promising direction is improving sampling speed, which could poten- tially be achieved through flow matching (Lipman et al., 2023) Appendix D.7. Despite these challenges, we believe our work represents a significant step towards leveraging diffusion models in real-world engineering applications.
10

---

Published as a conference paper at ICLR 2025
ACKNOWLEDGMENTS
M.L. and N.T. acknowledge the support of the European Research Council (ERC) Consolidator Grant SpaTe (No. CoG-2019-863850) The authors gratefully acknowledge the Gauss Centre for Supercomputing V. (www.gauss- cetesseeeteeeeeee eu) for providing computing time on the GCS Supercom- puter SuperMUC at the Leibniz Supercomputing Centre (www.Irz.de). M.L. also expresses gratitude to Qiang Liu for his valuable advice and discussions on diffusion models.
REFERENCES
Giancarlo Alfonsi.
Reynolds-averaged navier-stokes equations for turbulecce modeling. Applied Mechanics Reviews, 62(4):040802, 2009.
Jimmy Lei Ba, Jamie Ryan Kiros, and Geoffrey E Hinton. Layer normalization. arXiv preprint 1607:06007. 2016.
Yohai Bar-Sinai, Stephan Hoyer. Jason Hickey, and Michael P Brenner. Learning data-driven dis- cretizations for partial differential equations. Proceedings of the National Academy of Sciences, 116(31):15344-553) 2019.
Peter W Battaglia, Razvan Pascanu, Matthew Lai, Danilo Rezende, and Koray Kavukcuoglu. Inter- action networks for learning about objects, relations and physics. Advances Neural Information Processing Systems, 37, 2016.
Peter W Battaglia, Jessica B Hamrick, Victor Bapst, Alvaro Sanzhez- Vinicius Zambaldi, Mateusz Malinowski, Andrea Tacchetti, David Raposo, Adam Santoro, Ryan Faulkner, et al. Relational inductive biases, deep learning, and graph networks. crXhv:1806.01261, 2018.
Filipe de Avila Belbute-Peres, Thomas D. Economon, and J. Zico Kolter. Combining differentiable PDE solvers and graph neural networks for fluid flow prediction. In Proceedings of the 37th International Conference Machine Learning ICML 2020, 2020.
Chris D Cantwell, David Moxey, Andrew Comerford, Alessandro Bolis, Gabriele Rocco, Gianmarco Mengaldo, Danicle De Grazia, Sergey Yakovlev, J-E Lombard, Dirk Ekelschot, et al. Nektar++: An open-source spectral/hp element framework. Computer physics communications, 192:205- 219, 2015.
Yadi Cao, Menglei Chai, Minchen Li, and Chenfanfu Jiang. Efficient learning of mesh-based phys- ical simulation with bi-stride multi-scale graph neural network. In Proceedings of the 40th Inter- national Conference on Machine Learning, pp. 3541-3558. PMLR, 2023.
Lidia Caros, Oliver Buxton, Tsuyoshi Shigeta, Takayuki Nagata, Taku Nonomura, Keisuke Asai, and Peter Vincent. Direct numerical simulation of flow over a triangular airfoil under martian conditions. AIAA Journal, 60(7):3961-397) 2022. doi: 10.2514/1 00.2514/1.J061455
Prafulla Dhariwal and Alexander Nichol. Diffusion models beat gans on image synthesis. Advances in neural information processing systems, 34:8780-8794 2021.
Claudia Drygala, Benjamin Winhart, Francesca di Mare, and Hanno Gottschalk. Generative model- ing of turbulence. Physics of Fluids, 34(3), 2022.
Meire Fortunato, Tobias Pfaff, Peter Wirnsberger, Alexander Pritzel, and Peter Battaglia. Multiscale MeshGraphNets. In ICML 2022 Workshop on Al for Science, 2022.
Han Gao, Xu Han, Xiantao Fan, Luning Sun, Li-Ping Liu, Lian Duan, and Jian-Xun Wang. Bayesian conditional diffusion models for versatile spatiotemporal turbulence generation. Computer Meth- ods in Applied Mechanics and Engineering, 427:117033 2024.
Hongyang Gao and Shuiwang Ji. Graph U-Nets. In Proceedings of the 36th International Confer- on Machine Learning, pp. 2083-2092 PMLR, 2019.
Ian Goodfellow, Jean Pbudet-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farles Sherjil Ozair, Aaron Courville, and Yoshua Bengio. Generative adversarial networs. Advances in Neural Infor- mation Processing Systems, 27, 2014.
11

---

Published as a conference paper at ICLR 2025
Dan Guest, Kyle Cranmer, and Daniel Whiteson. Deep learning and its application to lhe physies. Anmual Review of Nuclear and Particle Science, 68(1)1111-88, 2018.
Hervé Guillard. Node-nested multi-grid method with Delaunay coarsening. Technical report, IN- RIA, 1993.
Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusion probabilistic models. Advances in neural information processing systems, 33:6840-688 2020.
Emiel Hoogcboom, Victor Garcia Satorras, Clément Vignac, and Max Welling. Equivariant dif- fusion for molecule generation in 3D. International conference on machine learning, pp- 8867-8887. PMLR, 2022.
Jiahe Huang. Guandao Yang. Zichen Wang, and Jeong Joon Park. DiffusionPDE: Generative PDE- Solving Under Partial Observation. arXiv preprint 1700::2:000..7663, 2024.
C Jané-Ippel, N Bempedelis, R Palacios, and S Laizet. High-fidelity simulations of wake-to- wake interaction in an atmospheric boundary layer over a complex terrain. In Journal of Physics: Conference Series, volume 2505, pp. 012033. IOP Publishing, 2023.
John Jumper, Richard Evans, Alexander Pritzel, Tim Green, Michael Figurnov, Olaf Ronneberger, Kathryn Tunyasuvunakool, Russ Bates, Augustin Zidek, Anna Potapenko, et al. Highly accurate protein structure prediction with alphafold. Nature, 596(7873):533-58) 2021.
George Karniadakis and Spencer Sherwin. Spectral/hp element methods for computational fluid dynamicy. Oxford University Press, 2013.
Byungsoo Kim, Vinicius C Azevedo, Nils Thuerey, Theodore Kim, Markus Gross, and Barbara Solenthaler. Deep fluids: A generative network for parameterized fluid simulations. In Computer Graphics Forum. Wiley Online Library, 2019.
Junhyuk Kim and Changhoon Lee. Deep unsupervised learning of turbulence for inflow generation at various reynolds numbers. Journal of Computational Physics, 406:109216 2020.
Durk P Kingma, Tim Salimans, and Max Welling. Variational dropout and the local reparameteri- zation trick. Advances neural information processing systems, 28, 2015.
Thomas N Kipf and Max Welling. 1611:160 2016.
Variational Graph Auto-Encoders.
arXiv preprint
Günter Klambauer, Thomas Unterthiner, Andreas Mayr, and Sepp Hochreiter. Self- neff-normalizing neural networks. Advances in Neural Information Processing Systems, 30, 2017.
Dmitrii Kochkov, Jamie A Smith, Ayya Alieva, Qing Wang, Michael P Brenner, and Stephan beaming-acceleratter Hoyer. Machine learning accelerated computational fluid dynamics. Proceedings of the National 2021.
Academy of Sciences, (18(1)):e10117811
Georg Kohl, Li-Wei Chen, and Nils Thuerey. Turbulent flow simulation using autoregressive condi- tional diffusion models. arXiv preprint arXiv:2309.01744 2024.
Remi Lam, Alvaro Sanchez-Gonzalez, Matthew Willson, Peter Wirnsberger, Meire Fortunato, Fer- ran Alet, Suman Ravuri, Timo Ewalds, Zach Eaton-Rosen, Weihua Hu, et al. Learning skillful medium-range global weather forecasting. Science, 382(6677)1141--422) 2023.
Tianyi Li, Alessandra S Lanotte, Michele Buzzicotti, Fabio Bonaccorso, and Luca Biferale. Multi- scale reconstruction of turbulent rotating flows with generative diffusion models. Atmosphere, 15 (1):60, 2023.
Marten Lienen, David Lüdke,, Jan Hansen- Palmus,, and Stephan Günnemann. From zero to tur- bulence: Generative modeling for 3d flow simulation. In Proceedings of the 12th International Conference on Learning Representations, 2024.
Mario Lino, Stathi Fotiadis, Anil Bharath, and Chris D Cantwell. Multi-scale rotation- -equiiin-qquvariiant graph neural networks for unsteady eulerian fluid dynamics. Physics of Fluids, 34(8), 2022.
12

---

Published as a conference paper at ICLR 2025
Mario Lino, Stathi Fotiadis, Anil A Bharath, and Chris D Cantwell. Current and emerging deep- learning methods for the simulation of fluid dynamics. Proceedings of the Royal Society A, 479 (2275):20330058, 2023.
Mario Lino, Tobias Pfaff, and Nils Thuerey. Squii-Equiaariant Diffusion Graph Nets: Synthesizing Flow Fields by Denoising Invariant Latents on Graphs. In ICML 2024 Al for Science Workshop, 2024.
Yaron Lipman, Ricky TQ Chen, Heli Ben Ben-Hamu, Maximilian Nickel, and Matthew Le. Flow match- ing for generative modeling. In In Proceedings of the 11th International Conference on Learning Representations, 2023.
Phillip Lippe, Bas Veeling, Paris Perdikaris, Richard Turner, and Johannes Brandstetter. Pde-refiner: Achieving accurate long rollouts with neural pde solvers. Advances in Neural Information Pro- cessing Systems, 36, 2024.
Qiang Liu and Nils Thuerey. Uncertainty-awar surrogate models for airfoil flow simulations with denoising diffusion probabilistic models. AIAA Journal, pp. 1-22, 2024.
Andreas Lugmayr, Martin Danelljan, Andres Romero, Fisher Yu, Radu Timofte, and Luc Van Gool. Repaint: Inpainting using denoising diffusion probabilistic models. In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pp. 11461-11471, 2022.
Nicola Maher, Sebastian Milinski, and Ralf Ludwig. Large ensemble climate model simulations: introduction, overview, and future prospects for utilising multiple types of large ensemble. Earth System Dynamics, 12(2):401-418. 2021.
Romit Maulik, Kai Fukami, Nesar Ramachandra, Koji Fukagata, and Kunihiko Taira. Probabilistic neural networks for fluid flow surrogate modeling and data recovery. Physical Review Fluids, 110):104401, 2020.
Parviz Moin and John Kim. Numerical investigation of turbulent channel flow. Journal of fluid mechanics, 11::341-377, 1982.
Jeremy Morton, Antony Jameson, Mykel J Kochenderfer, and Freddie Witherden. Deep dynamical modeling and control of unsteady fluid flows. In Advances in Neural Information Processing Systems, 2018
Alexander Quinn Nichol and Prafulla Dhariwal. Improved denoising diffusion probabilistic models. In International Conference on Machine Learning, pp. 8162-8171. PMLR, 2021.
OpenFOAM Foundation. OpenFOAM version 10. 2022. URL https://www.openfomm.comm
Jaideep Pathak, Shashank Subramanian, Peter Harrington, Sanjeev Raja, Ashesh Chattopadhyay, Morteza Mardani, Thorsten Kurth, David Hall, Zongyi Li, Kamyar Azizzadenesheli, et al. Four- castnet: global data-driven high-resolution weather model using adaptive fourier neural opera- tors. arXiv preprint arXiv:2202. /2 /2/4, 2022.
Tobias Pfaff, Meire Fortunato, Alvaro Banchez-Gonzaled and Peter W. Battaglia. Learning mesh- based simulation with graph networks. In 9th International Conference on Learning Representa- tions (ICLR 2021), 2021.
S. B. Pope. Ten questions concerning the large-eddy simulation of turbulent flows. New journal of Physics, 6(1):35, 2004.
Stephen B Pope. Turbulent flows. Cambridge University Press, 2000.
Robin Rombach, Andreas Blattmann, Dominik Lorenz, Patrick Esser, and Bjöm Ommer. High- resolution image synthesis with latent diffusion models. In Proceedings of the IEFF/CVF confer- ence on computer vision and pattern recognition, pp. 10684-10695, 2022.
Olaf Ronneberger, Philipp Fischer, and Thomas Brox. U-Net: Convolutional networks for biomedi- cal image segmentation. Proceedings of the International Conference on Medical image com- puting and computer-assisted intervention, pp. 234 241. Springer, 2015.
13

---

Published as a conference paper at ICLR 2025
Salva Rühling Cachay, Bo Zhao, Hailey Joren, and Rose Yu. Dyffusion: A dynamics-informed informed diffusion model for spatiotemporal forecasting, Advances in Neural Information Processing Sys- tems, 36, 2024.
Alvaro Ganchez-Gonzalez, Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovee, and Peter Battaglia. Learning to simulate complex physics with graph networks. In Proceedings of the 37th Interaatinnal Conference on Machine Learning, 2020.
Dule Shu, Zijie Li, and Amir Barati Farimani. physics-informed diffusion model for high-fidelity flow field reconstruction. Journal of Computational Physics, 478: 11972, 2023.
Chenyang Si, Ziqi Huang, Yuming Jiang, and Ziwei Liu. Freeu: Free lunch in diffusion u-net. arXiv preprint arXiv:2309. 1497, 2023.
Jiaming Song, Chenlin Meng, and Stefano Ermon. Denoising diffusion implicit models. In Pro- ceedings of the 9th International Conference on Learning Representatiovns, 2021.
Yang Song and Stefano Ermon. Improved techniques for training score-based generative models. Advances in neural information processing systems, 33: 12438-12445 2020.
Kimberly Stachenfeld, Drummond B Fielding. Dmitrii Kochkov, Miles Cranmer, Tobias Pfaff, Jonathan Godwin, Can Cui, Shirley Ho, Peter Battaglia, and Alvaro Gonzalez.. Learned coarse models for efficient turbulence simulation. 1:X1v:2112.15275. 2021.
Brian Trippe, Jason Yim, Doug Tischer, David Baker, Tamara Broderick, Regina Barzilay, and Tommi S Jaakkola. Diffusion Probabilistic Modeling of Protein Backbones in 3D for the motif- scaffolding problem. In In Proceedings of the 1 11th International Conference on Learning Repre- sentations, 2023.
Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Eukasz Kaiser, and Illia Polosukhin. Attention all you need. Advances in neural informa- fion processing systems, 30, 2017.
Siddhartha Verma, Guido Novati, and Petros Koumoutsakos. Efficient collective swimming by harnessing vortices through deep reinforcement learning. Proceedings of the National Academy doi: 1..0773/paa.180002333
of Sciences, 11542355544-55844 2018.
Clément Vignac, Igor Krawezuk, Antoine Siraudin, Bohan Wang, Volkan Cevher, and Pascal Frossard. Digress: Discrete denoising diffusion for graph generation. In Proceedings of the 11th International Conference on Learning Representations, 2023.
Haomin Wen, Youfang Lin, Yutong Xia, Huaiyu Wan, Qingsong Wen, Roger Zimmermann, and Yuxuan Liang. DiffSTG: Probabilistic spatio-temporal graph forecasting with denoising diffusion models. In Proceedings of the 31st ACM International Conference on Advances in Geographic Information Systems, pp. 1-12, 2023.
DC Wilcox. Turbulence Modeling for CFD. DCW industries, La Canada, 1998.
Kevin E Wu, Kevin K Yang, Rianne van den Berg, Sarah Alamdari, James Y Zou, Alex Lu, and Ava P Amini. Protein structure generation via folding diffusion. Nature communications, 15(1): 1059, 2024.
Minkai Xu, Lantao Yu, Yang Song, Chence Shi, Stefano Ermon, and Jian Tang. GeoDiff: A Geo- metric Diffusion Model for Molecular Conformation Generation. In In Proceedings of the 10th International Conference on Learning Representations, 2022.
Kai Yi, Bingxin Zhou, Yiqing Shen, Pietro Liò, and Yuguang Wang. Graph denoising diffusion for inverse protein folding. Advances in Neural Information Processing Systems, 36, 2024.
14

---

Published as a conference paper at ICLR 2025
A DENOISING DIFFUSION PROBABILISTIC MODELS (DDPMS)
DDPMs can learn to synthesise a sample =" from a data distribution ((=)) ) by progressively denois- ing sample z" drawn from an isotropic Gaussian distribution (Song & Ermon, 2020; Ho et al., 2020). This denoising process mirrors the inverse of a fixed Markov chain of length R, known as the forward process, and it consists of R denoising steps. The forward (or diffiusion) process produces variables 2 through z by sequentially injecting them with Gaussian noise:
g(="|=">> == (2";; VI -,,--1, ,,II)
where 3, E (0. Thanks the reparametrization trick (Kingma et al., 2015), it is possible to shortcut in the forward process and directly sample x", at any diffusion step r. according to
q<=1 = Nx Vā (1-6) (1-ā,))) ā,)1),
(9)
this is.
z" === Vā, + ā,E
(10)
Here, ar := Br, a. == IT 0₂ and c ~ (O I) (Ho et al., 2020). If R is large enough and the 3-schedule selected properly, the forward process leads to R (nearly) following a Gaussian distribution. Typical values for R are around thousands, and linear and cosine ß-schedules the most common (Nichol & Dhariwal, 2021).
Given the data distribution g(="), the transitions of the reverse (or denoising) process (==")1²²)) can be approximated using a neural network parametrised by 0,
Ea(z", rr). p=(="= | =- N' (="; po((", r),
(11)
Although the mean (= r) could be parametrised in several ways, Ho et al. (2020) found best to first predict the noise with the neural network and then compute pm :)) according to
po(2", .r)
Var 2 V- 6o <0(x".r)) .r)
_
3r
(12)
To simplify training, Ho et al. (2020) avoided learning the variance and fixed it to En(=", 1 with either
of = or of - 11 (1- á,),,,
which correspond to the upper and lower bounds on the reverse process entropy, respectively. Fol- lowing Nichol & Dhariwal (2021), we opt to parameterise Ez(z", r) as an interpolation between the lower, Br. and upper, B,,, bounds on the reverse process entropy in the log domain:
Eo(=", .)) e [o(=" x) log (1-(01v0(=",))) log 3.
(13)
where vo(2".1) r)) is predicted by the neural network. Nichol & Dhariwal (2021) demonstrated that learning En(2", .)) improved the likelihood of DDPMs and reduced the number of diffusion steps re- quired for sampling from the learned g( ((2)) distribution at inference, which we considered especially relevant for flow field synthesis due to the time constraints design optimisation tasks.
To train a network to approximate cog (2",, rr) and va(2", r), we minimise the loss function L, defined as the sum of two terms:
L Lsimple + Avhh
(14)
where Lximple and Lvbb are defined as
Csimple En,zo, [lle - [[lc-co(z",n)) cc((",r))111
(15)
T
Lup = log po (0 (2 + Σ DKL (q(x1-1 20) || po(2)-1120))
The term Lvbb represents the variational lower bound. Its first term is the negative log-likelihood of Gaussian distribution, while the remaining terms are the Kullback-Leible (KL) divergence between two Gaussians. The gradients of co(2", .r) are only backpropagated through Esimlle, whereas the Lyb term is used to optimise vo(=", r). This training strategy was proposed by Nichol & Dhariwal (2021) and has also been proven successful in subsequent work Dhariwal & Nichol (2021).
To reduce gradient noise, we employ importance sampling and dynamically weight each loss term based on the loss history (Nichol & Dhariwal, 2021).
15

---

Published as a conference paper at ICLR 2025
B ADDITIONAL MODEL DETAILS
The implementation of our models and baselines, including their weights, and demonstration scripts available at https://gihub..co//ung hotp://ooommubbmdmnccttt
B.1 (L)DGN ARCHITECTURE AND TRAINING
The message-passing layers the (L)DGN's processor follow the general framework described in Battaglia et al. (2016) and Battaglia et al. (2018). The edge- and node-update functions modeled by MLPs with one hidden layer, containing Fh neurons, and with F} output features. Their activation functions are SELU functions with their standard parameters (Klambauer et al., 2017) as are all the activations in the model. These MLPs are preceded by layer normalization (Ba et al., 2016). Specifically, the edge-update, edge-aggregation, and node-update steps follows:
Cij t Welis + MLP (LN ((ey)os/u))))
V(ijj E
(16)
6 Σ Cijy
Vj E V.
(17)
VEN;
Bj < Wou + MLP" (LN ((ee(l,;)))
Vj E,
(18)
The (L)DGN processes information at L length scales by creating rewou-resolution graphs and prop- agating the node and edge features across them. These low-resolution graphs have fewer nodes and edges than the original mesh graph, allowing a single pesssage-passing layer to propagate attributes over longer distances more efficiently. In the (L)DGN's ' processor, information is first distributed and processed the high-resolution graph through two sequential message- passinng layers. is then passed to the immediately lower-resolution graph via graph pooling (equations (21) and (20)), and, here, the features are again processed through two message- saatgggg-aiiingggg layers. This process repeated L- times in total. Once the lowest-resolution features processed at G they are passed back to the scale immediately above through a graph-unpooling layer (equation (22)), which also takes input the node-feature values this scale. The features successively passed through pairs of message-passing paassing layers and unpooling layers until the information processed again at the original mesh graph.
For a fair comparison all DGN, LDGN and baselines models in Appendices B.5, D.2 and D.3, trained for the same task, follow similar architecture, with the same number of message-passing layers and a comparable number of weights. In the ELLIPSE task, all models operate across 4 levels of resolution (14 message-passing layer in total) and have approximately 3.5M parameters. In the ELLJPSEFI.Ow task, they span 5 levels of resolution (18 pressage-passing layer in total) and have 4.5M parameters, while in the WING task, they cover 6 levels of resolution (22 message-passing layer in total) with 5.5M parameters. In the LDGN models, the VGAE applies two graph-pooling layers before denoising (Figure 8), and the denoising model applies the remaining graph-pooling layers until the target resolution is achieved. In the ELLIPSE task (ID meshes), this VGAE provides approximately 4-fold compression, whereas in the ELLIPSEFLOW and WING tasks (2D meshes), roughly offers 16-fold compression.
2x graph pooling in VGAE
physical-space mesh
latent-space mesh
Figure 8: In our experiments, the LDGN models are applied to the latent space of two-scale VGAEs.
The DGN and LDGN models were trained using the node- wise mean of the loss function described in equation (14). To reduce gradient noise and improve training stability, employed the impor-
16

---

Published as a conference paper at ICLR 2025
tance sampling technique from Nichol & Dhariwal (2021). During each training epoch, single state was sampled from each training trajectory.
The initial learning rate set to 10 and was reduced by factor of 10 when the training loss plateaued for n consecutive epochs. In the ELLIPSE and ELLIPSEFLOW tasks, we used n 50, while in the WING task, n set 250 due to the shorter length of the training dataset. Training continued until the learning rate reached 10 This training strategy was consistently applied to both the DGN, LDGN, and baseline models.
B.2 VGAE ARCHITECTURE AND TRAINING
The LDGN's VGAE consists of condition encoder, a node encoder, and node decoder, as illus- trated in Figure 9 and outlined in Section 3.2. It was trained to reconstruct the input node features with low-weighted KL term between the latent distribution of each node, Çs, and a standard normal distribution:
(- 2
LVAAE =
Σ LEV
2/1?
+
+10
--10
6
(-217
Σ iEV
(1
log (o:) - - -o,)
(19)
V
(a) VGAE's condition encoder
(b) VGAE's node encoder
(c) VGAE's node decoder
x
E!
Unzar
2 MP
Uinar Z
Linear
MP
MP
graph pooling
graph pooling
graph unpooling
2
V
E
V
Unears
2
V2 E
2 MP
graph pooling
graph pooling
graph unpooling
Linear
2 MP
3
MLP. f4: a;
2
MP
Limoaraa
Figure 9: Diagram of the VGAE architecture. It consist of (a) condition encoder, (b) node encoder, and (c) a node decoder.
B.3 GRAPH POOLING AND UNPOOLING
Our DGN, LDGN, and baseline models (except for the VGAE, since there are no skip connections between the encoder and decoder branches) follow a U-Net-like architecture (Ronneberger et al., 2015), where message passing is applied across L levels of resolution (Gao & Ji, 2019; Lino et al., 2022: Fortunato et al., 2022; Cao al., 2023). Thus, in addition to the primal graph, G = (V, E), L - 1 lower-resolution graphs are employed. Here, for each resolution level & from 2 to L, we denote their graphs as g == (Ve. E²), while the primal graph is denoted as G = g' := ( To generate each lower-resnlttion graph G", we apply Guillard's coarsening algorithm (Guillard, 1993), originally developed for multi-grid methods in CFD. This algorithm removes subset of nodes from the input mesh (or bi-directional graph) as outlined in Algorithm For ID meshes, the compression ratio is approximately two, and for 2D meshes, around four.
Unlike the methods in Gao & Ji (2019) and Lino et al. (2022), our coarsening approach preserves the original relative distribution of nodes in the coarser graphs. Furthermore, for the WING task, the voxel-grid coarsening used in Lino et al. (2022) introduces new nodes and edges that do not lie on the WING surface. In contrast, Guillard's coarsening ensures that the coarser nodes and edges remain strictly on the Wing surface, thereby avoiding undesired connections between nodes that may be spatially close in the world-space but distant the mesh-space (e.g., a pair of nodes located on the upper and lower surfaces of the Wing). Additionally, unlike the method in Fortunato et al. (2022), our approach operates automatically without relying on an external mesh generator. Another valid and similar alternative could have been bi-stride pooling (Cao et al., 2023).
In our implementation, once De+1 C ve is obtained from G², each node E De is assigned a parent node vi+1, denoted P, Each node i is assigned the parent node that is the fewest
17

---

Published as a conference paper at ICLR 2025
Algorithm Guillard's ss coarsening algorithm (Guillard, 1993)
1: mask < (V
Vector of size [V9 filled with ones
2: for node i V do if maskk[ then
Iterate node-by-node
3:
If first visit to node then this node is not dropped
4:
for node j E No do mask(y] <- 0 end for end
5:
The incoming neighbours are dropped
6:
7:
8:
end for
9: y++ t [mask]
hops away. multiple candidates are equidistant, the parent node chosen based on the smallest Euclidean distance. Nodes that share the same parent node E ye+1 are referred to as its child nodes, denoted by Chy C V The set of edges gh+ is constructed to preserve the connectivity of the child nodes in ye. Specifically, there is an edge from E vé+1 to j y if at least one edge exists between their child nodes in V All these operations are performed as pre-processing steps and do not occur during model evaluation.
To implement a U-Net-like architecture, we use graph pooling and unpooling transition between higher and lower-resolution graphs. Following Lino et al. (2022) and Fortunato et al. (2022), both pooling and unpooling operations are based on message passing.
Graph pooling Graph pooling from V to D++1 1 performed by message passing from each node
i E Chy to its parent node j E yf+1, following equation (20). Edge features g++1 are assigned via linear projection of the relative positions between parent nodes, as shown in equation (21).
Σ
MLP (LN ([Linear(æ) xx)));))),
Vj D4+1
(20)
KEChy
Gjk + Linear( -xk),
V(,k) E g(+)
(21)
Graph unpooling Similarly, graph unpooling from Ve to D is performed by message passing to each node i V from its parent node P. V as defined in equation (22):
Uy - MLP (LN ([Linear(x; x - (u, (x])),
Vi De-1 and j=P,
(22)
The features of g-- come via a skip connection from the encoder branch, just before the graph- pooling layer from ve-1 to Ve applied.
B.4 TREATMENT DIRICHLET BOUNDARY CONDITIONS
In the ELLIPSEFLLOw task, the velocity values at the inlet and the ellipse wall known priori, defining Dirichlet boundary conditions. To ensure that these boundary conditions are ex- actly satisfied in the DGN's samples, we apply deterministic diffusion transitions the nodes on these boundaries. Specifically, for boundary nodes, we propagate the values deterministically as of = Va, :. and instead of learning their denoising transitions, we directly reverse the diffu- sion process, using o.- = /
During training, the loss function is not evaluated at these boundary nodes, and during inference, their noisy node features are sampled not from Gaussian distribution but are instead assigned their "diffused" values, i.e., d.e Van vf. While more sophisticated methods exist for handling bound- ary conditions, such as those proposed by Lugmayr et al. (2022), our approach remains efficient and straightforward.
In the LDGN and baseline models, we directly replace the output at these boundary nodes with their right values.
18

---

Published as a conference paper at ICLR 2025
Table
Model size and myyper-paammeee
Task
Model
# Parameters
Fi
Femb
Fi
Fr.
# Scales
DGN
3.51
512
LDGN
3.52 M
128
512
126
1
2 2
Vanilla GNN
3.53
154
ELLIPSE
Bayesian GNN
3.56 M
138
4
GR-GNN
3.52
154
GM-GNN
3.53 M
154
4
VGAE
3.48
DGN
4.50 M
128
512
5
LDGN
4.51 M
128
512
126
2+3
ELLIPSEFLOW
Vanilla GNN
4,50 M
153
5
GM-GNN
4.50
153
VGAE
4.51 M
126
32
5
DGN
5.48 M
128
512
6
WING
LDGN
5.49 M
128
512
126
2 +
GM-GNN
5.45 M
152
6
B.5 BASELINE DETAILS
We compared the DGN and LDGN models with several GNN-based probabilistic baseline models: a Vanilla GNN, Bayesian GNN, a Gaussian Regression GNN, Gaussian Mixture GNN, and VGAE. Details these baselines are provided below.
The so-called Vanilla GNN is deterministic GNN that mirrors the architecture of the DGN, except for the absence of the diffusion-step encoder and the Temb conditioning layers. It takes the conditioning features V. and E. as input and returns constant (for the given conditions) node features, 2/. which represent the mean learned from the training data. It was trained to maxi- mize the expected data, this is
Vanilla GNN
1
Lvanilla
E1 2/11 M
(23)
iEV
where z; are the states sampled from the training trajectories. Vanilla GNNs were trained for the ELLIPSE and ELLIPSEFLOW tasks. Since it was trained on short trajectories (10 snapshots per trajec- tory), its predictions resembles single state rather than the mean of the complete distribution. This explains the relatively high sample accuracy of these models, while their distributional accuracy is, of course, low (Tables to 8).
Bayesian GNN
In Bayesian neural networks, prediction variability arises from modeling the pa- rameters random variables with associated probability distributions. We compared our DGN and LDGN models to Bayesian GNN on the ELLIPSE task. To train the Bayesian GNN efficiently, used variational inference and approximated the posterior distribution of the parameters with Gaussian. The architecture of the Bayesian GNN mirrors that the DGN, except the absence of the diffusion-step encoder and the Temb conditioning layers. During training, Evidence Lower Bound (ELBO), assuming a Gaussian prior distribution over the parameters. We performed grid search to find the optimal weight for the KL-divergence term in the ELBO loss (Liu & Thuerey, 2024), which was determined to be 0.1 for maximizing distributional accuracy on the ELLIPSE-INNIST dataset.
maximized the
Gaussian Mixture GNN (GM-GNN) A Gaussian mixture model (GMM) represents learned distribution as combination multiple Gaussian distributions (Maulik et al., 2020). Its architec- ture mirrors that of the DGN. Through grid search, selected three Gaussian components for the
19

---

Published as a conference paper at ICLR 2025
ELLIPSE and ELLIPSEFLOW tasks, and five components for the WING task. Once the GNN is evalu- ated for a given conditional input, sampling is efficient: a Gaussian distribution is selected based on the returned weights, followed by sampling from that distribution.
Gaussian Regression GNN (GR-GNN) We also considered Gaussian Regression model, re- ferred to as GR-GNN, which is special case of the GM-GNN that represents the learned distribu- tion single Gaussian.
VGAE As VGAE baseline, we considered a model mirroring the architecture of the VGAE proposed for the LDGN, but with the number of message-passing and graph-pooling layers matching the total used in the DGN models. For the ELLIPSE task, we set the size of the latent features to FL and performed grid search to find the optimal weight for the KL term (Kingma et al., 2015), which was determined to be 0.001 to maximize distributional accuracy on the ELLIPSE-INDIST dataset. For the ELI.IPFFI..w task, found the best KI. weight be 10-8, and the optimal size of the latent features is Fz- = 32 to maximize distributional accuracy on the ELL.PSSEFL.O-INDDSS dataset.
The number of hidden features per node and edge after encoding is denoted as FR. and the size of the diffusion-step embedding represented by Femb. In the LDGNs' VGAE, the intermediate hidden features are denoted by Fae. and the number of latent features by FL. The values of these hyperparameters for each model, as well as the total model size, are summarized in Table
C EXPERIMENTAL AND DATASET DETAILS
The datasets used in our ELLIIPSFFI.W and ELLIPSE experiments are available at https//hugginggaace cofdatastt/mmariolnoo//liii and the datasets used in our WING experiments are available at https: //huggingface.cmommm //buggingaae..o/daasse/mmiiiiii wooddatseee/mmarionnoo//ii
To the best of our knowledge, no previous studies have tackled large-scale, unsteady, and turbulent dynamics on detailed unstructured meshes, featured in our WING datasets. This motivated the creation of these new datasets to provide a realistic and challenging application context.
All the simulations in our experiments are governed by the incompressible Javier-Stoker equations, which describe the motion of fluids and capture the balance of mass and momentum in incompress- ible flows. In their non-dimensional form, these equations are expressed as follows:
Continuity equation (mass conservation):
Ou
dv
+
+
Oy
0z
Momentum conservation equations along the x, y. and z directions, respectively:
Ou
Ou
Ou
Ou
8
0 u
0Lu
+
+
Ot
dr
Oy
0z
dx
Re
Ox
0y2
0z
0v
du
de
dv
0p
8pv
00
+ u
+ 11
+ w
+
Of
Dr
Dy
0z
Oy
Re
Oz
Oy²
023
dw
Ow
Ow Oy
Ow
Op
8Pw
8
8Pw)
+ u
+ v
+ w
Øt
dr
0z
Re
Ox.2
Oy²
022
Here, u(t, z). v(t, y, 2), and w(1, 2r, y. z) represent the velocity components along the y, and z directions, respectively; p(t, v. y, 2) denotes the pressure, and Re the Reynolds number, dimensionless quantity measuring the ratio of inertial to viscous forces, which characterizes the flow regime (Pope, 2000).
Since these equations lack an analytical solution, understanding of fluid behavior under differ- ent conditions relies heavily on experiments and numerical simulations. The flow regime can be classified based the Reynolds number. For low Reynolds numbers, the flow is laminar, exhibit- ing smooth and predictable layers with minimal disruption. For high Reynolds numbers, the flow becomes turbulent, characterized by chaotic, irregular fluctuations. Although the Navier-Stokes
20

---

Published as a conference paper at ICLR 2025
equations are deterministic. turbulent flows chaotic due to the high nonlinearity of these equa- tions, which renders the system sensitive to even minor perturbations. Despite this chaos, in fully developed turbulence, the energy introduced (e.g., from shear pressure gradients) is balanced by energy dissipation through viscous forces, allowing the turbulence statistics to reach a steady state over time (Wilcox, 1998).
ELLIPSEFLOW task This task involves a canonical fluid dynamics problem, inferring the velocity and pressure p fields in simulations of flow around an elliptical cylinder where the free-stream locity is parallel to the ellipse's major axis (Figure 3b). The flow is laminar and unsteady, exhibiting periodic vortex shedding due the formation of unstable low-pressure wake behind the elliptical cylinder (Figure 3b). We compute probability density functions (PDFs) and statistical measures by sampling states within a periodic cycle.
In the training dataset (ELLIPSEFLOW-TRE TRAIN), the Reynolds number (Rt) ranges from 500 to 1000, and the ellipse's 's relative thickness ranges from 0.5 to 0.8. Test datasets evaluate model perfor- mance for in-distribution Re and thickness TELLIPSEFLOW-INDISSS Re (ELLInsEFFLOW-LLWE RE and ELLIPSEFLOW-HHHHRE: respectively) and lower and higher thickness (ELL.PSFFLOWWTHI and ELLIPSEFFI. THO-KKKKK respectively). These datasets are low-resolution ver- sions of the datasets from Lino et al. (2022). Here, each simulation consists of 101 consecutive time-steps of fully developed 2D flow around ellipse. The original simulations were generated using Nektar++ (Cantwell et al., 2015). We reduced the resolution by approximately a factor of two. The ellipse's ' major axis length is one (non-dimensional unit), and the domain height ranges from five to six. The training dataset contains 5000 simulations, while each test dataset comprises 50 simulations sampled from the original datasets in Lino et al. (2022). Additionally, the test dataset contains 24 simulations for ellipses with an angle-of-attack (the angle between the ellipse's 's major axis and the free-stream velocity) of 10 degrees. Parameter values for each dataset are detailed in Table 2.
as well as for lower and higher
ELLIPSEFLOWWAG
To condition probabilistic models, each node encodes the Reynolds number and a one-hot vector w; indicating its type (inlet, wall, or inner node), while each edge (i,j) encodes the relative position between nodes, Cij := Zy Iy. The output is the velocity u. ut(tt) and pressure Ps(t") each node i and at an arbitrary time t
Table 2: Values of the parameters for the ELLIPSE andELLIPSEFLON systems in the training and test datasets
Reynolds number R
Minor axis b
Angle attack (deg)
Average # nodes ELLIPSE
Average nodes FI..IPSFFL.w
Dataset
Purpose
BLLIPSE(FLOO)-TRAIN
[500, 1000)
[0.5 0.8]
71
2340
Training
BLLIISS(PLO)-INNDSSTT
[500, 1000]
(..5
2328
BLLIPSE(FLO)--LoWRE
[400, 500]
[0.5 0.8]
2376
Test
ELLIPSE(FLOO)-HHHHR::
[1000, 1100]
[0.5 0.8]
2312
Test
BLLIPsSEFFLOW)-THII
[500, 1000)
[0.45 0.5]
2286
Test
BLIIPSE(FLOW))THIIC
[500, 1000)
[0.8 0.9)
76
2282
Test
ELLIPSE(FLOW)-AOA
[500. 1000]
[0.5 0.8]
10
74
2514
Test
ELLIPSE task It involves inferring pressure p on the wall of the ellipse from the previous task
(Figure 3a). The pressure distribution on the ellipses in each ELLIPSEFLOw dataset was extracted to create these datasets. The Reynolds number is provided as an input condition for the probabilistic models at each node, and the relative position between adjacent nodes, along with the projection of the free-stream velocity along the edges are provided as input conditions for each edge. The models output the pressure p. (l*) at each node i. The distribution of simulation parameters (Re and relative thickness) in the ELLIPSE datasets mirrors that in the ELLIPSEFLOW datasets (Table 2).
In both ELLIPSE and ELLIPSEFLOW simulations, the temporal variance of the system's ''s state increases primarily with the relative thickness the ellipses due to the formation of larger vortices and, to lesser extent, with Re. The effect of these parameters temporal variance is shown in Figure 12a.
21

---

Published as a conference paper at ICLR 2025
WING task These experiments involve wing in turbulent flow (Re 2 x 10%), characterized by 3D vortices that spontaneously form and dissipate at various locations on the wing surface (Figure 3c). Each node encodes unit outer normal vector îu and each edge (i.j) encodes G. = 2 Zi, Uso tij. the n4. Uso biyl. where by is unit vector parallel to the edge and by = big x ú. The output is (c()"). Unlike the previous systems, this flow chaotic, and its time statistics converge
long time spans.
We generated the training and test datasets for the WING systems using the PISO solver of OpenFOAM with the Spalart-Allmaras Delayed Detached Eddy Simulation turbulence model (OpenFOAM Founda- tion. 2022). The simulation and dataset meshes were generated with mappyHexMeshh The free-stream velocity is 100 km/h, and the kinematic viscosity is 1.5 x10-5 m²/s. A symmetry boundary condition was used on the vertical plane intersecting the wing root. The wing sections are 24XX NACA airfoils with constant relative thickness (indicated by XX), and the quarter-chord is horizontal (i.e., zero dihe- dral angle). The wing has length of and root chord length of m, with a root-section angle- of-attack of 20 degrees. The chord length decreases linearly toward the tip, and the angle of attack also varies linearly toward the tip. The geometric param- eters of the wing that vary across the training and test datasets include the relative thickness, the taper ratio (ratio between the tip and root chords), the sweep an- gle (angle between the quarter-chord line and a line perpendicular to the wing root), and the twist angle (angle of attack at the tip minus the angle attack the root). The values of these parameters for each dataset are listed in Table 3. Time-wise variance pri- marily depends the wing's is relative thickness, with lower values resulting in higher variance due to more detached flow. Variance also increases with larger twist angles and smaller sweep angles.
0.032
0.030
Cummulative Std.
0.028
0.026
0.022
0.020
Ground-truth states Predicted states
0.018
500
1000
1500 2000 2500 3000 # states
Figure 10: Node-wise mean of the time/sample-wise
standard deviation computed over an increasing number of states. In the training trajectories, only the first 250 time-steps after the transient stage are used, and their standard deviation is not close to its ground-truth value, which computed over 2,500 time-steps.
The training dataset consists of 1,000 simulations, and the test dataset consists of 16 simulations. We simulated the first 4 seconds (transient regime) starting from a RANS solution, then recorded the pressure on the wing every 0.002 seconds for 0.5 seconds (250 time-steps)
the training dataset, and for 5 seconds (2,500 time-steps) in the test dataset. The first 250 time-steps are not sufficient represent the variability of the states, as illustrated in Figure 10, but we aim to reduce the high cost of running these simulations by relying on the probabilistic models' ability to learn common patterns from an ensemble of simulations.
Table 3: Values of the parameters for the WING systems in the training and test datasets
Dataset
Thickness
Taper ratio
Sweep (degrees)
Twist (degrees)
time steps
Average # nodes
Purpose
WING-TRAIN
[10, 14]
[0.3, 0.7]
[20, 40]
[-5, 5]
250
6852
Training
WING-TEST
11, 13
0.4, 0.6
25, 35
-2.5, 2.5
2500
6828
Testing
The conditional features provided as input, as well as the predicted output for each system, are listed in Table 4.
D SUPPLEMENTARY RESULTS
Tables and 6 collect our measures of sample accuracy, the coefficient of determination, on the test datasets; and Tables 7 9 collect our measures of distributional accuracy, the graph-level
22

---

Published as a conference paper at ICLR 2025
Table 4: Conditional inputs and predicted outputs for each system type
Node conditions
Node outputs
Edge conditions
System
C;j
ELLIPSE
Re
has
Pt
ELL.PSEFILOW
Re, W.
37 Zii
Ui V/s Pt
WING
Wi. Usx
Lij. Uno
Tais. Uno bis
Pt
Vanilla GNN
GR-GNN
GM-GNN
VGAE
DGN
LDGN
0.0806
[ystum 0.6914
IV 9.3344
0220
0 3714
22
03
02
4.5
422
For a system in dataset ELI.IPSE-INDDST, probability density function from the DGN, LDGN, and baseline models trained on 100-time-step trajectories (full distribution). Ground truth in Figure 4b.
Figure 11:
Wasserstein-2 distance. The correlation between the variance of the ground-truth trajectories and the accuracy of the LDGN's predictions is illustrated in Figure 12. Figures 13 and 14 showcase ad- ditional examples of results generated by DGN, LDGN, and baseline models for the ELLIPSEFLOw and WING tasks.
(a) de in test sets
(b) R² (sample accuracy)
(c) Wgrash (distribution accuracy)
Pressure distribution ellipse
1100-
0.20
1100
1.000
1100
1000
900-
0.15
1000
0.995
1000
0.8
820
0.6
0.990
700-
0.10
0.985
600
0.05
500
0.5
0.5
0.7
0.9
0.9
0.5
0.8
0.9
0.980
400
0.5
0.6
0.7
0.8
0.9
1300-
0.06
11.00
0.98
1100-
around an ellipse
1000-
1000
1000-
Flow Fields
900
900
0.97
900
800-
0.04
800
300
700.
702
700
0.03
0.96
0.02
0.6
0.7
0.7 0.8
0.5
0.7
Figure 12: For eachach simulation the ELLIPSE (top row) and ELLIPSEFLOW (bottom row) test datasets with a O-degree angle-of-attack: (a) node-wise mean of the time-wise standard deviation of the ground truth fields. (b) average R² between the LDGN outputs and their corresponding ground truth, and (c) graph-wise W₂, distance between the learned distribution of the fields and their ground truth. Higher standard deviation results in more difficulty generating accurate fields and learning their probability distribution.
23

---

Published as a conference paper at ICLR 2025
Table 5: Mean coefficient of determination (R2) on the ELLIPSE datasets
ELLIPSE
-INDIST
-Low Re
-HIGRRE
-THIN
-THICK
-AOA
Model
Vanilla GNN
0.991
0.992 0.003
0.964 0.025
0.948 0.029
Bayesian GNN
0,962 0.025
0.967 0.026
0.952 0.021
0.952 0.047
0.945 0.020
0.931 0,050
GR-GNN
0,944 0.034
0.971 0.018
0.902 0.039
0.966 0.017
0.900 0.033
0.907 0.031
GM-GNN
0.951 0.030
0.973 0.019
0.912 0.033
0.970 0.013
0.901 0.032
0,922 0.031
VGAE
0,981 # 0.038
0.989 0.026
0.971 0.047
0.996 0.006
0.947 = 0.077
0.951 = 0.049
DGN
0.994 + 0.006
0.997 + 0.001
0.988 0.015
0.994 0.002
0.992 0.007
0.968 0.026
LDGN
0.995 = 0.007
0.998 = 0.002
0.986 a 0.019
0.997 = 0.001
0.991 = 0.009
0.966 = 0.028
Table 6: Mean coefficient of determination (R²) on the ELL.PSEFI.Ow datasets
ELLIPSE
FLOw
-INDIST
-LowRe
-HIGHRE
-THIN
-THICK
-AOA
Model
Vanilla GNN
0.979 + 0.015
0.983 0.011
0.972 0.016
0.984 0.007
0.962 = 0.022
0.972 + 0.019
GM-GNN
0.954 0.022
0.964 0.012
0.940 0.021
0.967 0.008
0.923 0.032
0.947 + 0.024
VGAE
0.981 0.014
0.983 0.012
0.974 0.016
0.984 0.008
0.966 0.021
0.976 0.014
0.007
0.989 0.009
0.991 0.005
0.987 0.014
LDIGN
0.986 0.011
0.981 0.016
D.1 PERFORMANCE
Providing a truly hardware-agnosse performance comparison is challenging. Most of the neural network components in our models are optimized for GPUs, while official OpenFOAM distributions lack GPU support. Additionally, performance results can vary depending on the target statistic (e.g., mean flow requires shorter simulations and less sampling compared RMS calculations).
In Table 10, we present breakdown of the runtime for the DGN and LDGN models on our WING experimental domain. The runtimes were measured on a CPU, limited to threads, and on a sin- gle RTX 3080 GPU. The time for estimating the distribution accounts for the generation of 3,000 samples. This number of samples is sufficient for the RMS of the pressure to converge. When running the GPU, maximized concurrency to achieve 100% GPU utilization (50 to 60 concur- rent samples). The ground -trutt simulation was performed using OpenFOAM''s PISO solver on the same workstation, limited to 8 CPU threads. Excluding pre-pot proccescssing and post-processing steps, sece-decond simulation (required for the RMS to become statistically stationary) takes 50 hours.
Table 7: Mean |ii-dimensional (i.e., graph-wise) Wasserstein-2 distance (Wyrpph) on the ELLIPSE datasets
ELLIPSE
-INDIST
-LoWRe
-HIGHRE
-THIN
-THICK
-AoA
Model
Vanilla GNN
0.96 0.52
0.70
1.32 + 0.49
0.35 0.10
2.15 0.27
1.31
Bayesian GNN
0.83 + 0.35
0.66 + 0.25
1.05 + 0.34
0.49 + 0.04
1.87 + 0.22
1.28 * 0.26
GR-GNN
1.09 + 0.55
0.79 + 0.42
1.49 + 0.51
0.44 + 0.12
2.32 + 0.24
1.51 * 0.53
GM-GNN
1.07 0.52
0.77
1.44 1 0.51
2.29
1.38 0.46
VGAE
0.44 +: 0.25
0.32 4: 0.21
0.59 + 0.21
0.13 - 0.03
1.13 4 0.20
0.80 + 0.16
DGN
0.29 + 0.15
0.21 + 0.09
0.42 + 0.18
0.16 0.02
0.56 + 0.14
0.58 0.1
LDGN
0.23 0.12
0.17 0.08
0.42 0.18
0.10 0.02
0.57 0.15
0.59 0.11
24

---

Published as a conference paper at ICLR 2025
Table 8: Mean Vlldimensional (i.e., graph-wise) Wasserstein-2 distance (W2mrh))
the ELLIPSE-
FLOW datasets
ELLIPSE
FLOw
-INDIST
-LowRe
-HIGHRE
-THIN
-THICK
-A0A
Model
Vanilla GNN
6.23 = 1.93
5.87 # 1.61
6.88 = 1.80
4.10 + 0.68
9.73 1.49
7.302 1.700
GM-GNN
7.24 + 2.12
6.76 + 1.67
8.29 + 2.04
5.05 * 0.83
10.77 * 1.69
8.393 1.767
VGAE
6.10
5.60 +: 1.54
6.70 1.72
4,00 0.68
9.36 1.51
7.077 1.607
DGN
4.72 2.10
4.04 + 1.74
5.48 = 2.01
3.20 =
7.76 2.39
5.964 2.125
LDGN
3.07 # 0.93
2.53 + 0.71
3.84 + 1.24
2.81 + 0.59
3.62 0.91
3.709 0.864
Table
Mean |V]-dimensional (i.e., graph-wise) Wasserstein-2 distance (Wggrph on the WING-TEST dataset
Metrie
GM-GNN
DGN
LDGN
W W
4.32 0.86
2.12 + 0.90
1.95 + 0.89
Compared to this baseline, we observe a one to two two three order-of-magniin redgn-of-manniuue speedup on the GPU. For single-sample generation, the speedup is indicated with respect to the time- timestepping runtime of the numerical solver. This speedup would appear even more significant if we also accounted for the simulation the transition regime.
speedup on the CPU and
The LDGN model is 8x more efficient than the DGN due to performing the denoising process in compressed latent space. However, this efficiency is only observed when the CPU/GPU cores are fully utilized.
D.2 INFLUENCE OF THE NUMBER OF SCALES IN DGN MODELS
A key distinction of our DGN architecture, compared to previous GNN-based DDPM models, is the of a U-Net-like architecture with graph pooling and unpooling layers specifically designed for mesh graphs (see Appendix B.3). Previous models were limited to small graphs, where applying message passing directly on the input graph was sufficient (Xu al., 2022; Hoogeboom al., 2022; Trippe et al., 2023; Vignac et al., 2023; Wu et al., 2024; Yi et al., 2024; Wen et al., 2023). However, DDPMs are known to require hierarchical and global denoising transitions, as these allow the model to capture and process information across multiple spatial resolutions (Dhariwal & Nichol, 2021; Si al., 2023). Below, we compare our DGN with a single-scale model. For a fair comparison to heterogeneous previous work, we have considered single-scale version of our DGN, with the
(a) OOD test case with b
(b) OOD test case
1100
Ground
Vanilla
GM-GNN
Ground truth
Vanilla GNN
GM-GNN
0.9307
0.9E10
R" 0.550
DGN
LDGN
VGAE
DGN
LDGN Af 8.940
0.9095
0.9468
0.9974
0.9958
Figure 13: Samples from the DGN, LDGN, baseline models, and ground truth for (a) a simulation from dataset ELLIPSEFLOW-THICS and (b) simulation from dataset GLLIPSEFLOW-HIGHHE The DGN and LDGN achieve the highest sample accuracy, with the DGN showing good accuracy but retaining some nigh-frequency noise.
25

---

Published as a conference paper at ICLR 2025
(a) NACA2411, taper 0.6, sweep 35 deg, twist 2.5 deg
(b) NACA2413, taper 0.4, sweep 25 deg, bwist 2.5 deg
Ground truth
GM-GNN
Ground truth
GM-GNN
Figure 14:
Samples from the DGN, LDGN, GM-GNN, and ground truth for two simulation from the WING-TEST dataset.
Table 10:
Performance comparison for sample generation and distribution estimation on the WING task, along with the speedup achieved over the numerical solver
CPU min/distribution (speedup)
GPU min/distribution (speedup)
CPU s/sample
GPU s/sample
Model
DGN
6.81
(x8.8)
0.59
(x1)2)
LDGN
0.98
49 (x61)
0.20
2.43 (x 1235)
same denoising parametrization, number and formulation of message-passe passsing layers, diffusion-step conditioning, and training settings.
In the ELLIPSE task (with approximately 70 nodes per graph), four-scale DGN significantly outper- formed its single-scale counterpart (with the same total number of pessaag-passing layers) both sample and distributional accuracy, as shown Table 11. In the ELLIPSEFLOw task (approximately 2.3k nodes per graph) and the WING task (approximately 6.8k nodes per graph), the single-scale DGN failed to produce feasible solutions, as illustrated Figure 15. Instead, the model learned transition noise close to zero, resulting in the amplification of the initial Gaussian noise during the denoising process, as described by equation (2). In the ELLIPSEFLOW task, performance improved as the number of scales increased (while maintaining the same total number of pessage-passing lay- ers), reaching optimal results with the five-scale model, as also shown in Table 11. This five-scale DGN was also the most computationally efficient: it was 25% faster than the two-scale model, 12% faster than the three-scale model, and 7% faster than the four-scale model.
Sample accuracy (R²), distributional accuracy (W:aph)) and inference time for GNN mod- els with varying numbers of scales
Table 11:
Dataset
# scales
Rta
Inference time ms /sample
in model
0.964 0.043
0.60
204
ELLIPSE-INDIST
4
0.995 0.006
0.29 0.15
160
2
0.957 0.026
5.37 1.34
436
3
0.966 0.034
5,48 + 2.04
376
ELLIPSEFLOWIINDA
4
0.990 0.010
5.01 +.10
353
5
0.990 0.010
4.72 + 2.11
328
26

---

Published as a conference paper at ICLR 2025
(a) Denoising with single-scale model
R 1000
R = 500
R = =
(b) Samples from single-scale multi--cael models
1-scale model
2-scale model
3-scale model
4-scale model
5-scale model
-478125.600
0.9958
Figure 15: Single-scale models failed to learn the ELLIPSEFLOW and WING tasks due to the limited
range of their receptive fields. (a) Instead of removing noise between transitions, single-scale model amplified the initial Gaussian noise.
(b) Samples from DGN models with one-, two-, three-, four-, and five-scale architectures for simulation from the ELLIPSEFLOW-INDIS dataset.
D.3 ALTERNATIVE VGAE ARCHITECTURE
Our VGAE incorporates condition encoder block (Figure 9a) that generates encoding of Va and E. the latent graph G which is then passed as conditioning input to the LDGN model. These encodings cannot be directly generated by the VGAE's node encoder (Figure 9b), as the node encoder also requires Zt) as input, which is unavailable during inference since we directly sample Gaussian noise in the latent space.
An alternative approach representation of the system provided by G To test this, we compared our LDGN with another version that lacks the condition encoder in the VGAE (but has similar number of weights), where node and edge conditional features are defined directly on G We evaluated this the ELLIPSE task and found that our proposed LDGN architecture provides better sample and distributional accuracy, as shown in Tables 12 and 13. This improvement is likely because our conditional encodings. although compressed, still retain information about relevant high-frequency features.
could be define the LDGN conditions directly within the coarse
Table 12: Comparison of the coefficient of determination (1i2) for an LDGN using our VGAE and one using the alternative VGAE
ELLIPSE
INDIST
-LowRe
-HIGHRE
-THIN
-THICK
Model
Alternative
0.993 + 0.009
0.996 + 0.003
0.985 + 0.017
0.994 + 0.004
0.990 0.010
Ours
0.996 = 0.007
0.998 = 0.002
0.988 = 0.018
0.998 0.001
0.993 0.008
Table 13: Comparison of the graph-wise Wasserstein-2 distance (Wgraph) for LDGN using VGAE and one using the alternative VGAE
ELLIPSE
-INDIST
-LowRe
-HIGHRE
-THIN
-THICK
Model
Alternative
0.22 0.08
0.44 0.17
0.16 + 0.03
0.60 + 0.16
0.62 0.09
Ours
0.17 0.07
0.41 + 0.17
0.10 + 0.02
0.54 + 0.16
0.57 + 0.10
27

---

Published as a conference paper at ICLR 2025
D.4 RESULTS FOR HIGH-FREQUENT ERROR ANALYSIS
To quantify the high-frequency noise present in the DGN and LDGN samples, we computed the graph Fourier transform of the predicted fields and assessed the energy associated with the 10 highest graph eigenvalues for all systems in the ELLIPSEFLOW-INDE as a measure of the high-frequency noise present in the predicted fields. The enror of this high- frequency energy relative to its ground truth, in percent, is collected in Table 14.
and WING- TEST datasets. This serves
We note that the high-frequency noise in the samples generated by DGN models is consistently higher than in the ground-truth fields. In contrast, the noise LDGN samples differs only slightly.
Table 14: Relative error (%) between the high-frequency energy in the predicted samples and the ground truth
Dataset Magnitude
DGN
LDGN
7
3.38
1.87
ELLIPSEFLOW-INDDE
V
8.97
3.49
p
11.09
3.11
WING-TEST
P
2.84
0.97
D.5 EFFICIENCY IN LEVERAGING TRAINING TRAJECTORIES
While the results in Section demonstrate the efficacy of our approach, there is no theoretical guar- antee that a model can learn full distributions from collection of short trajectories. This ability relies on the model's capacity to recognize common physical patterns across different trajectories and leverage this knowledge enhance its understanding of each individual trajectory. Such inter- polative capacity depends on two key factors: (i) having sufficient number of training trajectories enable smooth interpolation across the parameter space, and (ii) the model's ability effectively capture and generalize these dependencies.
To disentangle the effects the number of training tra- jectories and the model's design, we conducted additional experiments on the ELLIPSE task. Specifically, we trained DGNs, LDGNs, and baseline models on datasets contain- ing varying numbers of short trajectories, ranging from 500 to 5,000 samples. As shown in Figure 16, the dis- tributional accuracy the DGN and LDGN models im-
Wgraph (graph-wise dist. error)
1.0
GM-GNN
VGAE
0.8
DGN
LDGN
proves with an increasing number of training trajecto- ries, demonstrating their ability to interpolate effectively across trajectories associated with different input condi- tions and accurately reproduce the underlying distribu- tion. In contrast, the perfoomiin baseline model (VGAE) plateaued in performance beyond 2,000 trajecto- ries, highlighting its limited interpolation capacity despite having the same number of learnable parameters.
0.6
0.4
1000 2000 3000 4000 5000 Size of the training dataset (N)
The DGN, and especially the LDGN, models more efficient than the VGAE and GM-GNN at learn- ing common features from collection of short trajectories, resulting in signifi- cantly lower distributional error in large datasets.
Figure 16:
This advantage likely stems from the diffusion- based approach of DGNs and LDGNs, which, unlike models that rely on sampling from compressed latent space (e.g., VGAE) or impose assumptions about the underly- ing probability distributions (e.g.. GM-GNN), minimizes overfitting by gradually denoising. This approach en- ables them to infer missing dynamics by leveraging sim- ilar states encountered during training. Furthermore, the LDGN variant outperforms the DGN by learning distributions compressed latent space that is less sensitive small-scale fluctuations, resulting in superior interpolation performance.
28

---

Published as a conference paper at ICLR 2025
Regarding sampling accuracy, verified that it remains unaffected by the number of training tra- jectories. This behavior aligns with our expectations, as accurate samples can still be generated even without the ability to complete the system's 's PDF by feature interpolation.
D.6 FLOW STATISTICS DERIVED FROM LEARNED DISTRIBUTIONS
With our DGN and LDGN models, we aim to learn the probability density function (PDF) of fully- developed unsteady flows. In certain cases, however, learning the mean flow alone using deter- ministic GNN trained directly on mean-flow data might suffice. Nevertheless, such deterministic models may encounter challenges when the available training trajectories are too short to properly capture the mean flow. Furthermore, the distributions generated by the DGN and LDGN models can be used to derive other relevant statisties, such as the turbulent kinetic energy (TKE) or the Reynolds shear stress.
To evaluate the advantage provided by DGN and LDGN models over deterministic baselines in the mean-flow prediction task, we trained deterministic model, referred to as the Mean-Flow flow fields from the ELLIPSEFLOW-TRE TRAIN dataset. We compared
GNN, on wetme-step-averaaeeeee its performance to that of the models discussed in Section 5: the Vanilla GNN (another determin- istic model trained to maximize the likelihood; see Section B.5), the DGN, the LDGN, and other probabilistic baselines. The coefficient of determination for mean-flow predictions, evaluated the ELLIPSEFLLOW test datasets, is summarized in Table 15. Notably, the LDGN significantly out- performs both the deterministic and probabílistic baselines in predicting the mean flow, achieving coefficient of determination values close to one despite being trained on only 10 consecutive states per trajectory.
Figure 17 further illustrates these findings. Predictions from the two deterministic baselines (Mean- Flow GNN and Vanilla GNN) exhibit diffused vertex distributions corresponding the 10-time-step (training simulation length) average fields. This diffusion indicates that during training, these models are overly biased toward the conditions of individual samples and fail to generalize the true average flow across training samples. The DGN slightly mitigates this problem through its gradual denoising approach, which minimizes overfitting to individual training samples. The improvement more evident in the LDGN, which, thanks to its more meaningful latent space and latent input conditions, better leverages the information present across systems that are close in parameter space.
(a) Mean field for an In-distribution test case
Ground truth
Mean-Flow GNN
Vanilla GNN
DGN
LDGN
(b) Mean u field for generalisation test case with AoA 10 deg
Ground truth
Mean-Flow GNN
Vanilla GNN
DGN
LDGN
0.7735
0.8504
0.9263
R 0.9705
Mean field predicted by the Mean-Flow GNN and the Vanilla GNN (deterministic GNNs), as well as obtained from the distributions of the DGN and LDGN, for a test case from (a) dataset FLLIPPSFFL.OW-IND
Figure 17:
and (b) dataset ELL.PSEFL.Ow-AAA.
To further demonstrate the effectiveness the proposed diffusion models, we computed the turbu- lent kinetic energy (TKE). This quantity is defined as:
TKE ((u'v (u'n') (v'v'))
where (.) denotes ensemble averaging, and u/ == u (w) and o' = e- (u) are velocity fluctuations. The coefficients of determination for the TKE presented in Table 16. Reflecting the superior distributional accuracy the LDGN, this model significantly outperformed the bascline models, with the DGN emerging as the second-best performer. The LDGN's s advantage in these predictions
29

---

Published as a conference paper at ICLR 2025
Table 15: Mean coefficient of determination (R²) of the mean-flow on the ELLIPSEFLOW datasets
ELLIPSE
FLOW
INDIST
-LOWRR
-HIGHRE
-THIN
-THICK
-AoA
Model
Mean Flown-Flowwao GNN
0,868 0.014
0.872 0.012
0.865 0.015
0.883 0.002
0.843 0.008
0.863 0.013
Vanilla GNN
0.869 0.013
0.871 0.012
0.864
0.014
0.883 0.003
0.845 0.009
0.864 0.013
GM-GNN
0.867 0.015
0.870 + 0.012
0.857 0.016
0.883 0.003
0.841 0.009
0.860 = 0.014
VGAE
0.926 0.014
0.931 0.012
0.921 0.015
0.942 0.003
0.902 0.011
0.921 0.0155
DGN
0.941 a 0.020
0.946 0.015
0.936 0.021
0.954 = 0.006
0.912 = 0.028
0.930 = 0.025
LDGN
0.976 0.004
0.978 0.002
0.975 0.005
0.978
0.001
0.974 = 0.003
0.975 = 0.003
further illustrated in Figure 18. It is important to highlight that all these models were trained on only 10 consecutive time steps from each training simulation (see Section 5), making the results even more remarkable.
Ground truth
GM-GNN R 0.5552
LDGN
R 0.9007
0.9552
(a) 1/2 (suhe> + <rv's)
Ground truth
GM-GNN
DGN
LDGN
(b) Reynols shear stress <u V
R -0,0018
R 0.7822
R 0.8025
Figure 18: (a) Turbulent kinetic energy and (b) Reynolds shear stress obtained from the distributions predicted by the GM-GNN, DGN, and LDGN for test case from the BLLIPSEFLLOW-INDIS dataset.
Table 16: Mean coefficient determination (R²) of the TKE on ELLIPSEFLOW datasets
ELLIPSE
FLOW
-INDIST
-LowRe
-HIGHRE
-THIN
-THICK
-AoA
Model
GM-GNN
0.669 0.141
0.675 1 0.155
0.646 = 0.120
0.881 # 0.043
0.393 = 0.091
0.564 # 0.2
VGAB
-0.201 + 0.040
-0.182 + 0.040
-0.218 = 0.050
-0.145 = 0.011
-0,253 = 0.038
-0.197 0.043
DGN
0.783 0.215
0.874 0.141
0.755 0.241
0.863
0.485 0.334
0.654 0.279
LDGN
0.911 0.226
0,947 0.064
0,876 0.234
0.895 0.071
0.901 +: 0.046
0.895 0.156
D.7
FAST SAMPLING WITH FLOW MATCHING
Flow matching (Lipman et al., 2023) has recently emerged as promising generative model frame- work, offering significant efficiency gains over diffusion models by defining direct linear mapping between samples of Gaussian distribution (or any other distribution, in theory) and the target dis- tribution. This enables sampling to be performed with fewer steps.
Our DGN and LDGN architectures can be seamlessly adapted to the flow-matching training frame- work, benefiting from faster sampling. We observed that flow-matching GNNs (FM-GNNs) and latent-flow -ontttii--mttching GNNs (LFM-GNNs) outperform their diffusion- based counterparts when the number of denoising steps is limited to 10 fewer. However, for ~0 more denoising steps, diffusion models demonstrate superior performance.
This behavior is illustrated in Figure 19a, which compares the standard deviation of the pressure predicted by the DGN and the FM-GNN with 10 denoising steps for a sample from the WING- TEST dataset (ground truth shown in Figure 6). The shift in performance for 20 more denoising steps is
30

---

Published as a conference paper at ICLR 2025
shown in Figure 19b, which presents the sample accuracy and distributional error function of the number denoising steps for the DGN, FM-GNN, and their latent variants on the ELLIPSE-INDIST dataset.
From these results, we observe that the LFM-GNN achieves compelling balance between accuracy and sampling speed. However, further experimentation is necessary to fully explore its potential, which we leave for future work.
(a) 10-step denoising diffusion vs. flow matching
(b) Accuracy vs. # denoising steps in EllipseVal dataset
0.425
DGN
FM-GNN
0.994
0.400
0.933
0.375
Campe
0.901
DSN
R 0.8152
R 0.8912
LFMI-GNN
LOGN
10
70
80
56
10
70
30
60
denoising steps
. denolsing steps
Figure 19: flow-matching-basee models outperform diffusion- models when only small number of denoising steps employed for sampling. This was tested (a) the WING-TEST dataset and (b) the ALLIPSE-INDISS dataset.
31


Reviews:
[{"review_id": "IJUmRsYnfC", "reviewer": "Reviewer_5xrZ", "summary": "For many fluid dynamics problems, the exact trajectory dynamics are not relevant; instead, domain scientists are more concerned with quantities of interest that arise from steady state distributions. The paper proposes a multi-scale GNN architecture contextualized within the DDPM framework. Similar to prior and current state-of-the-art diffusion models, the work uses a latent representation, instead of performing diffusion over the original data space. In order to differentiate between both contributions, experiments are performed on both DGNs (Diffusion Graph Networks) and LDGNs (Latent Diffusion Graph Networks).\n\nExperiments are performed on 3 experimental settings: 1) pressure on the wall of an ellipse in 2D laminar flow 2) velocity and pressure around the ellipse 3) pressure on a 3D wing. The baselines include a Bayesian GNN, a Gaussian mixture model, a Gaussian regression model, and a VAE.", "questions": "- How does the proposed method differentiate itself from [1] (cited in the related works) which also uses a UNet style GNN? Similarly, how is the proposed GNN different / improved upon from [2]?\n- Are performance improvements over the baselines a result of the diffusion model formulation or the GNN design? Have the authors considered using a latent DiT or different sampling approach (e.g., DDIM)?\n\n[1] Haomin Wen, Youfang Lin, Yutong Xia, Huaiyu Wan, Qingsong Wen, Roger Zimmermann, & Yuxuan Liang. (2024). DiffSTG: Probabilistic Spatio-Temporal Graph Forecasting with Denoising Diffusion Models.\n\n[2] Lino, M., Fotiadis, S., Bharath, A., & Cantwell, C. (2022). Multi-scale rotation-equivariant graph neural networks for unsteady Eulerian fluid dynamics. Physics of Fluids, 34(8), 087110.", "limitations": "", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 2, "contribution": 3, "strengths": "- The studied problem is of particular interest to the broader physics-inspired ML community and domain scientists. \n- Both DGN and LDGN clearly outperform the tested baselines, and show promising results.\n- The mathematical and architectural details are well-written.\n- To my knowledge, this is one of the first works using diffusion models + GNNs in the space of ML4PDEs.", "weaknesses": "- The compared baselines are not widely used or state-of-the-art methods. At the same time, the three datasets, while representative or real-world tasks, are newly generated. This makes the results within the paper hard to contextualize within the broader literature.\n- The related works section on DDPMs for graph structure needs to further differentiate the proposed method from existing works.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Yxbuqtx4L5", "reviewer": "Reviewer_praC", "summary": "This paper presents two graph-based diffusion models for generating equilibrium states of fluid dynamics systems: one employing diffusion in physical space and another in the latent space generated by graph encoders. The presentation is clear and the methodological choices appear well-reasoned, particularly the use of multi-scale GNN as an encoder to support the latent diffusion model. The approach shows promise for generalization to broader applications where steady state prediction is needed.\n\nWhile technically sound, the primary limitation lies in the methodology's originality. The work largely assembles existing methods rather than introducing novel techniques, though this integration is done throughoutly.\n\nI assumed two important references should also be included and addressed to help with the literature review of this paper.\n\n1. \"Direct Prediction of Steady-State Flow Fields in Meshed Domain with Graph Networks\" should be included as a baseline:\n    - it directly addresses steady-state flow prediction using GNNs. This leads me to question when we should prefer deterministic versus generative approaches.\n    - I can see their generative method working well when \"steady states\" show variations (like the periodic wakes in von Kármán vortices), but for cases with strictly stationary steady states - like this paper's ellipse and airfoil examples - wouldn't a deterministic approach make more sense?\n    - I'd strongly encourage the authors to include this baseline and run comparative experiments and discuss these trade-offs based on results\n\n2. \"Efficient Learning of Mesh-Based Physical Simulation with Bi-Stride Multi-Scale Graph Neural Network,\" which shares similar ideas about multi-scale GNN architecture and more importantly an automatic graph hierarchy construction.\n    - Frankly, compared to a refered paper in the current manuscript: (Fortunato 2022) which requires manual drawing of coarsen layer, this paper's automatic graph hierarchy construction is closer pipeline-wise. Hence I think it should be acknowledged in the related work, and the related method paragraphs.\n    - While I don't think they need to conduct ablation studies here, as it's not the main focus, mentioning it would strengthen their literature review.\n\nI recommend a borderline accept for this paper now. I would consider improving my score if the authors address these missing references, and if they add the suggested baseline and the comparative analysis between deterministic and generative approaches, and thoroughly discuss when each method works best.", "questions": "- NA", "limitations": "", "rating": 8, "confidence": 5, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- Clear writing\n- The components of the whole pipeline proposed are rather reasonable\n- The experiment details, the network input/output structure, and the architectures are clear and hence the reproductbility should be easy", "weaknesses": "- A bit lack of innovation of the components of the whole method", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "VdEMUpA8FL", "reviewer": "Reviewer_j9He", "summary": "This paper uses Denoising Diffusion to sample states from the equilibrium distribution of fluid flow tasks. Instead of operating on a 2D grid, it employs a graph-based approach combined with a multi-scale Graph Neural Network (GNN) to process information. Using only short target trajectories, the model can learn the full distribution of solutions. Additionally, by combining the diffusion approach with a variational graph auto-encoder, a latent version that operates on a reduced latent graph with fewer nodes can, in certain scenarios, improve results.", "questions": "- Do you plan to release the code in the future, or are there specific reasons it hasn’t been shared yet? If not, could you provide further details on any proprietary components or dependencies that may impact reproducibility? In the paper, the authors only mentioned the release of the  datasets.", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 4, "strengths": "- Novelty and Innovation: The paper presents a novel and compelling approach by combining multi-scale Graph Neural Networks (GNNs), Denoising Diffusion Probabilistic Models (DDPM), and Variational Graph Auto-Encoders (VGAEs), pushing the boundaries in fluid flow modeling.\n\n- Clear and Well-Structured Presentation: The authors provide a thorough and well-organized presentation, clearly explaining the approach and the motivation behind it, making it evident why this method is both practical and impactful.\n\n- Relevant and Thoughtful Application: The paper demonstrates a meaningful application of the algorithm, showcasing how the method can effectively address specific challenges in fluid flow tasks.\n\n- Strong Results and Evaluation: The results are impressive, backed by a robust evaluation that validates the approach's effectiveness and potential advantages over existing methods.", "weaknesses": "- Lack of Code Submission: The absence of a code submission limits the reproducibility and accessibility of the results for further research and practical implementation.\n\n- Minor Formatting Issues: Acronyms are not consistently introduced upon first appearance, which could impact readability. For instance, \"DDPM\" is used without prior definition (e.g., on line 76).\n\n- Missing Citation: The original GAN paper by Goodfellow (2014) is not cited (e.g., on line 74).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "VENomttbkG", "reviewer": "Reviewer_84S9", "summary": "The paper aims to predict possible steady states of complex unsteady dynamical systems given their mesh discretization and physical parameters. With a novel graph-based (latent) diffusion model to predict equilibrium states of dynamical systems, the proposed model can handle unstructured mesh.", "questions": "- What is a statistical equilibrium?\n- What is the generality of the model? Does it generalize to different trajectories, different physical parameters, or even different domain shapes or sizes?\n- Generally, how complex a dynamical system is the model able to handle? Some dynamical systems such as the Lorenz system are chaotic, so are they still predictable with the proposed model?\n\n---\n**Update**\n\nThank you for your detailed response. I think I have a better understanding of the concepts, and I believe it is meaningful work. Thus, I decided to raise my score.", "limitations": "", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "- Integrating GNN with a diffusion model is an interesting idea, as it allows for better generality and great flexibility in unstructured mesh.\n- The paper claims to predict a full distribution of possible equilibrium states, although trained on states of a short trajectory.", "weaknesses": "Some critical information is missing. \n- The underlying physics problem needs to be explained more. I would expect some rigorous mathematical definitions such as PDEs to formulate the problem.  \n- The related work should cover more classical approaches to the problem of predicting equilibrium states of dynamical systems. \n- The learning problem should also be better formulated.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "ulaaQaFUwo", "reviewer": "Reviewer_uXMJ", "summary": "This paper proposes using graph-based diffusion methods to sample steady states in fluid dynamics problems.  Impressively, the training does not need simulations run all the way to convergence.  Behind the scenes, the proposed model uses a multiscale graph network, downsampled using strategies that are typical in numerical PDE.", "questions": "Line 107, “defined their discretization mesh” (typo?)\n\nCan the authors provide intuition on why the problem described on line 142 is well posed?  In particular, if the learning method exclusively sees short trajectories that end before reaching equilibrium, why should this information be sufficient to train a method that directly predicts converged states it hasn’t seen in the training data?  The results seem to indicate this is possible, but it seems counterintuitive to me that this task is reasonable.\n\nLine 169, “interpolates between the bounds of the process’ entropy” --- I wasn’t sure what this means.  Also, (Ho et al., 2020) isn’t written for graph-based models --- does this make a difference?\n\nThe receptive field of a “vanilla” DGN (line 180) is relatively small --- typically a vertex only “talks to” its k-ring, where k is the number of layers.  Does this affect performance of the proposed model, even after incorporating the latent diffusion method?\n\nLine 207:  What does “perceptually equivariant” mean?  How do you verify?\n\nLine 249:  no need for a hyphen in “bottleneck”\n\nLine 384:  Do statistical limitations get in the way of estimating W_2^graph accurately?  It seems like this would be a high-dimensional optimal transport problem, which itself can be difficult to estimate reliably from samples.  How many samples are used here?\n\nLine 451:  If I understand properly, the “3x faster” claim here would be substantially dampened if you had to include training time of the diffusion model.  Is this right?", "limitations": "", "rating": 8, "confidence": 2, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "The results in this paper are impressive, and the algorithm makes sense.  I am particularly glad to see relatively realistic test cases as well as reasonable engineering justification for the design behind the proposed architecture.\n\nAlso, thank you for preparing a polished, easy-to-read paper submission.  Relative to the other ICLR papers I had to read, this one had a thoughtful and clear discussion.", "weaknesses": "From an algorithmic perspective, the building blocks of the proposed method are unsurprising and reflective of larger trends in the fields (graph networks with a multiscale aspect, diffusion sampling, and so on).  So, I would not consider this work to be too exciting from a methodological perspective, although perhaps a reader with more experience specifically in numerical fluid dynamics would better be able to evaluate this aspect.  I am willing to overlook this, however, as the method seems solidly engineered (why invent a new approach to ML just to claim novelty?) and the results look great.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Oral)