When a neural surrogate cannot accelerate a solver: runtime share, closed-loop drift, and the economics of uncertainty gating in a stiff coupled simulation
Authors: L. Thümmler, T. Kuroda
Organizations: Department of Mathematics, ETH Zürich, Rämistrasse 101, 8092 Zürich, Switzerland · Max Planck Institute for Gravitational Physics (Albert Einstein Institute), Am Mühlenberg 1, D-14476 Potsdam, Germany
Learned surrogates for expensive inner solver blocks are a widely pursued route to faster multiphysics simulation. We report a controlled, end-to-end negative result and identify three structural barriers, none of them a deficiency of the network we trained. The testbed is the implicit Newton solve coupling energy-dependent neutrino radiation to matter in a general-relativistic radiation-hydrodynamics code, its most expensive physics routine per call. First, per-call cost and share of runtime are different quantities, and only the second bounds acceleration. An exclusive self-time profile puts the target block at 16.9% of critical-rank wall clock, capping any surrogate at ~1.2x by Amdahl's law. A surrogate 5.8x cheaper per call merely ties the solver, and the configuration stable enough to run without fallback reaches only parity. Second, offline accuracy cannot rank surrogates for deployment: across fourteen networks the pooled Spearman error-versus-survival correlation (rho = +0.73) is a between-family confound that vanishes under control (rho = -0.04). Third, a correct out-of-distribution gate cannot accelerate a loop that leaves its training distribution. We give the break-even deferral fraction in closed form: because the visited states sit 73x off the data manifold, the gate defers 96.8 to 99.7% of cells, almost invariant to surrogate quality. Including its own cost, the gated loop is a 0.94 to 0.96x slowdown. We further separate stability from fidelity: a never-crashing gated run accumulates a linear -19.9% density bias over 6000 steps. The error is a directed, ballistically accumulating bias, not the variance-driven divergence the autoregressive literature targets.
Figures & tables
Routine (disjoint leaf)
Pre-bounce
Bounce
Post-bounce
Ghost exchange, mostly idle wait
30.8
27.8
34.6
PPM advection
20.3
27.6
21.7
Ghost exchange, genuine MPI
13.5
9.1
12.0
Loop overhead residual
11.0
10.1
11.8
Implicit coupling solve
7.5
16.1
10.5
Periodic I/O
9.8
2.5
3.1
Table 1: Exclusive self-time breakdown of one time step, as a percentage of mean wall clock across MPI ranks, for the three phases of the collapse. Leaves are disjoint; the six named categories are the ones large enough to discuss individually, and the residual row recovers the remainder so that each column sums to 100% . The lower block gives the solver’s share in two other framings and the mean number of Newton iterations per cell, which is the physical reason the solver cost varies across phases.
Figure 1: Per-step accuracy of the wide model on a held-out split of the training distribution. (a) Predicted against true electron fraction with the identity line ( R2=0.982 , mean absolute error 3.4×10−3 ). (b) Held-out R2 per predicted quantity. These in-distribution figures are optimistic, and the remainder of this section shows they carry no ranking signal for deployment.
Figure 2: Offline error does not rank surrogates by closed-loop survival. Fourteen independently trained networks, all scored on one common held-out set (verified excluded from every training set), against closed-loop steps survived. The points separate into two clusters, narrow specialists (cyan, upper right) and broad models (red, lower left), and the apparently strong pooled correlation in the left panel is produced by that separation: within either cluster there is no trend, and the partial rank correlation controlling for cluster membership is ρ=−0.04 ( p=0.89 ). Restricting the error to the dense core (centre), where the broad models actually fail, does not help, nor does the multi-step rollout error growth factor (right). Circles are the five seeds of each configuration; stars are earlier single-seed models.
Configuration
Steps survived
Density
Full identity (network inert)
300 (clean)
stable
Network drives Ye +QSN+entropy only
∼170
stable, then Ye drift
Network drives everything (including τ,m1 )
2
collapse 100×
Table 2: Identity ablations isolate two obstacles. The framework itself is sound ( 300 clean steps with an inert network). Predicting the conserved energy and momentum causes an acute density collapse at step 2 ; with those held fixed, a slow electron-fraction drift terminates the run near step 170 .
Figure 3: Signal-to-noise of the energy channel. The true per-step energy change is ∼5500× smaller than the network error: the increment is below the learnable resolution.
Figure 4: The closed-loop drift is covariate shift. Network-visited states (red) lie far off the real-data manifold (blue): 73× the real-to-real nearest-neighbour distance, 61% beyond the 99 th-percentile Mahalanobis shell, with unphysical entropies ( ≤62 versus a real maximum of 7.8 ) that no real trajectory occupies.
Figure 5: Closed-loop horizon (steps to a genuine crash) versus aggregation iteration. Naive aggregation oscillates without converging (best 286 , output model worse); keep-best plus a recency window reaches 557 steps, 5.5× the baseline ( 102 ), with Ye held physical.
Figure 6: Closed-loop horizon for five independently seeded trainings of the stacked configuration and of its reduced-width counterpart, against the single-run figures quoted elsewhere in this paper (stars). The original 2625 -step figure falls outside its own five-seed distribution; the reduced-width single run ( 2768 ) falls inside its own.
Configuration
Horizon (steps)
Speedup / defer
Note
Naive surrogate (predicts τ,m1 too)
∼ 2
n/a
acute sub-resolution energy collapse ( 100× density)
+ thermodynamic energy
102
pure replacement
acute obstacle removed; chronic Ye drift remains
Aggregation, naive iteration
102 – 286 (oscillates)
pure replacement
does not converge; last model = best
Aggregation, keep-best + recency
557 ( 5.5× )
pure replacement
Ye physical; delays but does not remove the crash
+ Jacobian contractivity
2625† ( 26× )
pure replacement
levers stack; contractivity alone gives ∼ 36%
Per-cell Mahalanobis gate
>1000
95 – 99% solver
stable, but ≈ 0.95 × once gate cost is included
Table 3: Closed-loop horizon by method: steps the live simulation survives before a genuine crash, with the resulting speedup or solver-deferral fraction. Horizons are not directly comparable across rows, as they probe different phases and stabilisations; all are in the early-infall regime. The un-aided baseline for the “ × ” factors is the thermodynamic-energy model without drift coverage ( 102 steps).
Figure 7: Detection quality of the Mahalanobis gate. (a) ROC curves for three definitions of a “bad” step: the gate separates large blow-up errors cleanly (err >3×10−2 : AUC 0.996 ) but the small accumulating errors that constitute the drift only weakly (err >10−3 : AUC 0.658 ). (b) Gate AUC as a function of the error threshold, rising from 0.66 to 0.996 : the detector is sharp exactly for the failures that are not the binding ones.
Figure 8: Break-even deferral fraction as a function of gate cost, from Eq. ( 2 ). Cutting the gate’s cost by 24× (deployed g=0.391 to idealised diagonal g=0.016 ) only moves the break-even from 0.53 to 0.98 : the formula’s dependence on g is far shallower than its dependence on d near d=1 . The measured deferral, 0.97 – 1.00 , sits at or beyond all three.
Figure 9: Central density of a gated run relative to a solver-only reference from the same starting state, pre-bounce. After a quiescent window of ∼600 steps the error grows linearly (fitted slopes of −3.6 to −3.7×10−5 per step over steps 1000 – 3000 , 3000 – 6000 , and the full window, R2>0.993 in each case), reaching −19.9% by step 6000 despite 99.88% of cells being handled by the exact solver throughout.
Figure 10: The directed bias generalises across progenitors. The frozen model is evaluated on the 20 ms post-bounce data of six independent progenitor collapses ( 12 to 40M⊙ ); the directed electron-fraction bias ∣bias∣/MAE stays at 0.60 to 0.85 , comparable to the 0.60 baseline (dashed), while the entropy-increment bias is smaller on the independent masses. Bootstrap over ten temporal regions.