TurtLES: A Large-Scale Benchmark for Turbulent 3D Neural PDE Surrogates
Authors: Jarno Platenburg, Armand Kassaï Koupaï, Paola Cinnella, Patrick Gallinari
Organizations: Sorbonne Université, CNRS, ISIR, Paris, France · Delft University of Technology, Faculty of Aerospace Engineering, Delft, The Netherlands · Sorbonne Université, CNRS, Institut Jean Le Rond d’Alembert, Paris, France · Criteo AI Lab, Paris, France
Machine-learning surrogates are increasingly used to accelerate computational fluid dynamics, yet progress is limited by the lack of benchmarks capturing realistic, time-dependent turbulent flows. We introduce TurtLES, a 14 TB dataset of 3D turbulent cylinder wakes in the shear-layer transition regime, computed via high-fidelity implicit large-eddy simulations. The dataset contains 400 trajectories of 400 time steps each, with 3-11 million points per frame, spanning variations in geometry, Reynolds number, and angle of attack. Unlike existing datasets, TurtLES combines (i) physically realistic turbulence with an active energy cascade driven by vortex stretching, (ii) dense irregular meshes, and (iii) long-horizon temporal dynamics. We benchmark state-of-the-art neural operators on autoregressive rollouts and identify key failure modes, including temporal instability, loss of small-scale energy, and challenges in learning on large irregular domains. TurtLES provides a new benchmark for evaluating surrogate models on 3D turbulent wakes at industrially relevant mesh resolution.
Figures & tables
Figure 1: TurtLES visualized through Q-criterion isosurfaces and its dataset parameters.
Figure 2: Overview of the parameter sampling space. (a) Definition of the super-elliptical geometry. (b) The Ns=400 samples visualized in three dimensions, including projections on the x,y,z planes in gray for reference. (c-f) Histograms of the sampled parameters. Note that the colors used for different values of n in (f) correspond with those in (a) and (b).
Figure 3: Schematic overview of the domain. (a) Three-dimensional view illustrating the physical boundaries, including the inlet ( A ), outlets ( B – D ), object surface ( E ), and periodic side-walls ( F , G ). (b) xy -plane of the computational mesh, which is extended into three dimensions via a structured extrusion in the span-wise ( z ) direction.
Figure 4: Overview of the flow case data structure. For each case, a single set of coordinates is stored with the solutions at 400 time-steps (e.g. t0,t1,...,tf , where tf=399 ). A solution comprises the three-dimensional pressure ( p ), density ( ρ ), and velocity ( u , v , w ) fields.
Field error
Spectral error
Model
Rel- L2↓
VRMSE ↓
Low ↓
Mid ↓
High ↓
UPT
7.01e-2
8.94e-2
4.69e-2
2.33e-2
4.42e-2
GINO
1.32e-1
1.21
1.32e-1
7.31e-2
4.48e-2
Transolver++
1.08e-1
1.06
8.02e-2
2.89e-2
4.41e-1
CALM-PDE
6.80e-2
7.26e-2
4.66e-2
2.38e-2
4.44e-2
Table 1: Next frame results on the test set. Lower is better for all metrics. Spectral errors are RMS log-power discrepancies computed over low-, mid-, and high-frequency spanwise wavenumber bands.
Rollout VRMSE ↓
Model
1 – 50
51 – 150
151 – 300
301 – 399
Cost [GPU-s] ↓
Solver (PyFR)
–
–
–
–
13160
UPT
6.51e-1
1.56
3.31
5.12
835
GINO
3.18
10.32
24.31
33.12
4097
Transolver++
2.33
8.32
19.33
31.24
1148
CALM-PDE
4.61e-1
5.46e-1
6.52e-1
7.89e-1
1073
Table 2: Rollout stability and computational cost. Rollout columns: VRMSE averaged over windows of rollout steps. Cost: GPU time to generate (solver) or predict (surrogates) one full trajectory, i.e. wall-clock time × GPUs used (four for the solver, one for each surrogate) on NVIDIA H100 80 GB.
Figure 5: Spectral fidelity of autoregressive rollouts. Panel (a) shows reference and predicted spanwise energy spectra averaged over test trajectories. Panel (b) shows the log-power discrepancy across wavenumbers.
Appendix figures & tables19 assets
Supplementary material from the paper’s appendix.
Appendix
Figure 6: Schematic overview of the computational domain. (a) Three-dimensional view illustrating the physical boundaries, including the inlet ( A ), outlets ( B – D ), object surface ( E ), and periodic side-walls ( F , G ). (b) 2D side view showing characteristic dimensions relative to the object thickness D . Note that black lines have no physical implications but indicate regions used for mesh refinement.
Category
Symbol (in Figure 6 a )
Boundary condition
Inlet
A
Riemann invariant
Outlet
B – D
Riemann invariant
Wall
E
Adiabatic no-slip wall
Left side
F
Periodic
Right side
G
Periodic
Appendix
Table 3: Boundary conditions of the simulation. Surface symbols refer to boundaries indicated in Figure 6 a
Figure 7: Two-dimensional ( xy -)plane of the computational mesh. The full three-dimensional domain is a 30-layer extrusion in the z direction. (a) Full domain. (b) Zoomed-in view of the geometric object.
Figure 8: Comparison between present simulations (solid black markers) and validation data from the literature (open and hatched markers) for (a) Strouhal number, (b) drag coefficient, and (c) recirculation length. Validation data are taken at discrete Reynolds numbers from the following sources: Re=1000 ( Jiang and Cheng, 2017 ) ; Re=1783 ( Bai and Alam, 2018 ) ; Re=2000 (hatched) ( Sharma and Dutta, 2022 ) and (open) ( Brun et al., 2008 ) ; Re=3000 ( Cimarelli et al., 2018 ) ; Re=3900 ( Lehmkuhl et al., 2013 ) .
Figure 9: Results for the coarse, medium and fine computational grids: (a) pressure coefficient distribution over the cylinder surface as a function of the angular position θ , defined from the front stagnation point ( θ=0∘ ) and increasing along the cylinder surface in the clockwise direction; and (b) mean streamwise velocity along the centerline ( y/D=0 ) downstream of the cylinder. Validation data are taken from Lehmkuhl et al. (2013) .
Figure 10: Results for the coarse, medium, and fine computational grids: mean wake velocity profiles of (a) streamwise component u and (b) transverse component v at downstream locations x/D=1.06 , 1.54 , and 2.02 . For clarity, profiles at x/D=1.54 and 2.02 are shifted downward. Validation data are taken from Lehmkuhl et al. (2013) .
Figure 11: Results for the coarse, medium, and fine computational grids: mean fluctuation profiles of (a) streamwise component u′u′ and (b) transverse component v′v′ at downstream locations x/D=1.06 , 1.54 , and 2.02 . For clarity, profiles at x/D=1.54 and 2.02 are shifted downward. Validation data are taken from Lehmkuhl et al. (2013) .
Figure 12: Comparison of velocity power spectra for (a) stream-wise and (b) cross-stream components at position (x/D,y/D)=(0.71,0.66) . Spectra are shown for both the full temporal resolution and the sub-sampled resolution Δts . The characteristic frequencies Stvs , Stsh , and StKH are annotated for reference alongside the theoretical −5/3 turbulence decay slope.
Figure 13: Comparison of velocity power spectra for (a) stream-wise and (b) cross-stream components at position (x/D,y/D)=(3.0,0.0) . Spectra are shown for both the full temporal resolution and the sub-sampled resolution Δts . The theoretical −5/3 turbulence decay slope is indicated for reference.
Figure 14: Comparison of (a) mean pressure coefficient distribution, (b) root-mean-square pressure coefficient fluctuations, and (c) skin-friction coefficient over the rectangle surface. Validation data are taken from Cimarelli et al. (2018) .
Figure 15: Comparison between the present study (top row) and Cimarelli et al. (2018) (bottom row) of the mean flow field above the rectangular surface. Streamlines of the mean velocity field are shown in (a), and the mean pressure field in (b). In (a), green lines indicate the primary vortex , red lines the secondary vortex , and black lines the wake vortex .
Figure 16: Mean velocity profiles in (a) from near to far downstream of the object, and Reynolds stress components: (b) u′v′ close to the cylinder and (c) u′u′ in the far wake region. Validation data are taken from Sharma and Dutta (2022) .
Figure 17: Numerical convergence criterion χ for (a) streamwise and (b) vertical velocity at probes P4 , P6.5 , and P9.5 for the [Re,α,a,n]=[1796,11,0.424,1] case. For v , χ is computed using standard deviation only. Stabilization indicates the evacuation of numerical transients and the onset of statistically steady flow.
Size
Shape family (% of split)
Split
Trajectories
Fraction (%)
n=1
n=2
n=5
n=∞
Train
320
80.0
25.0
25.0
25.0
25.0
Validation
40
10.0
25.0
25.0
25.0
25.0
Test
40
10.0
25.0
25.0
25.0
25.0
All
400
100.0
25.0
25.0
25.0
25.0
Appendix
Table 4: Size and geometric composition of the TurtLES train/validation/test split. Each subset reproduces the composition of the full dataset.
Re [-]
α [deg]
a [-]
Split
mean ± std
range
mean ± std
range
mean ± std
range
Train
2456 ± 839
[1001, 3889]
22.6 ± 13.2
[0, 45]
0.603 ± 0.365
[0.167, 1.499]
Validation
2434 ± 847
[1077, 3868]
22.5 ± 13.3
[2, 45]
0.614 ± 0.393
[0.172, 1.476]
Test
2431 ± 859
[1089, 3832]
22.4 ± 13.3
[0, 45]
0.604 ± 0.368
[0.178, 1.399]
All
2451 ± 839
[1001, 3889]
22.5 ± 13.2
[0, 45]
0.604 ± 0.367
[0.167, 1.499]
Appendix
Table 5: Parameter coverage of each split. Entries report the mean ± standard deviation and the [min,max] range of the Reynolds number Re , the angle of attack α , and the super-ellipse semi-axis parameter a .
Figure 18: Coverage of the parameter space by the train/validation/test split. The first three panels show the sampled cases projected onto the (Re,α) , (Re,a) and (α,a) planes, with training cases in gray, validation cases as blue squares and test cases as red triangles. The right-most panel shows the fraction of each subset belonging to each shape family n .
Geometry tier
Families n
T100
T200
T300
T400
G4
{1,2,5,∞}
3.5 TB
7.0 TB
10.5 TB
14 TB
G3
{1,2,∞}
2.6 TB
5.1 TB
7.7 TB
10.2 TB
G2
{2,∞}
1.8 TB
3.6 TB
5.3 TB
7.1 TB
G1
{2}
0.8 TB
1.6 TB
2.4 TB
3.3 TB
Appendix
Table 6: Approximate download size of each mini-TurtLES subset. Geometry tiers reduce the number of shape families; temporal tiers truncate each trajectory.
Figure 19: Visualization of eight TurtLES cases. For each case, 3D, side and top views show isosurfaces of the Q-criterion ( Hunt et al., 1988 ) ( Q∗=QD2/U∞2=1 ) at an instantaneous snapshot, colored by the non-dimensional velocity magnitude ∥U∥/U∞ .
Hyperparameter
UPT
GINO
Transolver++
CALM-PDE
Trainable parameters
61.96 M
3.10 M
21.79 M
45.14 M
Hidden dim
384
64 (FNO)
512
384
Heads
6
–
8
6
MLP ratio
4
–
4
4
Depth
2/4/2 (enc/lat/dec)
4 (FNO)
8
8
Tokens / latent grid
K=16,384 , L=4,096
123
128 slices
163=4,096
Appendix
Table 7: Summary of baseline hyperparameters. Counts are per model.
Evaluating neural operators for 3D turbulent flow requires validated datasets with physical benchmarks. We present a reproducible pipeline generating training data for 3D channel flows around generated geometries at Re=1,000-10,000. Our lattice Boltzmann solver with cumulant collision operators is rigorously verified against experimental measurements (Strouhal number, drag coefficients, turbulent fluctuations) with comprehensive grid convergence studies at resolution 1024x512x512. Building upon an established framework, this validated pipeline enables standardized surrogate model comparison. We outline planned systematic evaluation of Fourier Neural Operator and U-Net variants on forecasting, super-resolution, and error correction tasks, using physics-informed metrics to assess turbulent energy cascade representation. Future work will compare computational efficiency between numerical solvers and neural surrogates, exploring practical application. We seek community feedback on our validation approach, planned benchmark methodology, and evaluation priorities for neural operators in turbulent flows.
Lukas Schröder, Shubham Kavane, Harald Köstler
Chair of Computer Science 10 (System Simulation) Friedrich-Alexander-Universität Erlangen-Nürnberg
Turbulence is a central testbed for machine learning on physical dynamics because its governing laws are known exactly. However, most existing studies remain in 2D, while 3D turbulence has fundamentally different physics and is far more costly to simulate. Existing 3D resources also typically provide only one realization per configuration, making it difficult to distinguish learning the dynamics from fitting the statistics of a single flow. In this paper, we introduce TIDE (Turbulent Incompressible DNS Ensembles), a 256^3 DNS corpus and benchmark for 3D incompressible turbulence, with 15 configurations on eight controlled axes, independent ensembles, pressure fields, and equation-level verification. The benchmark includes five tasks, standardized learned baselines, controlled generalization splits, and physical-fidelity metrics alongside pointwise error. Across the main forecasting configurations, current learned models barely outperform persistence and still make about twice the error of a spectral solver given the true equations. Moreover, lower pointwise error can coincide with severely distorted small-scale dynamics, showing that accuracy alone does not ensure physical fidelity. Generalization results further show that most regime shifts reflect limited training coverage, whereas forced-to-decay transfer exposes a missing conditioning variable: operators trained under forcing continue to predict driven evolution when the external drive is removed. Closing these accuracy, fidelity, and conditioning gaps is the central open problem made measurable by TIDE.
Yilong Dai, Yiming Sun, Yiheng Chen +4
University of Alabama Tuscaloosa, USA · University of Pittsburgh Pittsburgh, USA
Neural surrogates are trained to predict 3D turbulent flows in place of direct numerical simulation (DNS). For chaotic flows, the goal is short-term pointwise accuracy followed by long-term physical and statistical fidelity. However, small prediction errors can carry a surrogate away from the flow's attractor. Off-attractor states are poorly represented in training data, leaving their evolution weakly constrained. The learned dynamics can then amplify deviations and lead to blow-up, freezing, or statistical drift. We propose off-attractor supervision (OAS) to supervise neural surrogates where they leave the attractor. OAS teaches the model how the true Navier-Stokes dynamics would evolve from these states. Each selected state is paired with its own future computed by DNS. Three generators select a few hundred states for relabeling. The first collects states from the surrogate's own rollouts. The second uses surrogate attacks to target freezing, excessive amplification, and violations of incompressibility and energy balance. The third perturbs training states along an amplified direction and a strongly damped random direction of the dynamics. All attacks run on the surrogate alone, and DNS relabeling is performed offline once per selected state. Experiments on 1283 turbulence show that OAS increases the median time to failure from 21 to 721 steps. The compared baselines achieve medians of at most 110 steps, and the advantage holds across training seeds. OAS also achieves the lowest pointwise error at step 15 and the best long-horizon statistics among the compared methods. OAS integrates physical models into neural simulation by extending supervision from fixed reference trajectories to states where the surrogate is likely to fail. This principle can guide the development of more reliable scientific surrogates when deployment takes models beyond the coverage of their training data.
Yilong Dai, Shaswata Mitra, Raj Patel +7
The University of Alabama · Rutgers University · University of Pittsburgh