ONE CYLinder: A Benchmark for Graph-Based Surrogate Modeling of Unsteady Bluff-Body Flows

Introduction

• Accurate long unsteady trajectories on unstructured meshes remain challenging because autoregressive errors accumulate, especially as wake dynamics become more complex.Source: §1 · p. 2
• Existing benchmarks often cover limited Reynolds-number ranges, simplified geometries, or specific applications, and rely mainly on pointwise reconstruction errors.Those evaluations provide limited insight into physical fidelity and long-term autoregressive stability.Source: §1 · p. 3
• 450 high-fidelity simulations across laminar, transitional, and high-Reynolds-number regimes form ONECYL, with time-resolved fields, geometry and mesh data, and integrated drag and lift quantities.Each regime contains 150 randomized cylinder configurations.Source: §1 · p. 3
• The benchmark combines full-field rollout errors, virtual probes, and integrated force predictions to evaluate numerical accuracy and physical consistency.Source: §1 · p. 3
• The physics-constrained Graph Transformer predicts velocity and pressure autoregressively on the mesh and supports investigation of geometry, input features, and physics-based losses across regimes.Its architecture combines sparse graph attention with divergence, convection, and velocity-gradient constraints.Source: §1 · p. 3

2. Solver methodology

• VMS finite-element simulations cover Reynolds numbers from 102 to 10^4 while varying cylinder diameter and position, storing velocity, pressure, drag, and lift over time.The simulations use the CIMLIB_CFD solver and provide the reference solutions for ONECYL.Source: §2 · pp. 3–4

2.1. Governing equations

• The incompressible Navier–Stokes equations govern the Newtonian fluid flow, and their weak formulation underlies the stabilized finite-element discretization.The equations express mass and momentum conservation.Source: §2.1 · p. 4

2.2. Variational Multiscale formulation

• The VMS formulation separates resolved and unresolved scales and adds stabilization terms derived from the governing equations while preserving Galerkin consistency.For linear P1 elements, second-order subscale derivatives vanish within each element; time integration treats convection semi-implicitly and viscous, pressure, and incompressibility terms implicitly.Source: §2.2 · pp. 4–5
• At high Reynolds numbers, VMS acts as an implicit LES closure, representing unresolved turbulent scales through its multiscale decomposition and stabilization without an explicit subgrid-scale model.Source: §2.2 · p. 5

2.3. Boundary conditions

• The canonical setup prescribes uniform inlet flow, zero-gradient outlet conditions, slip at the upper and lower boundaries, and no-slip at the cylinder.The setup is intended to make wake dynamics depend on Reynolds number and cylinder geometry.Source: §2.3 · p. 5

2.4. Meshing

• Gmsh generates body-fitted triangular meshes with boundary-layer refinement near the cylinder and progressive coarsening away from it; meshes contain approximately 1750±170 nodes.A nearly constant mesh complexity helps keep simulation differences tied primarily to flow physics rather than spatial-resolution variation.Source: §2.4 · pp. 5–6, Figure 1 · p. 6

3.1. Geometry, parameters, and data generation

• The two-dimensional incompressible-flow domain uses a cylinder with varied position and diameter, uniform inflow, no-slip cylinder boundary, and a traction-free downstream boundary.The rectangular domain measures 50 × 15 dimensionless units.Source: §3.1 · p. 7

3.1.1. Simulation parameters

• Simulations run for 120 seconds at 0.01-second steps, retaining every twentieth frame for 601 stored snapshots per case.The free-stream velocity is fixed at 1 m/s, with three kinematic viscosities generating the Reynolds-number families.Source: §3.1.1 · pp. 7–8

3. The ONECYL benchmark

• 450 simulations across three Reynolds-number regimes yield approximately 270000 flow snapshots spanning laminar, transitional, and turbulent wake dynamics.The benchmark reuses 150 randomized cylinder geometries at Re = 10^2, 10^3, and 10^4.Source: §3.1.2 · p. 8

3.1.2. Case variability and distributions

• Sampled geometries uniformly cover admissible region B without noticeable spatial clustering or geometric bias; diameter sampling is uniform across the domain, and the same set is reused across Reynolds-number regimes.Across Re2, Re3, and Re4, effective Reynolds numbers are 142, 1420, and 14200, respectively, with wakes progressing from laminar vortex shedding to increasingly complex turbulent dynamics.Source: §3.1.2 · pp. 8–9

3.2. Features

• ONECYL provides time-dependent velocity, pressure, drag, and lift alongside node type, Reynolds number, and a level-set descriptor of cylinder distance.The level-set gives each mesh node a continuous geometric representation beyond discrete node-type labels.Source: §3.2 · pp. 9, 11
• Drag and lift histories show progressively richer temporal dynamics across the three regimes, from regular periodic shedding to more irregular wake behavior.Source: §3.2 · p. 11

3.3. Computational infrastructure

• The VMS-stabilized finite-element simulations use CIMLIB_CFD and run on a 32-node high-performance computing cluster.Each node has two AMD EPYC 7502 processors and 256 GB of RAM.Source: §3 · p. 6, §3.3 · p. 11

4.1. Graph representation

• The mesh is represented as a graph, and each node becomes a Transformer token; the proposed model omits edge attributes and uses adjacency to encode connectivity.Node features include flow variables, coordinates, distance-to-object, Reynolds number, and node type.Source: §4.1 · pp. 12–13, Encoder and decoder · p. 14

4.2. Transformer-based model

• The model replaces the conventional message-passing processor with Transformer blocks, making one-step predictions that are recursively rolled out into flow trajectories.The architecture incorporates geometry through node coordinates and connectivity, without explicit edge attributes during feature propagation.Source: §4.2 · p. 13

Encoder and decoder

• A two-layer MLP encoder projects node features into latent embeddings, and a two-layer MLP decoder maps the final embeddings back to predicted physical quantities.Both networks use RMSNorm followed by ReLU activations.Source: §4.2 · p. 13, Encoder and decoder · p. 14

4.3. Transformer processor

• The processor stacks Transformer blocks that alternate masked multi-head self-attention with a gated feed-forward network, using residual connections and RMS normalization.The attention layer uses adjacency as a sparse mask, restricting nodes to graph neighbors and excluding self-loops.Source: §4.3 · p. 14, Figure 7 · p. 14, Masked multi-head self-attention · p. 15

Masked multi-head self-attention

• Sparse attention restricts interactions to mesh neighbors, reducing computational and memory costs while preserving mesh locality.Random long-range connections can shorten graph distances between distant regions with limited added computational overhead.Source: Masked multi-head self-attention · p. 15

Gated MLP and normalization

• The gated feed-forward branch complements graph-constrained attention with channel-wise nonlinear transformations, while training minimizes one-step MSE and adds noise to dynamic features for rollout robustness.The final training objective supplements the data-driven loss with physics-based constraints.Source: Gated MLP and normalization · p. 16, Learning · p. 16, Noise · p. 16

4.4. Physics-constrained GNN

• The physics-constrained objective combines data loss with continuity, convection, and gradient penalties computed using differentiable operators adapted to irregular meshes.Continuity penalizes incompressibility violations, convection penalizes differences in advective dynamics, and gradient loss targets spatial-gradient discrepancies.Source: §4.4 · pp. 17–18

5.1. Model assessment metrics Autoregressive rollout predictions

• Evaluation combines global rollout errors at short, medium, and long horizons with nine fixed wake probes that assess local flow dynamics and time-series behavior.The probes sample the near wake, shear layers, and downstream recovery region; values are taken from the nearest mesh node rather than interpolated.Source: Rollout errors · p. 20, Sensor points · pp. 20–21, Sensor-based evaluation · pp. 27–28

Integrated surface forces

• Integrated drag and lift predictions complement field and probe errors by testing the physical consistency of the predicted flow through global surface forces.The forces are computed at every rollout step by integrating predicted pressure and velocity over the cylinder surface and compared with high-fidelity CFD reference values.Source: Integrated surface forces · p. 21

5.2. Initializer–iterator dual-model approach

• A two-model initializer–iterator strategy separates the initial-boundary-condition prediction from subsequent autoregressive flow evolution, improving long-horizon rollout accuracy and stability over a single trajectory-trained model.The initializer predicts the first physical state, and the iterator propagates the remaining rollout; this separation lets each model specialize in its task.Source: §5.2 · pp. 21–22, Figure 10 · p. 22

5.3. Initializer performance

• Adding level-set geometry reduces mean initializer error versus VPN-l by approximately 45% (Re2), 43% (Re3), and 20% (Re4), while VPLN-ld lowers test-case variability.The level-set field helps reconstruct strong spatial variations near the solid boundary, particularly pressure; representative errors are localized near the cylinder, where gradients are strongest.Source: §5.3 · pp. 24–25, Table 6 · p. 25

Rollout errors

• Divergence regularization degrades Re2 rollouts but benefits long-horizon Re3 and Re4; near-equivalent one-step errors conceal stability differences, making complete autoregressive rollouts essential to evaluation.The constraint is increasingly useful where stronger advection, vortex shedding, and flow complexity challenge prediction.Source: Rollout errors · pp. 25–26
• Across 600-step rollouts, VPLN-l lowers Re2 complete-rollout error from 73.40 × 10−3 to 51.87 × 10−3, while VPLN-ld reaches 132.63×10−3 (Re3) and 170.26×10−3 (Re4).VPLN-ld is the lowest-error variant in Re3 and Re4; in Re4, VPLN-l performs better over the first 50 steps, so divergence regularization primarily benefits longer-term behavior.Source: Rollout errors · pp. 25–26, Table 7 · p. 26

Flow-field visualizations

• The models remain numerically stable over 600 steps, but increasingly lose wake phase and amplitude accuracy as flow complexity rises, while preserving the principal wake topology.For the Re3 VPLN-ld case, predicted wake and pressure structures remain close to references at steps 50 and 600 without obvious non-physical structures.Source: Figure 13 · p. 27, Flow-field visualizations · p. 27, Figure 14 · p. 28

Sensor-based evaluation

• Sensor predictions capture transient development and regime-dependent oscillations, while long-horizon pointwise discrepancies mainly reflect phase and amplitude errors rather than instability.Oscillation frequency and amplitude increase with Reynolds number; level-set variants generally improve sensor agreement, especially in Re3 and Re4, while divergence regularization benefits some higher-Re cases but not consistently lower-Re predictions.Source: Sensor-based evaluation · pp. 28, 30

Predicted-force evaluation

• Predicted forces retain physically consistent mean levels and fluctuation magnitudes, with level-set variants generally more accurate than VPN-l, especially at higher Reynolds numbers.In Re4, lift is less accurate than drag because its unsteady asymmetric wake dependence makes it more sensitive to phase and amplitude errors; force trends persist across test geometries and Reynolds numbers.Source: Figure 19 · p. 30, Predicted-force evaluation · p. 31

5.5. Combined initializer–iterator performance

• 66.88 × 10−3 instead of 71.14 × 10−3 (Re2), 131.10 × 10−3 instead of 132.63 × 10−3 (Re3), and 168.50 × 10−3 instead of 170.26 × 10−3 (Re4) were mean full-rollout RMSEs with learned versus ground-truth initialization, not systematic improvements.The initializer increases short-horizon error, but the difference does not grow across the 600-step rollout; Gaussian noise during iterator training may help it recover from perturbed inputs.Source: Table 8 · p. 33, §5.5 · p. 33
• The combined pipeline produces complete trajectories from prescribed initial and boundary conditions, preserving principal flow structures, sensor dynamics, and integrated forces without CFD-computed initialization.Its learned initial state remains close to the ground-truth-initialized rollout, including force trends; minor phase and amplitude differences do not become unstable or divergent.Source: §5.5 · pp. 33–36, Figure 25 · p. 35

5.6. Limitations and future work

• Limitations include narrow geometry coverage, weak cross-regime generalization, imperfect transient timing and phase metrics, greater Re4 difficulty, and CFD resolution constraints at still higher Reynolds numbers.Future work calls for broader, more uniformly sampled geometries and physical parameters, dedicated phase-sensitive and spectral temporal metrics, and finer meshes with improved near-wall or adaptive refinement to capture increasingly complex flow scales.Source: §5.6 · pp. 36–37

6. Conclusion

• The 450-simulation, 270,000-snapshot ONECYL benchmark spans three flow regimes and supplies fields, geometric descriptors, Reynolds numbers, and aerodynamic quantities.Its data support evaluation from full-field reconstruction to integral quantities such as drag and lift.Source: §6 · p. 37
• The physics-constrained Graph Transformer predicts velocity and pressure autoregressively on unstructured meshes, with stable rollouts over 600 time steps.It combines sparse multi-head self-attention and physics-based regularization while accommodating arbitrary mesh topologies.Source: §6 · p. 37
• Level-set geometry encoding consistently improves accuracy and generalization to unseen cylinders, while divergence regularization helps increasingly as Reynolds number rises.Its benefits are strongest at medium and long horizons; divergence regularization offers limited improvement in the laminar regime but improves stability and reduces variability in transitional and turbulent flows.Source: §6 · p. 37
• Field, sensor, and force evaluations show dominant wake dynamics persist over long rollouts, while phase and amplitude errors have limited impact on reconstructed aerodynamic forces.The remaining temporal discrepancies occur in the statistically stationary oscillatory regime rather than through numerical instability or error accumulation.Source: §6 · p. 37
• The initializer–iterator pipeline starts from prescribed initial and boundary conditions without a CFD-computed initial flow field, and its iterator absorbs initialization perturbations during rollout.Long-term predictions remain comparable to those initialized with the exact first state.Source: §6 · pp. 37–38
• Higher-Reynolds-number extensions require finer or adaptive CFD meshes, architectures for multi-resolution discretizations, and metrics targeting phase accuracy, transient timing, and oscillatory dynamics.Broader geometric variability and explicit Reynolds-number conditioning are also identified as directions for unified multi-regime models.Source: §6 · p. 38

Appendix A. ONECYL dataset splits

• Case IDs retain the same cylinder geometry across the three Reynolds-number subsets, so each matched ID has different Reynolds-number flow data.The first digit identifies the single-cylinder dataset and the remaining four digits give the shared case number.Source: Appendix A · p. 48
• The full split uses 120 training, 15 validation, and 15 test cases; the mini split selects 50 training and 5 validation cases while retaining all 15 test cases.Both splits follow a 90/10 train/validation split.Source: Appendix A · p. 48

Appendix B. Force computation from flow fields

• Drag and lift are computed from pressure and viscous stresses along each cylinder boundary, using the discrete velocity and pressure fields on the triangular mesh.Boundary edges are identified using mesh connectivity, nodetype flags, and object bounds.Source: Appendix B · p. 50
• Edge forces are integrated around each object to produce force time series for comparison with ground truth or calculation of average force coefficients.The level-set gradient helps orient the outward normal, including for complex geometries.Source: Appendix B · pp. 50–51

Appendix C.1. Iterator model predictions: flow visualization

• Re4 predictions differ from ground truth in wake velocity amplitude and vortex-shedding period, although the principal wake structure remains identifiable.The visualizations compare predicted and ground-truth fields at steps 50 and 600, with predictions shown above ground truth.Source: Appendix C.1 · p. 51, Figure C.28 · p. 52
• The Re2 velocity panels show predicted and ground-truth fields at rollout steps 50 and 600 for case 10070 (Re = 142).The predicted fields appear in the top row and ground truth in the bottom row.Source: Figure C.26 · p. 51
• Pressure panels compare predictions with ground truth at steps 50 and 600 for Re2 (Re = 142) and Re4 (Re = 14, 200).The Re4 panels correspond to case 10070.Source: Figure C.27 · p. 52, Figure C.29 · p. 52

Appendix C.2. Iterator model predictions: sensor based evaluation

• Sensor plots compare velocity norm and pressure time series at nine points across Re2, Re3, and Re4 test cases.Ground truth is black; VPN-l, VPLN-l, and VPLN-ld predictions are blue, orange, and green, respectively.Source: Appendix C.2 · p. 53
• Compare feature-loss combinations across the three regimes using the same sensor-based velocity and pressure signals.Figures C.30–C.35 cover Re2, Re3, and Re4; zoomed views show six of the nine sensors.Source: Appendix C.2 · p. 53
• The Re2 figure shows velocity-norm time series at nine sensor points for different feature-loss combinations.The plotted quantity is ∥u(t)∥2 for case 10070 (Re = 142).Source: Figure C.30 · p. 53

• Mean and standard deviation of drag and lift are reported for iterator models across the Re2, Re3, and Re4 test sets.The statistics are computed from predicted flow fields.Source: Appendix C.3 · p. 53
• CFD ground-truth forces are gray, while VPN-l, VPLN-l, and VPLN-ld predictions are shown in blue, orange, and green, respectively.Source: Appendix C.3 · p. 53
• Force statistics use the last 400 rollout steps to avoid initial transient effects.Source: Appendix C.3 · p. 53

• Sensor time series compare feature-loss combinations for Re2 pressure at nine sensors and Re3/Re4 velocity and pressure at the first six of nine sensors.Source: Figure C.31 · p. 54, Figure C.32 · p. 54, Figure C.33 · p. 55
• Mean and fluctuating drag and lift derived from predicted fields are shown across the Re2, Re3, and Re4 test sets for different feature-loss combinations.Source: Figure C.36 · p. 56, Figure C.37 · p. 56, Figure C.38 · p. 57
• Re2 and Re4 field visualizations compare combined initializer-iterator predictions with truth-iterator predictions and CFD ground truth at rollout steps 50 and 600.The combined pipeline initializes the VPLN-ld iterator with a one-step prediction from the VPLN-l initializer.Source: Appendix C.4 · pp. 53, 57, Figure C.39 · p. 57
• Velocity-norm and pressure fields are shown for Re2 and Re4 test cases, with ground truth in the bottom row, truth-iterator predictions above, and combined predictions in the middle row.Source: Appendix C.4 · p. 57, Figure C.39 · p. 57

• Mean and standard deviation of drag and lift are computed from combined initializer-iterator predictions across the Re2 and Re4 test datasets.The iterator is initialized with a one-step prediction from the VPLN-l initializer.Source: Appendix C.5 · p. 57
• For Re2, derived forces from the VPLN-ld combined pipeline are compared with ground-truth and truth-initialized forces across the test set.Source: Figure C.43 · p. 59
• For Re4, derived forces from the VPLN-ld combined pipeline are compared with ground-truth and truth-initialized forces across the test set.Source: Figure C.44 · p. 59