Measured on EuroHPC MeluXina · A100-SXM4-40GB

Meshes that don’t fit
on any single GPU

The whole SIMPLE loop stays resident on the GPUs — mesh, fields, assembly, turbulence and the pressure multigrid. Nothing is copied back to the host between iterations.

Every figure below is traceable to the Slurm job that produced it.

147,449,586

cells solved across 16 GPUs

2.9e-10

agreement with single-GPU reference

29.9 s

to partition a 148 M-cell mesh

2.4%

mean drag error across four published cases

0.18%

drag variation between 8 and 16 GPUs

Capability: a 40 M-cell mesh no single A100 can hold

One A100-40GB cannot load this case at all — it fails allocating the field pool. Adding GPUs is what makes it solvable, and then what makes it faster.

1 GPU
out of memory — 40 GB exceeded
2 GPUs
3,278.8 ms
4 GPUs
2,145.3 ms
8 GPUs
2,497.1 ms
16 GPUs
1,564.3 ms
milliseconds per SIMPLE iteration — lower is betterjob 4983803
Decomposition scales sublinearly with mesh size

Distributed graph partitioning across 16 GPUs. An 74× increase in cells costs 16× the partitioning time. Partitions are cached per mesh and GPU count, so repeat runs skip this entirely.

2 M
1.8 s
20 M
6.3 s
40 M
10.9 s
80 M
17.0 s
148 M
29.9 s
seconds to partition, 16 GPUs — lower is betterjob 5031623
DrivAerML 147 M cells: where the first 150 seconds go

A real industrial automotive case — 147,449,586 cells, 59 boundary patches, meshed in ANSA. Startup is paid once; the run then proceeds at roughly 30 s per iteration.

mesh read
25.1 s
decompose
29.6 s
build AMG
11.6 s
build solver
76.7 s
first step
4.5 s
seconds, one-time startup on 16 GPUsjob 5031974

Correctness is gated before any timing

No performance number is reported from a build that has not first reproduced the single-GPU reference solution, field by field, on a converged case.

ConfigurationVelocity (L2 relative)Pressure (L2 relative)
2 GPUs vs 1 GPU2.912e-103.943e-09
4 GPUs vs 1 GPU3.955e-102.256e-09

Reproduced independently on two toolchains (CUDA 12.8 + NVSHMEM 3.4.5, and CUDA 13.3 + NVSHMEM 3.6.5).

Validated against a public benchmark

Every section above shows brae reproducing itself — multi-GPU output matching single-GPU output to a relative L2 difference of 3e-10. This section shows it reproducing published physics across four independent cases from the DrivAerML open dataset, an open reference for automotive external aerodynamics. Both reference drag coefficients are published by the dataset authors, so the comparison is fully auditable.

CaseCellsbrae CdPublished CdDifference
DrivAerML run_1147,449,5860.3075040.31092-1.1%
DrivAerML run_2137,102,1330.2765880.29529-6.3%
DrivAerML run_3143,870,6570.2487550.25239-1.4%
DrivAerML run_4137,424,2990.2723450.27416-0.7%
Mean absolute deviation2.4%

All four runs at 6000 SIMPLE iterations.

At 6000 iterations brae reproduces the published drag ordering exactly: run_1 > run_2 > run_4 > run_3. Every case is under-predicted — by −1.1%, −6.3%, −1.4% and −0.7% — a systematic bias of about −2.4%, not scatter. That is the expected direction for steady RANS compared with the time-averaged scale-resolving simulation the reference used: RANS suppresses the separated wake and consequently underestimates pressure drag. A systematic offset is more useful than scatter because it is understandable and can be accounted for when interpreting results.

How we know these are converged

Force coefficients are evaluated at every checkpoint, so each run produces a drag history rather than a single end-of-run number.

On run_2, carried to 15,000 iterations, Cd settles from 4000 onward to a mean of 0.27623 within a band of 0.27486–0.27769 — roughly ±0.5%. The 6000-iteration value used in the table above sits within 0.13% of that converged mean.

Residual convergence is not force convergence. On another case the velocity residual was marginally worse at 6000 iterations than at 3000, while the drag coefficient moved 9%. Judging convergence by residuals alone is what produced the earlier, wrong numbers.

run_1 lift and side force (steady RANS vs SA-DDES)

QuantitybraePublishedDifference
Lift Cl0.10920.0694+57%
Side Cs0.08380.0488+72%

The same answer on 8 GPUs and on 16

GPUsDrag CdLift Cl
80.3080460.117852
160.3075040.112735
difference0.18%4.5%

Changing the GPU count changes the decomposition entirely: different partitions, different halo topology, different treatment at every partition boundary. This is a real test of the distributed path, not a formality.

Drag differs by 0.18%. A separate 15,000-iteration study on this case shows the drag coefficient oscillating about its converged mean by roughly ±0.5% from one checkpoint to the next. That means 0.18% is inside normal run-to-run scatter — the two GPU counts are statistically indistinguishable, and presenting 0.18% as a precision figure would be misleading.

Lift differs by 4.5%, more than drag, because it is dominated by the rear wake — the least settled part of the flow and the part where steady RANS is most sensitive to boundary placement.

Drag agrees; lift and side force do not — and that split is expected. Cd sits within 0.4% of the reference. Cl and Cs are over-predicted by 57–72%. That is the expected signature of steady RANS: drag is dominated by front stagnation and overall body pressure, which RANS predicts well; lift is dominated by the rear wake and separation zones, which steady RANS predicts badly. A result where Cl also matched to within 1% would be more suspicious, not less.

These are different turbulence models. The reference values are time-averaged SA-DDES — transient and scale-resolving. brae ran steady SA RANS. This is a sanity check of the kind found routinely in the RANS-vs-DES literature, not a claim that brae reproduces DES.

Reference convention matters. The comparison uses the constant reference convention (Aref 2.17 m², lref 2.786 m), which is the pair brae’s case setup used. DrivAerML also publishes coefficients under a run-specific reference; comparing against that pair instead would shift Cd by 2.4%.

Run conditions. 147,449,586 cells · 16× NVIDIA A100 · 3000 SIMPLE iterations · final residuals Ux 5.6×10−4 and p 1.7×10−4 · 6861 s wall time.

How it scales from 8 GPUs to 64

One case (DrivAerML run_1, 147,449,586 cells), one turbulence model (Spalart-Allmaras), one iteration count (2000), one binary. The GPU count is the only variable. Startup — mesh read, decomposition, distribution — is measured separately from the solve and is not folded into the per-iteration figure.

GPUsNodesCells/GPUStartupms/iterSpeedupEfficiency
8218.4 M227.8 s3163.61.00×
1649.2 M201.5 s1800.21.76×88%
3284.6 M188.3 s1333.32.37×59%
64162.3 M260.1 s1120.02.83×35%

Time to a converged 6000-iteration solve

GPUsTime to 6000 iterations
85.3 h
163.0 h
322.2 h
641.9 h

Scaling is near-linear to 16 GPUs (88% efficiency), useful at 32 (59%), and clearly diminishing at 64 (35%). There is no single headline figure that honestly represents the whole range — the curve has a knee, and knowing where that knee is matters more than an average that conceals it.

The mechanism is visible in the cells/GPU column. At 64 GPUs each partition holds 2.3 M cells, and halo exchange between partitions starts to dominate the arithmetic work. This is the expected shape of a strong-scaling curve for a distributed solver; diminishing returns at high rank counts follow from the ratio of communication to arithmetic, not from a defect in the implementation.

Startup does not improve monotonically with more GPUs. It falls from 228 s at 8 GPUs to 188 s at 32, then rises to 260 s at 64 — because decomposition and inter-node distribution become more expensive as the rank count grows. Startup is reported separately precisely so this cost is not hidden inside the headline solve rate.

Concurrency caveat. All four runs were executed concurrently on the cluster and therefore shared the interconnect. The 64-GPU run spans 16 nodes and is the most exposed to that contention. Its 35% efficiency should be read as a lower bound — some of the loss may be network contention rather than a property of the solver.

The four decompositions agreed on the drag coefficient to within 0.11% (0.306222 to 0.306559). This is agreement at a fixed iteration count between the four runs. At 2000 iterations the case is not converged, so that number is not an accuracy claim and cannot be compared with published data.

Two turbulence models, on the same case

Both runs use DrivAerML run_1 (147,449,586 cells, 6000 iterations). The published reference for both is Cd 0.31092, Cl 0.0694.

ModelGPUsCdvs publishedClvs published
Spalart-Allmaras160.307504-1.1%0.112735+62%
k-epsilon320.295209-5.1%0.108892+57%

Both models under-predict drag. Spalart-Allmaras is closer to the reference (−1.1% against −5.1%), which is the expected result for external aerodynamics with largely attached flow: SA handles the attached boundary layer well, while k-epsilon over-predicts turbulent mixing, delays separation, and consequently loses more pressure drag.

Both models over-predict lift by a similar and large margin — 57–62%. That margin is a property of steady RANS, not of either model in particular: the lift on this body is dominated by the rear wake, which a steady formulation cannot resolve. Neither lift number should be trusted at that level of detail, regardless of which model produced it.

The k-epsilon run used 32 GPUs (8 nodes), exercising the solver at that scale as well. Neither model is presented as the right choice — they are different tools with different known biases, and that is precisely what the results confirm.

Speedup arrives with size

Extra GPUs buy nothing on a mesh that already fits one card, and a great deal on a mesh that does not. Both of these are measured on the same solver.

Mesh8 → 16 GPUsSpeedupEfficiency
8 M cellsfits one GPU1.04×essentially none
147 M cells10573 s → 6861 s1.54×77%

At 8 M cells the per-iteration cost is dominated by blocking global reductions that do not shrink with subdomain size, so extra GPUs buy nothing. At 147 M cells each GPU holds about 9.2 M cells — enough work to amortise those reductions. Both rows are the converged, drag-validated DrivAerML case. The honest summary: multi-GPU buys capability at every size, and buys speed once the mesh is large enough to need it.

What we do not claim

Benchmark pages that report only good news are not worth reading. These are the limits we have measured in our own system.

More GPUs does not mean proportionally faster

On a mesh that already fits one GPU, 16 GPUs return about 1.04× — essentially nothing. The per-iteration cost is dominated by blocking global reductions in the Krylov solver, and that cost does not shrink with subdomain size. Multi-GPU here buys capability, not speedup. That figure is measured at 8 M cells; the 147 M-cell case scales at 77% efficiency from 8 to 16 GPUs — see "Speedup arrives with size" above.

Krylov iteration counts grow with GPU count

The distributed pressure preconditioner is one-level block-Jacobi, so its conditioning degrades as subdomains multiply: 61 iterations per step on 8 M cells becomes 273 on 147 M across 16 GPUs. A two-level coarse space cuts the count but has so far lost on wall-clock time.

These are not validated aerodynamic predictions

The published DrivAerML setup is transient SA-DDES. We run steady SIMPLE. The result demonstrates that the machinery handles an industrial mesh at full size — it is not a claim about drag coefficients.