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.
Neural surrogates promise large speedups over classical solvers for physical dynamics but fail silently at sharp dynamical events such as shocks, fronts, and contact. We present hybrid neural world models for physical dynamics: a recipe for training and deploying multi-horizon surrogates in physical state space, where a single network with continuous horizon conditioning is trained with direct supervision against textbook reference solvers to predict any future state at horizon T in one forward pass. Although no part of the training data, loss function, or architecture supervises discontinuity location, the trained surrogate encodes it implicitly, recoverable from its forward passes alone as a per-trajectory error map that concentrates on shocks, fronts, and contacts, and stays small elsewhere. The map is competitive with or better than standard label-free baselines including deep ensembles, learned error heads, gradient-magnitude indicators, and locally-adaptive conformal prediction, while using only a single trained network and requiring no calibration set or governing-equation knowledge. The recipe supports two operating points. Mode 1 runs the surrogate alone for maximum throughput, with same-hardware CPU speedups of 26x to 72x against textbook solvers on the PDE environments. Mode 2 uses the error map to gate a reference-solver fallback, deferring uncertain trajectories and roughly halving the surrogate's residual error at the default operating point. The recipe applies without modification across reaction-diffusion, compressible Euler, and rigid-body collision dynamics.
Neural PDE surrogates are trained on numerical solver outputs that contain both physical evolution and solver-specific discretization errors. Because surrogates are also evaluated against held-out trajectories from the same solver, standard benchmarks cannot distinguish fidelity to the exact evolution from imitation of the numerical scheme. We introduce an empirical Fourier-symbol diagnostic that probes a trained surrogate's linearized one-step operator with individual Fourier modes and compares it with both exact-evolution and training-scheme references. To address architectural spectral bias, we train identical networks on schemes with orthogonal dissipative and dispersive signatures and compare their learned operators. In linear advection, the learned surrogates reproduce the training schemes' amplitude and phase errors, with the twin-scheme difference reaching more than 99.8% of the analytically predicted full-imitation ceiling. The same behavior occurs for a non-local Fourier neural operator and at the operator level for nonlinear Burgers dynamics. These results show that agreement with solver-generated test data does not by itself establish fidelity to the exact evolution. Fourier-symbol measurements provide a direct diagnostic of numerical provenance.
Ridham Patel
Computer Science & Engineering Indian Institute of Technology Gandhinagar Gandhinagar, India
Neural surrogates for Partial Differential Equations (PDEs) are usually trained on synthetic data from numerical simulators. As models grow and scaling laws demand larger datasets, the cost of data generation often outweighs the training cost and becomes the bottleneck in the scaling quest. Unlike data in other domains, however, each PDE sample has a tunable cost that depends on its numerical fidelity, with the resulting cost-accuracy trade-off governed by known convergence rates. We exploit this property by formalizing neural surrogate training as a function of both data quantity and fidelity under a fixed data generation compute budget, treating fidelity as a fine-grained discrete spectrum rather than a binary distinction. Building on this, we show that multi-fidelity training enables Pareto-optimal error scaling against data generation compute. By using simulations at lower fidelities as an additional training signal, training on the optimal fidelity mixture matches the accuracy of high-fidelity-only training, with up to 6.3× less data generation compute across our settings. We validate these findings across six PDE problems spanning structured grids and unstructured meshes, two distinct error sources (discretization and iterative-solver error), and four surrogate architectures, totaling over 7,000 training runs.
Paul Setinek, Fabian Paischer, Nils Thuerey +2
Institute for Machine Learning, JKU Linz · Mistral AI · Technical University of Munich +1