A Two-Phase Deep Learning Framework for Adaptive Time-Stepping in High-Speed Flow Modeling
Authors: Jacob Helwig
Sai Sreeharsha Adavi
Xuan Zhang
Yuchao Lin
Felix S. Chim
Luke Takeshi Vizzini
Haiyang Yu
Muhammad Hasnain
Saykat Kumar Biswas
John J. Holloway
Narendra Singh
N. K. Anand
Swagnik Guhathakurta
Shuiwang Ji
Venue: ICLR 2026 Poster (2026)
Keywords: Partial Differential Equations
PDEs
CFD
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)

Abstract: We consider the problem of modeling high-speed flows using machine learning methods. While most prior studies focus on low-speed fluid flows in which uniform time-stepping is practical, flows approaching and exceeding the speed of sound exhibit sudden changes such as shock waves. In such cases, it is essential to use adaptive time-stepping methods to allow a temporal resolution sufficient to resolve these phenomena while simultaneously balancing computational costs. Here, we propose a two-phase machine learning method, known as ShockCast, to model high-speed flows with adaptive time-stepping. In the first phase, we propose to employ a machine learning model to predict the timestep size. In the second phase, the predicted timestep is used as an input along with the current fluid fields to advance the system state by the predicted timestep. We explore several physically-motivated components for timestep prediction and introduce timestep conditioning strategies inspired by neural ODE and Mixture of Experts. We evaluate our methods by generating three supersonic flow datasets, available at https://huggingface.co/divelab. Our code is publicly available as part of the AIRS library (https://github.com/divelab/AIRS).

Published as a conference paper at ICLR 2026
A TWO-PHASE DEEP LEARNING FRAMEWORK FOR ADAPTIVE TIME-STEPPING IN HIGH-SPEED FLOW
MODELING
Jacob Helwig'
Sai Sreeharsha Adavi²
Xuan Zhang'
Yuchao Lin' Muhammad Hasnain
Felix S. Chim
Luke Takeshi Vizzini
Haiyang Yu!
Saykat Kumar Biswas John J. Holloway
Narendra Singh²
N. K. Anand
Swagnik Guhathakurta Shuiwang Ji
Department of Computer Science and Engineering, Texas A&M University
Department of Aerospace Engineering, Texas A&M University
JJ. Mike Walker '66 Department of Mechanical Engineering, Texas A&M University
{jacob.c hel wiig, swagnik, sji ) aamataaaaee. edu
ABSTRACT
We consider the problem of modeling high-speed flows using machine learning methods. While most prior studies focus on low-speed fluid flows in which uni- form time-stepping practical, flows approaching and exceeding the speed of sound exhibit sudden changes such as shock waves. In such cases, it is essential to use adaptive time-stepping methods to allow temporal resolution sufficient resolve these phenomena while simultancously balancing computational costs. Here, we propose a two-phase machine learning method, known as ShockCast, model high-speed flows with adaptive time-stepping. In the first phase, we propose to employ a machine learning model to predict the timestep size. the second phase, the predicted timestep used as an input along with the current fluid fields to advance the system state by the predicted timestep. We explore several physically-motivttes components for timestep prediction and introduce timestep conditioning strategies inspired by neural ODE and Mixture of Experts. We evaluate our methods by generating three supersonic flow datasets, available https://huggingffcc.oo//// co/divellal Our code is publicly available as part of the AIRS library
INTRODUCTION
Learning fluid dynamics aims to accelerate fluid modeling using machine learning models (Li et al. 2021 Zhang al. 2023). Because this is emerging area of research, most current works fo- cus on low-speed scenarios in which flows are assumed to be incompressible. In such cases, the time scale of dynamics is relatively stable, enabling the use of time-stepping schemes with uniform step sizes without substantially affecting solution quality or the required computational effort. In contrast, time scales vary greatly for high-speed flows such that uniform time-stepping is longer a tractable strategy, For example, supersonic flow occurs when fluid moves faster than the lo- cal speed of sound (Anderson) 2023 2020). The speed of such flows is typically characterized by the Mach number MZ defined as the ratio of the flow velocity u to the local speed of sound as M v/a. Flows in the supersonie regime (commonly 1 M < 5) exhibit distinct phenomena with small time scales, including shock waves, expansion fans, and significant compressibility ef- fects (Anderson) 2020). Hypersonic flow refers extremely high-speed fluid flows, conventionally defined by Mach numbers greater than 5. In the hypersonic regime, encountered in the design of spacecraft, missiles, and atmospheric reentry vehicles, flows exhibit unique and complex behaviors such heating, strong shock wave interactions, and chemical reactions.
"qqual senior authorship

---

Published as a conference paper at ICLR 2026
For both supersonic and hypersonic flows, the time scale required accurately resolve these phe- nomena is much smaller than other parts of the dynamics. Therefore, uniform time-stepping is no longer practical, it would require the of the smallest time scale for all steps, inflating the required computation prohibitively. Instead, high-speed flow solvers employ adaptive time-stepping schemes which dynamically adjust the timestep size such that smaller steps are taken in the presence of sharp gradients. Adaptive time-stepping can also benefit neural solvers through more balanced objectives. When using uniform step size, the amount of evolution the model is required to learn can vary greatly between states with sharp gradients and smoother states, which is an especially pertinent consideration for high-speed flows. By instead inversely scaling the step size according to the rate of change, the difficulties of different training pairs are more evenly distributed.
However, because neural solvers achieve speedup through use of coarsened space-time meshes, classical approaches for determining the timestep size are not applicable. We therefore develop ShockCast, a machine learning framework for temporally-adaptin ShockCast consists of two phases: neural CFL phase, where the timestep is predicted, and neural solver phase, where the flow field is evolved forward in time by the predicted timestep size. We investigate the effect of physically-motivated motivated components in neural CFL model, and develop several novel timestep conditioning strategies for neural solvers inspired by neural ODE (Chen et al.) 2018) and Mixture of Experts (Shazeer et al.] 2017). To evaluate our framework, new high-speed flow datasets modeling a circular blast (maximum Mach numbers vary from 0.49 to 2.97 across cases), multiphase coal dust explosion (Mach numbers of the initial shock vary from 1.2 to 2.1), and shock sweeping over an airfoil (Mach numbers of the initial shock again vary from 1.2 2.1). Our work takes steps towards developing machine learning models for high-speed flows, where there is great potential for neural acceleration due to the immense computational requirements of classical methods.
modeling of high-speed flows.
generate three
2 BACKGROUND
In Section 2.1 we briefly introduce how PDEs are solved numerically before describing the role that the CFI. condition plays in this task in Section 2.2 We then overview several prominent machine learning approaches for solving PDEs in Section 2.3
2.1 SOLVING TIME-DEEPEDDEN PARTIAL DIFFERENTIAL EQUATIONS
Time-dependent PDEs are common in engineering, with some of the most prominent applications arising fluid dynamics. In two spatial dimensions, they typically equate first derivative in time to some operator H of spatial derivatives for an unknown solution function u(I, t) E IR as
0,u- = H(u, O₂ u, Urr 0, Oyy Ony ...),
(1)
with boundary conditions and initial conditions imposing additional constraints on u. To solve a PDE, we need to identify function u satisfying these constraints. In most real-world applications of PDE modeling, including fluid dynamics, an analytical form of this solution is intractable to obtain, and so we instead rely on producing u in numerical form, that is, obtaining point-wiss wise evaluations on a discrete set of collocation points in srace-time. This is done by evolving the PDE forward in time by first approximating the spatial derivatives on the right-hand side Eqquatio (1) using finite difference methods, finite volume methods, finite element methods (Reddy et al. 2022).
or spectral methods (Gottlieb Orszag] 1977] Gottlieb & Hesthaven, 200T Canuto et al. 2007
Kopriva 2009). Plugging these quantities into H. gives (tu(t), which is then time-integrated to advance the solution in time by a step size of At.
2.2 THE COURANT-RREDDICIS-LLLNE CONDITION
Numerical time integrators are very sensitive to the timestep size At. When u(t) is changing rapidly, more formally, when [|0,u(t)] grows large, too large At can lead to divergence of the numerical solution (Anderson 2023) or nonphysical solutions containing negative densities or pres- sures (Patkar et al. 2016). In low-s owespeed flows, the time scale does not vary drastically, which is to say that the magnitudes of temporal derivatives do not fluctuate substantially. In such cases, At can be chosen as fixed value to match the smallest time scale, simplifying the numerical solution
2

---

Published as a conference paper at ICLR 2026
process. On the other hand, time scales vary greatly in high-speed flows. For example, shock wave interactions in supersonic and hypersonic flows produce extremely sharp spatial gradients that can only be resolved with small timestep sizes. Following the dissipation of these phenomena, the lution can become smoother such that the time scale is substantially larger. Uniform time-stepping schemes, which must maintain a timestep small enough to resolve sharp gradients even in smooth gions, therefore impose greater computational burden compared to adaptive time-stepping methods, where timestep sizes are dynamically adjusted according the rate of change of the solution.
Adaptive time-stepping methods employ the Courant-Friedrich--Leyy (CFL) Condition (Courant al., [967) to determine the timestep size. The CFL condition is a necessary condition on the timestep size to attain convergence of the numerical solution (Bartels, 2016). For single-phase flow in two spatial dimensions and a target Courant number E (0, 1), the condition requires that
C Amax
At
min (Ar, Ay), 2.9
where ming.y(a,, Ay) is the minimum cell height and width in the spatial discretization and the maximum wave speed Amax is defined as
max X(x, v) 2.9
Mxx = max (lu(x, v\l vl.,,)), lo v)1 + a(x, v)),
(3)
max
where u and u denote the x and y components of the velocity, and a(x, y) = RT (z, y) is the local sound speed defined by the ratio of specific heats y, the specific gas constant R and the temper- ature T. Intuitively, the CFL condition restricts information flow for stability such that information propagates no more than one cell in any direction per timestep, which can be seen from the scaling by the minimum cell size and inverse scaling by the wave speed in Equation (2)
2.3 NEURAL SOLVERS
As previously discussed. numerically solving PDEs is a computationally intensive process. Deep surrogate models which can accelerate the solution of PDEs are therefore of great interest. Various approaches have emerged over the last decade, including Physics-Informed Neural Networks (Raissi
ct 2019). hybrid solvers (Kochkov et al. [202 Kaneda et al. |003) (Sun et al. (023), and operator learning (Kovachki et al., (2023). Speedups over classical methods are primarily achieved by the
ability of neural solvers learn solution mappings on coarsened grids in space-time (Stachenteld
et al. 2021; Kochkov et al. (02T). To maintain stability, classical methods require computational grids to be sufficiently fine in time, as specified by the CFL condition, as well as in space. On the other hand, neural methods can learn to map between solutions spaced hundreds of classical solver steps apart on much lower-dimensional spatial discretizations, thereby realizing substantial speedups. Additionally, classical methods require that all D flow variables comprising the PDE be evolved, whereas neural solvers can learn to explicitly model only a subset of variables of interest while implicitly learning incorporate the effect of the omitted variables.
METHODS
3.1 LEARNING HIGHI-PPEED FLOWS
High-speed flows one of the most resource-intensive applications of PDE modeling and therefore stand to benefit greatly from the speedup offered by neural solvers. In such settings, time-adaptive meshes present a more balanced one-step objective for neural solvers. When using a uniform step size At, the difference between inputs and targets || + At) u(t) tends to grow with the rate of change of inputs with respect to time ||0 (t)]l. which follows from the definition of derivative d.u(t)= == lim (u(t + h) /h. This implies that inputs with sharp gradients are less corre- lated with targets such that the mapping for charp-gradient inputs is more difficult to learn compared to smoother ones. By instead following adaptive scheme to inversely scale At according to the rate of change of (t)), the difference [u(t) - (ttttutttutt At] is more uniformly distributed across inputs with varying degrees of sharpness in gradients, thereby reducing variance in the training objective. This is a particularly relevant consideration for high-speed flows, where the sharpness of gradients varies greatly throughout time due phenomena such as shock waves. Furthermore, the ability
3

---

Published as a conference paper at ICLR 2026
Inputs Ground Truth 10
ShockCate and
Targis
Adaptive Time Simp-Sieppigg
Neural CFL
time
dif Ad)
år Ar)
wil Ae
Neural
Neural
Neural Solver
Nousal Selver
Nearal Solver
Training
Auforegressive Rollout
Figure I: Overview of the ShockCast framework for time-adaptive modeling of high-speed flows. Left: Training pipeline. The neural CFL model and time-conditioned neural solver are conditioned on the current flow state and predict the corresponding timestep size At and flow state At ahead, respectively. Right: Inference pipeline. ShockCast autoregressively alternates between predicting the timestep size given the current flow state using the neural CFL model and evolving the flow state forward in time by the predicted timestep size using the neural solver model. Note that the example data are from the circular blast dataset we generated in this work.
to train neural solvers time-adaptive meshes allows direct of solutions produced from high- speed flow solvers without introducing error by interpolating to a uniform grid or modifying solver codes to save solutions at uniform timesteps.
Although well-motivated, the use of adaptive temporal meshes presents several challenges for neural solvers. Because the step size determined by the solution, it is not known ahead of time and must instead be computed on-the-fly during autoregressive rollout. Importantly, the At at inference time must be aligned with those from training to avoid a test-time distribution shift. While the CFL condi- tion is used to determine the step size for classical solvers during data generation, it cannot produce At matching those in the training data due to the use coarsened computational meshes and only modeling a subset of the variables comprising the PDE. Specifically, while neural solvers are trained to advance time by hundreds of classical solver steps with single forward pass, the CFL condition used to compute the size of single solver step, and will therefore suggest step size orders of magnitude smaller than the neural solver encountered during training. Furthermore, because neural solvers learn coarsened spatial meshes, the cell sizes Ar and Ay appearing in Eqquatio (2) will not match those used to compute At in the training data. [qquttion (2) is the condition for a single phase flow, whereas multiphase setting that features, for example, solid phase interacting with a liquid phase has a much more complicated form involving large number of field variables.
Direct use of this form would require the neural solver learn to evolve all of these variables, thereby reducing the model's 's capacity to capture fields of interest accurately.
These challenges motivate us develop ShockCast, a two-phase framework consisting of conditop-o neural solver and a neural CFL model which can emulate the timestep sizes in the training data on a coarsened space-time mesh using only a subset of the field variables. At infer- ence time, each unrolling step utilizes each of the phases in turn. In the first phase, the neural CFL model predicts the timestep size which used by the neural solver in the second phase to evolve the current flow field forward in time by the predicted timestep. We investigate approaches for bet- ter aligning the neural CFL model with the CFL condition, and introduce several novel timestep conditioning strategies for the neural solver.
3.2 A TWO-PHASE FRAMEWORK
Our datasets D == (U.).V consist of numerical solutions to the compressible Navier-Stokes equations produced by a classical high-speed flow solver. Each solution U = {uy} consists of a series of n snapshots on a temporal grid T = {t where u, == u(t,) E R D denotes the solution at time point t, sampled on spatial discretization with M mesh points and D fields. Notably, T coarsened relative to the grid used by the classical solver by selecting every J-th
solution from the solver for coarsening factor J 2 100.
4

---

Published as a conference paper at ICLR 2026
For the first phase of our framework, train neural CFL model to minimize
ET~T.U~ (Ce (vb(w))
for the loss C.. where we take L. to be the MAE. In the second phase, we train a neural solver to map the solution at the current timestep u and the timestep size A, to the subsequent solution uj+t by optimizing the one-step objective given by
ETTTUUD [Cs A,), uz+]):
where we take C, to be the relative enror averaged over fields. While the neural solver is trained to emulate the behavior of the classical solver used to generate the data, the neural CFL model is trained emulate the process by which the timestep sizes are chosen while generating data. is important to again emphasize that due to the use of a coarsened computational mesh and reduced number of states being modeled, it not possible to directly use the CFL condition to deterministically predict the timestep size. At inference time, ShockCast predicts the full solution given the initial condition uo by alternating between the two phases autoregressively as
Åt (t))
útt Åt) - (ū((). At).
where út)) denotes the predicted state up until time This process, shown in Figure 1 repeated until a pre-specified stopping time is reached.
3.3 NEURAL CFL
Inspired by the classical CFL condition, we experiment with several modifications to the input fea- tures and internal structure of our neural CFL model. As adaptive time- timessteping schemes adjust At according to the sharpness of the gradients of ((), we include the spatial gradients Vu computed using finite differences for all fields in as inputs. From [Equation (2) we furthermore observe the dependence of the CFL condition on the max wave speed, computed by taking the maximum over the local wave speed A(z, y) in all computational cells. This operation can be viewed as a functional mapping À to the scalar value Xmaa via max pooling. This motivates us to employ max pooling as our spatial downsampling function. Finally, as previously discussed, the classical CFL condition is not directly applicable due to the use of mesh coarsening and modeling only a subset of the field variables. However, it is possible that the functions comprising the condition can be used to learn a surrogate condition. We therefore add CFL features inputs as the local wave speed X(ar, y), the velocity magnitudes |u(z, y and [0 y)l. and the local sound speed a(x, y).
3.4 TIMESTEP CONDITIONING FOR NEURAL SOLVERS
We now discuss our approaches for timestep conditioning for the neural solver phase.
Lime-Conditioned Layer Norm. Several prior works have considered training models to advance
time by multiples of a uniform step size (Gupta Brandstetter, 2023; Herde al. 2022). These models have utilized time-conditioned layer norm, technique originally introduced for condition- ing diffusion models on the diffusion time (Nichol & Dhariwal 2021 Dhariwal & Nichol 2021). Prior to each layer, the timestep size At is embedded into two vectors and b with sizes match- ing the hidden dimension dmodel of the feature map z. These are then applied as a scale and shift following each normalization layer LN(x) a) b.
Spatial-Spectral Conditioning. Many neural solvers perform convolutions in Fourier space and may not use normalization layers by default (Li et al. 2021 Tran al. 2023). For these models, the Spatial-Spectral conditioning strategy introduced by Gupta & Brandstetter (023) can be used to
perform timestep conditioning in the frequency domain. Under this scheme, the Fourier transform of the feature map is moint-wise multiplied with complex-valued embedding & of At as F(z)E. has different entries for each frequency of F(z), ) and so to maintain parameter-efficiency, Gupta & Brandstetter (2023) share & across all channels of z.
Euler Residuals. Neural solvers often employ residual connections He et al. 2016) as Z
21 + (=)),, where the 1-th solver layer F; includes spatial integration operations such convolution
5

---

Published as a conference paper at ICLR 2026
or attention. Many works have studied the relationship between residual connections and Euler integration (Lu et al. 2018 Haber & Ruthotto 2017 2017 Ruthotto & Haber 2020), and even extended
it to more general classes of integrators (Chen et al. Kidger 2022)
For function of time v,
Fuler integrators approximate time integration as
vt + At) = v(t) +
Qgu(T)dr ~ v(t) + Atd,v(t).
When viewing the evolution of the latent features H 21+1 from layer-to-layer as the evolution of some latent map z(t) E Rddaedd XM in time, the time integration of carried out numerically by the Euler integrator corresponds exactly to residual connections. In our scenario, we interpret the 1-th layer representation y corresponding to the input field u(t) and timestep size At as the latent form of some intermediate state u(t + adt), where E [0, 1J increases monotonically with depth in the network. This interpretation leads to Euler Residuals, in which the period of time integration executed by F) is related to the timestep size At as 2141 21 + (z)),, where is affine transformation of At as a = W + c for W.C [dmodde)
Mixture of Experts. Mixture of Experts (Shazeer al., 2017 Jacobs al., 1991 Fedus al.
2022) has emerged as an effective approach to scaling up models and has been applied in PDE modeling tasks (Hao et al. 2023). The Mixture of Experts (MoE) layer consists of gating network Gj and K experts Fi......, FLA-- The gating network is often simple MLP, while the experts can have more complex architectures. Based on layer inputs 2,, the gating network weighs the outputs
of the experts according to the gate (((z)) E IRK as 21+1 = 2 + Σx Fi,z The gating network can then learn to partition the latent space such that each expert specializes in particular area. As our task involves learning to evolve dynamics by variable time lengths, we gate experts according to the timestep Al
K
2111 21 + Σ Gi(At)x (L.FL.()))
where dg an affine transformation of At for the k-tth expert, and plays the same role as for the Euler residuals. This enables each layer to have experts specializing in short time integration periods where wtt) contains sharp gradients, as well as experts for handling longer timesteps where the dynamics behave more smoothly. Because our router weights are dense, we do not leverage the anddtionallcomputane
efficiency of sparse MoE architectures (Shazeer et al., 2017 Fedus et al., 2022). As result, the MoE conditioning mechanism trades increased model capacity for greater computational cost.
4 RELATED WORK
Learnable spatial re-meshing for PDEs has been an active area of research, with advances made in supervised Plaff et al. (202T); Song et al. (202)); Zhang et al. (2024a) and reinforcement learn- ing WWu et al. (202a)); Freymuth et al. (2023) Yang et al. (2023) frameworks. However, to the best of our knowledge, we the first to consider learning to temporally re-mesh by utilizing data with adaptable temporal resolution, which is vital for developing models for high-speed flows. Works in this direction have instead developed various schemes using temporally uniform data. |Wu et al.) 2025) recently introduced a pipeline wherein a timestep prediction model is trained using an unsu- pervised loss designed to avoid timesteps that are small. A second module then predicts temporal derivatives of various orders such that the solution can be queried using a Taylor expansion at any ground
point up until the predicted timestep. enabling supervision using the temporally-ume truth data. Following similar continuous-tim time strategy, Janny et al. (2024) propose to learn an in- terpolator such that arbitrary time points between temporally uniform training data can be queried. Hagnberger et al. (2024) employ conditional neural field to map from initial conditions to arbitrary query times. Other works learn to map forward in time by various multiples of a uniform step size, Lith|Liu et al. (2022a); Hamid et al. (2024) training different models for each step size and Cupta & Brandstetter (2023): y Herde et al. (2024) training one model to be shared across step sizes using time- conditioned layer norm. Importantly, Gupta & Brandstetter (2023); Herde al. (2024) consider the timestep size to be known a priori, with Gupta & Brandstetter (2023) treating this as a benchmark for probing the ability of neural solvers respond timestep conditioning and Herde et al. (2024)

---

Published as a conference paper at ICLR 2026
Neural CFI. One-Step Test Errors
Coal Dust Explosion
Circular Blast
0.425
0.19
00400
0.12
0.575
0,030
01025
Buse
Vu+CFL
Viu
Vu-Max Tu-Max-CUI
Vu+CIL
Vu
Vu-Max Vu-Max-CFL
Bass
Figure 2: One-step MAE Neural CFL models on At averaged over 3 training runs, where A/ normalized to have standard deviation 1. Error bars are + 2 standard emrors.
ShockCast: Unrolled At
Coal Dust Explosion, Shock Mach Number 1.85
Circular Blast, Max Mach Number 2.68
10
5.25
-9
5.00
10
4.75
80000
B.005
B.00
0.00155
0.0021
0.0025
0.0050
0.000
0.001
B.002
0,003
0.004
0.035
Predicted
True
Figure 3: Al predicted by autoregressive unrolling of ShockCast with F-rNO+Euler conditioning neural solver backbone for a selected solution.
using pre-training task. As detailed in (Section 3.1 this assumption is realistic in the setting we consider. Adaptive time-stepping extends beyond PDE modeling applications in ML re- search, with connections adaptive learning rates (Duchi al. |001] Tieleman 2012 Kingma & Ba 2015) and continuous-time generative models (Song et al. 2021 Lipman et al. 2023
EXPERIMENTS
5.1 DATASETS
We consider three settings of high-speed flows in the supersonic regime with two spatial dimen- sions, one of which we extend to three spatial dimensions. We discuss the generation of these cases in Appendix and visualize solutions in Appendix
Coal Dust Explosion. The first setting we consider is a multiphase problem containing both gaseous air and granular coal particles. The simulation begins with thin, uniform layer of coal dust settled on the bottom of a channel. Near the left boundary, we initialize normal shock. We vary the initial strength of the shock between Mach 1.2 and 2.1 along with the particle diameter from case to case for total of 100 cases, with 90 for training and 10 for evaluation. Once the simulation starts, the normal shock travels the right shown in Figure 10 where interacts with the dust layer, first compressing it and later generating instabilities at the gas-dust layer interface. These instabilities further grow with time into turbulent vortical structures, which raise the dust in the channel and mix them with the air. The amount of mixing depends both on the initial shock strength and the particle diameter. We train our models predict the velocity and temperature of the gas and the volume fraction describing the percentage of coal dust at each point.
Circular Blast. The second setting we consider is two-dimensional circular blast case, which represents a two-dimensional version the Sod's shock tube problem (Sod) 1978). We initialize circular region of high pressure such that the pressure inside the circle is substantially higher than its surroundings as shown in Figure 11 We vary the ratio these initial pressures from 1.99 to 50 to produce set of 99 cases split into 90 training cases and 9 evaluation cases. Once the simulation starts, a circular shock travels radially outward, while an expansion wave travels in the opposite direction. This continues until the outward moving shock reflects from the boundaries

---

Published as a conference paper at ICLR 2026
Circular Blast: Density
1,006 10
1,263
2,307
t= 3.779 10
4884 10
True
Predicted
Figure 4: Comparison of the ground truth (top) and predicted (bottom) density fields for the circular blast data. We obtain predictions with autoregressive unrolling of ShockCast.
ShockCast Results: Coal Dust Explosion
Correlation Time Proportion (1)
TKE Error (4)
Mean Flow Error (4)
0.100
0.62
0.12
0.60
0.000
0.10
0.085
F-FNO
U-Net
CNO Transolver
F-FN
U-Net
CNO Transolver
F-NO
U-Net
CNO
Transolver
Affine
Ealer
MoF
Figure 5: Coal dust explosion results averaged over three neural solver training runs. Error bars are + 2 standard errors.
and travels inwards toward the origin. The interaction of the reflected shocks with the post-shock gas generates instabilities which grow into turbulent structures. Once these reflected shocks reach the origin, they reflect again, thus propagating radially outward. This continues repeatedly, while with each reflection, the shocks lose strength. The maximum Mach number, which correlates with the initial pressure ratio, varied from 0.49 to 2.97 across cases. We train our models to predict the velocity, temperature, and density fields. To assess the scalability of ShockCast three spatial dimensions, we additionally generate and evaluate an extension of this setting to Spherical Blast. We additionally assess the stability of ShockCast across many time integrations by generating an alternative extension of this setting to Long Circular Blast by more than doubling the length of the simulation, resulting in trajectories with more than 100 unrolling steps.
Airfoil Shock. The third setting we consider is a shock-airfoil interaction problem involving a NACA 0012 airfoil immersed in compressible flow. We initialize normal shock inside the domain and vary its initial strength between Mach numbers 1.2 and 2.1. In addition, we vary the angle of attack of the airfoil between +-8 yielding total of 100 cases with 90 used for training and 10 for evaluation. Once the simulation begins, the shock travels from left to right and impinges on the airfoil, generating strong compression on the windward side followed by shock reflections and downstream vortical structures. The complexity of the flow increases with both the shock strength and the angle of attack. The irregularity of the airfoil geometry allows us evaluate ShockCast beyond regular domains with uniform spatial discretizations.
5.2 SHO BACKBONES
In the setting with regular domain and with uniform spacing, use the ConvNeXt architecture (Liu jet al. 20226) as the backbone architecture for our Neural CFL model, while in the Airfoil Shock settings, GNN with spatial downsampling between each message passing layer. Both
8

---

Published as a conference paper at ICLR 2026
ShockCast Results: Circular Blast
Correlation Time Proportion (f)
TKE Error (U)
Mean Flow Error (4)
0.99
0.0300
0.08
0.0273
0.98
0.06
0.0250
0.97
0.0225
0.04
0.0200
F-FNO
U-Net
CNO Transolver
F-FNO
U-Net
CNO Transolver
F-PNO
U-Net
CNO Transolver
Aftine
Buler
MoE
Figure 6: Circular blast results averaged over three neural solver training runs. Error bars are = standard errors.
ShockCast Results: Long Circular Blast
ShockCast Results: Airfoil Shock
Corr. Prop. (t)
TKE (4)
Mean Flow (4)
Com:. Prop. (t)
TKE (1)
Mcan Flow (1)
0.0466
0.95
6.080
0.1194
0.080
0.075
0.0464
0.1192
0.50
0.70
0.0462
D.066
0.1190
0.85
0.0160
0.020
0000
0,1188
PENO U-Net
F-ND U-Net
F-ND U-Net
MGN
DGN
MGN
DXGN
MGN
DGN
(a) Long circular blast results.
Airfoil shock results.
trained with the noise injection strategy from Gonzaez et al. (2020). We experiment with a variety of neural solver architectures. For the Coal Dust Explosion and Circular Blast settings, we experiment with the U-Net from Gupta & Braandsteterr (2023), Convolutional Neural Operator (CNO) (Raonic et al., 2024), Factorized FNO (Tran et al. 2023) (F-FNO), and Transolver (Wu al.) 2024). We refer to the Affine version of models the using either spatial-spectral conditioning in the case of F-FNO or time-conditioned layer norm in the case the remaining architectures. The remaining model variants add onto the Affine versions by additionally employing Fuler condition- ing, as well as MoE conditioning, both of which use Euler residuals. Finally, for the Airfoil Shock settings, we experiment with MeshGraphNets (MGN) (Pfaff et al.) (021) and Diffusion Graph Nets (DGN) from Valencia et al. 2025), which is based on the Graph U-Net (Lino et al. 2022 Gao & Ti 2019). Both GNN architectures use affine timestep conditioning. We discuss these architectures and training procedures in greater detail in Appendices B and respectively.
5.3
METRICS
To evaluate predicted solutions, we compute the Pearson's 's correlation coefficient for each field and at each timestep with the ground truth data. This requires the predicted fields to be on the same tem- poral grid as the ground truth data, for which we use linear interpolation, described in Anpppeddix
The correlation time (Kochkov al., 2021 Lippe al. 2023 Alkin al. 2024 for given field defined as the last time f before the correlation sinks below a threshold, which we take to be 0.9 here. We then average this time
all fields and report
percentage of the full simulation
Airfoil Shock: Density
3,405 x 10-
3,777 x 10-
7.365 10
(- 1,073 x 10-4
3,096 x 10-6
- 3,807 x 10-3
7.00 10
(- 1,077 x 10-1
Proticad
Figure 8: Comparison of the ground truth (top) and predicted (bottom) density fields for the Airfoil Shock data. We obtain predictions with autoregressive unrolling of ShockCast.

---

Published as a conference paper at ICLR 2026
time. We additionally the predicted fields to compute several of the primary physical quantities of interest to practitioners. The mean flow is computed by averaging each flow field over time, while the Turbulence Kinetic Energy (TKE) is calculated as the sum of the variances of the fluctuating part of the velocity field. These quantities are given by
1 = 7 J. u(t)dt
((((t)) TKE 2T / (u(1) - u)² + (v(t) - d)
(4)
respectively. We report the relative error averaged over each variable for the mean flow, and relative error of the TKE field, where the integrals approximated using the trapezoidal rule.
5.4 RESULTS
For the circular blast, single-phase problem where the CFL condition is determined entirely by the velocity and temperature fields, the base model yields the best performance, as the modeled variables by themselves sufficient accurately predict the timestep size. In contrast, the coal dust explosion has more complicated form of the CFL condition to account for both the solid phase describing the coal dust and the gas phase. The modeled variables primarily focus on the gas phase, with the exception of the volume fraction. In this more challenging setting, we observe substantial benefits from our physically-mottivtee motivated enhancements to the neural CFL model. Our best results are achieved when using max-poolee max-pooling, with the spatial gradient of the flow state Vu and CFL features added to inputs. In Figure 3] we examine the predicted At obtained through autoregressive unrolling of the ShockCast framework in each setting and observe a close match with the ground truth. We use the best neural CFL model in each of the settings for the full ShockCast models discussed next.
We examine the one-step MAE of several variants of the neural CFL model in Figure 2
We visualize unrolled predictions on the circular blast density field in Figure 4 In Figures 5 and 6 we examine the performance of ShockCast with each of the neural solver architectures and timestep conditioning strategies. For both the coal dust explosion and circular blast settings, ShockCast achieves the strongest performance in terms of correlation time when leveraging a U-Net backbone with moneicondtioone layer norm. On the coal dust explosion cases, MoE conditioning and Euler conditioning with a U-Net backbone achieve the first and second best performance in terms of TKE error, respectively. Similarly, the TKE error for the circular blast is lowest for Fuler and MoE conditioning strategies with a F-FNO backbone. Finally, while U-Net with con-ccondtiooee layer norm has the best mean flow error for the coal dust explosion setting, the F-FNO with MoE variant of ShockCast has the lowest circular blast mean flow eтror.
In Figure 9 we visualize unrolled predictions for spherical blast case using a 3D U-Net neu- ral solver backbone with time-conditionel layer norm and 3D ConvNeXt neural CFL backbone. In Figure 7a| we present results for the long circular blast setting for U-Net and F-FNO neural solver backbones with time-conditioned layer norm. The solutions remain correlated about 0.9 for approx- imately half of the simulation time, and the low TKE and mean flow errors demonstrate that the predicted solutions remain stable and in-distribution, which is further supported by the rollout visu- alizations shown in im[ppendix G Finally, we compare ShockCast with MGN and DGN neural solver backbones in Figure 7b] where find DGN to be more performant across metrics. We visualize ShockCast predictions for a test airfoil case in Figure 8 We present extended results in Appendix K
CONCLUSION
In this work, we develop machine learning methods for modeling high-speed flows with adaptive time-stepping. To this end, we propose ShockCast, a two-phase framework that learns timestep sizes in the first phase and evolves fluid fields by the predicted step size in the second phase. To evaluate ShockCast, we generate three new supersonic flow datasets. Results show that ShockCast effective learning temporally re-mesh and evolve fluid fields. Our work takes steps towards developing machine learning models for high-speed flows, where there is great potential for neural acceleration due the immense computational requirements of classical methods.
10

---

Published as a conference paper at ICLR 2026
ACKNOWLEDGMENTS
This work was supported in part by the National Science Foundation under grant IIS-2243850. The authors express their gratitude to Professor Ryan Houim of the University of Florida for providing access to HyBum, the computational fluid dynamics (CFD) code utilized for the simulations pre- sented in this study. The authors are grateful to Fernando Guzman for his assistance generating CFD data. The CFD calculations presented in this work were partly performed on the Texas A&M sesh-performnnc computing cluster Grace.
REFERENCES
Benedikt Alkin, Andreas Fürst, Simon Lucas Schmid, Lukas Gruber, Markus Holzleitner, and Jo- hannes Brandstetter. Universal physics transformers: A framework for efficiently scaling neural operators. In The eightt Annual Conference on Neural Information Processing Systems, 2024. URL httpps//ppeeeeeeveeeeeeffennnooo
John D. Anderson. Modern Compressible Flow: With Historical Perspective. McGraw-Hill Educa- tion, 4th edition, 2020.
John D. Anderson. Fundamentals of Aerodynamics. McGraw Hill, New York, 7th edition, 2023.
Dinshaw Balsara and Chi-Wang Shu. Monotonicity preserving weighted essentially non- oscillatory schemes with increasingly high order of accuracy. Journal of Computational Physics, 160(2):405-452, 2000.
Sören Bartels. Numerical approximation of partial differential equations, volume 64. Springer, 2016.
Kaifeng Bi, Lingxi Xie, Hengheng Zhang, Xin Chen, Xiaotao Gu, and Qi Tian. Pangu-weather: A 3d high-resolution model for fast and accurate global weather forecast. arXiv preprint arXiv:2211.02556, 2022.
Saykat Kumar Biswas and NK Anand. Three-dimensional laminar flow using physics informed deep neural networks. Physics of Fluids, 35(12), 2023.
Saykat Kumar Biswas and NK Anand. Interfacial conditioning in physies informed neural networks. Physics of Fluids, 36(7), 2024.
Boris Bonev, Thorsten Kurth, Christian Hundt, Jaideep Pathak, Maximilian Baust, Karthik Kashinath, and Anima Anandkumar. Spherical fourier neural operators: Learning stable dy- namics on the sphere. In Proceedings the 40th International Conference on Machine Learning, 2023.
Florent Bonnet, Jocelyn Ahmed Mazari, Paola Cinella, and Patrick Gallinari. Airfrans: High fi- delity computational fluid dynamics dataset for approximating reynolds-averaged navier-stokes solutions. In 36th Conference on Neural Information Processing Systems (NeurIPS 2022) Track on Datasets and Benchmarks, 2022.
Johannes Brandstetter, Daniel E. Worrall, and Max Welling. Message passing neural PDE solvers. In International Conference Learning Representations, 2022. URL https://openneview.nee// forum
Robert Bridson. Fast poisson disk sampling in arbitrary dimensions. SIGGRAPH sketches, 10(1):1, 2007.
Claudio Canuto, M Yousuff Hussaini, Alfio Quarteroni, and Thomas A Zang. Spectral methods: evolution to complex geometries and applications to fluid dynamics. Springer Science & Business Media, 2007.
Shuhao Cao. Choose a transformer: Fourier or galerkin. Advances in neural information processing systems, 34:24924-24999 2021.
11

---

Published as a conference paper at ICLR 2026
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. Advances in neural information processing systems, 31, 2018.
Woojin Cho, Kookjin Lee, Donsub Rim, and Noseong Park. Hypernetwork-based meta-learning for
low-rank physics-informed neural networks. Advances in Neural Information Processing Systems, 36, 2024.
Richard Courant, Kurt Friedrichs, and Hans Lewy. On the partial difference equations of mathemat- ical physics. IBM journal of Research and Development, 11(2):15-234, 1967.
Chengyuan Deng, Shihang Feng, Hanchen Wang, Xitong Zhang, Peng Jin, Yinan Feng, Qili Zeng, Yinpeng Chen, and Youzuo Lin. Openfwi: Large-scale multi-structural benchmark datasets for full waveform inversion. Advances in Neural Information Processing Systems, 35:6007-6020, 2022.
Prafulla Dhariwal and Alexander Nichol. Diffusion models beat gans on image synthesis. Advances in neural information processing systems, 11888888-794 2021.
Dimitris Drikakis, Marco Hahn, Andrew Mosedale, and Ben Thomber. Large eddy simulation using high-resolution and high-order methods. Philosophical Transactions of the Royal Society 2997, 2009.
A: Mathematical, Physical and Engineering Sciences,
John Duchi and Hongseok Namkoong. Variance- regularization with convex objectives. Jour- nal of Machine Learning Research, 20(68):1-55, 2019.
John Duchi, Elad Hazan, and Yoram Singer. Adaptive subgradient methods for online learning and stochastic optimization. Journal of machine learning reseurch, 12(7), 2011.
Anthony A Egeln, John C Hewson, Daniel R Guildenbecher, Ryan Marinis, Mare C Welliver, and Ryan W Houim. Post-detonation fireball modeling: Validation of freeze out approximations. Physics of Fluids, 35(6), 2023.
William A Falcon. Pytorch lightning. GitHub, 3, 2019.
T Farrukh, R Houim, D Guildenbecher, M Welliver, and S Balachandar. Particle and fluid time scales in a spherical multiphase blast flow. Shock Waves, pp. 1-19, 2025.
William Fedus, Barret Zoph, and Noam Shazeer. Switch transformers: Scaling to trillion parameter models with simple and efficient sparsity. Journal of Machine Learning Research, 23(120):1-39 2022.
Nico Fleischmann, Stefan Adami, and Nikolaus A Adams. A shock-stable modification of the hile riemann solver with reduced numerical dissipation. Journal of computational physics, 423: 109762, 2020.
Niklas Freymuth, Philipp Dahlinger, Tobias Daniel Würth, Simon Reisch, Luise Kärger, and Ger- hard Neumann. Swarm reinforcement learning for adaptive mesh refinement. In Thirly-seventh Conference Neural Information Processing Systems, 2023. URL https://ppeneeiee.ne//foonn id=:ZqRu8c4uc
Hongyang Gao and Shuiwang Ji. Graph u-nets. In international conference on machine learning, pp. 2083-2092. PMLR, 2019.
David Gottlieb and Jan S Hesthaven. Spectral methods for hyperbolic problems. Journal of Com- putational and Applied Mathematics, 128(1-2):83 1111113333-131 2001.
David Gottlieb and Steven Orszag. Numerical analysis of spectral methods: theory and applica- tions. SIAM, 1977.
Fernando F Grinstein, Len G Margolin, and William J Rider. Implicit large eddy simulation, vol- ume 10. Cambridge university press Cambridge, 2007.
Swagnik Guhathakurta. Effect of Radiative Heat Transfer on the Structure and Propagation of Layered Cost-Dust Explosions. PhD thesis, University of Florida, 2021.
12

---

Published as a conference paper at ICLR 2026
Swagnik Guhathakurta and Ryan W Houim. Influence of thermal radiation on layered dust explo- sions. Journal of Loss Prevention in the Process Industries, 72:104509, 2021.
Swagnik Gubathakurta and Ryan W Houim. Impact of particle diameter and thermal radiation on the explosion of dust layers. Proceedings of the Combustion Institute, 39(3):2905-299) 2914, 2023a.
Swagnik Guhathakurta and Ryan W Houim. Propagation and severity of coal-dust explosions and the effect of radiation in different channel lengths. In Proceedings of the 29th international col- loquium on the dynamics of explosions and reactive systems, 2023b.
Gaurav Gupta, Xiongye Xiao, and Paul Bogdan. Multiwavelet-base operator learning for dif- ferential equations. In A. Beygelzimer, Y. Dauphin, P. Liang, and J. Wortman Vaughan (eds.), Advances in Neural Information Processing Systems, 2021. URL https://ppeneevee..e////unm
id=LZDiWaC9CGL
Jayesh K Gupta and Johannes Brandstetter. Towards multi-spatiteemporal-saaal generalized PDE modeling. Transactions on Machine Learning Research, 2023. ISSN 2835-8856. URL. https:
Eldad Haber and Lars Ruthotto. Stable architectures for deep neural networks. Inverse problems, 34(1):014004, 2017.
Jan Hagnberger, Marimuthu Kalimuthu, Daniel Musekamp, and Mathias Niepert. Vectorized conditional neural fields: A framework for solving time-dependent parametric partial differ- ential equations. In Forty-first International Conference on Machine Learning, 2024. URL
Asif Hamid, Danish Rafiq, Shahkar Ahmad Nahvi, and Mohammad Abid Bazaz. Hierarchical deep learning-based adaptive time stepping scheme for multiscale simulations. Engineering Applica- tions of Artificial Intelligence, 133:108430, 2024.
Zhongkai Hao, Zhengyi Wang, Hang Su, Chengyang Ying, Yinpeng Dong, Songming Liu, Ze Cheng, Jian Song, and Jun Zhu. Gnot: A general neural operator transformer for operator learning. International Conference on Machine Learning, pp. 12556- 12569. PMLR, 2023.
Zhongkai Hao, Chang Su, Songming Liu, Julius Berner, Chengyang Ying, Hang Su, Anima Anand- kumar. Jian Song, and Jun Zhu. DPOT: Auto-regressive denoising operator transformer for large- scale PDE pre- trainngin In Forty-first International Conference on Machine Learning, 2024. URL
Joshua W Hargis, Anthony Egein, Ryan Houim, and Daniel R Guildenbecher. Visualization of post- detonation fireball flowfields and comparison to cfd modeling. Proceedings of the Combustion Institute, 40(1 (0(1-4):100 105230, 2024.
Amiram Harten, Peter D Lax, and Bram Leer. On upstream differencing and godunov-type schemes for hyperbolic conservation laws. SIAM review, 25(1 25(1)335-1, 1983.
Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recog- nition. In Proceedings of the IEEE conference on computer vision and pattern recognition, pp- 770-778, 2016.
Jacob Helwig, Xuan Zhang, Cong Fu, Jerry Kurtin, Stephan Wojtowytsch, and Shuiwang Ji. Group equivariant Fourier neural operators for partial differential equations. Proceedings of the 40th International Conference Machine Learning, 2023.
Jacob Helwig, Xuan Zhang, Haiyang Yu, and Shuiwang Ji. A geometry-aware message passing neural network for modeling aerodynamies over airfoils. arXiv preprint arXiv:2412.09399, 2024.
Maximilian Herde, Bogdan Raonić, Tobias Rohner, Roger Käppeli, Roberto Molinaro, Emmanuel de Bézenac, and Siddhartha Mishra. Poscidon: Efficient foundation models for pdes. arXiv preprint ar arXiv:2405.1 19101, 2024.
Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusion probabilistic models. Advances in neural information processing systems, 33:6840-6851 2020.
13

---

Published as a conference paper at ICLR 2026
Masanobu Horie and Naoto Mitsume. Physics-embedeee neural networks: Graph neural pde solvers with mixed boundary conditions. Advances in Neural Information Processing Systems, 35:23218- 23229, 2022.
Ryan W Houim and Kenneth K Kuo. A low-dissipation and tim-raccurate method for compressible multi-component flow with variable specific heat ratios. Journal of Compulational Physics, 230 (23):8527-8553, 2011.
Ryan W Houim and Elaine S Oran. A multiphase model for compressible granular-gaseous flows: formulation and initial tests. Journal of fluid mechanics, 789:166-2201 2016.
Robert A Jacobs, Michael Jordan, Steven Nowlan, and Geoffrey Hinton. Adaptive mixtures of local experts. Neural computation, 3(1):79-87, 1991.
Steeven Janny, Aurélien Bénéteau, Madiha Nadri, Julie Digne, Nicolas Thome, and Christian Wolf. EAGLE: Large-scale learning of turbulent fluid dynamics with mesh transformers. In Interna- tional Conference on Learning Representations, 2023. URL https//openreviewnett/onmm ?id=
miLX4QpsARI
Steeven Janny, Madiha Nadri, Julie Digne, and Christian Wolf. Space and time continuous physics simulation from partial observations. In The Twelfth International Conference on Learning Rep- resentations, 2024. URL http:://pppoeeveeeeaae/ooooddvdd
Ayano Kaneda, Osman Akar, Jingyu Chen, Victoria Alicia Trevino Kala, David Hyde, and Joseph Teran. A deep conjugate direction method for iteratively solving linear systems. In International Conference on Machine Learning, pp. 15720-15736. PMLR, 2023.
Tero Karras, Miik Aittala, Samuli Laine, Erik Härkönen, Janne Hellsten, Jaakko Lehtinen, and Timo Aila, Alias-free generative adversarial networks. Advances in neural information processing systems, 34:852-863 2021.
Patrick Kidger. On neural differential equations. arXiv preprint arXiv:2202.02433. 2022.
Diederik Kingma and Jimmy Ba. Adam: method for stochastic optimization. In International Conference on Learning Representations, 2015.
Dmitrii Kochkov, Jamie A Smith, Ayya Alieva, Qing Wang, Michael P Brenner, and Stephan Hoyer. Machine learning-acceleraaeed computational fluid dynamics. Proceedings of the National Academy of Sciences, (18(1)):21017844 18, 2021.
David Kopriva. Implementing spectral methods for partial differential equations: Algorithms for scientists and engineers. Springer Science & Business Media, 2009.
Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, An- drew Stuart, and Anima Anandkumar. Neural operator: Learning maps between function spaces with applications PDEs. Journal of Machine Learning Research, 24(89):1-97, 2023. URL
Remi Lam, Alvaro Sanchez-Gonzalez Matthew Willson, Peter Wirnsberger, Meire Fortunato, Alexander Pritzel, Suman Ravuri, Timo Ewalds, Ferran Alet, Zach Faton-Rosen, et al. Graphcast: Learning skillful medium-range global weather forecasting. arXiv preprint arXiv:22/2.12794 2022.
Hsiao- Cha Li. Flame Quenching Accelerating by Metal Foam in Square Channel Using Im- mersed Boundary Method. University of Florida, 2022.
Hsiao-Chi Li and Ryan W Houim. Pore-scale resolved simulation of quenching, acceleration, and transition to detonation hydrogen explosions by metal foams. Combustion and Flame, 259: 113118, 2024.
Zijie Li, Kazem Meidani, and Amir Barati Farimani. Transformer for partial differential equations' operator learning. Transactions on Machine Learning Research, 2023a. ISSN 2835-8856. URL
14

---

Published as a conference paper at ICLR 2026
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, An- drew Stuart, and Anima Anandkumar. Neural operator: Graph kernel network for partial differ- ential equations. arXiv preprint Stt:::20000034855, 2020a.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Andrew Stuart, Kaushik Bhat- tacharya, and Anima Anandkumar. Multipole graph neural operator for parametric partial differ- ential equations. Advances in Neural Information Processing Systems, 33:6755-6766 2020b.
Zongyi Li, Nikola Borislavov Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhat- tacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential equations. In International Conference on Learning Representations, 2021. URL
http:://ppnneevie.ee///oommmm/ http://rrrunniiiooP/o//mmmmm
Zongyi Li, Daniel Zhengyu Huang, Burigede Liu, and Anima Anandkumar. Fourier neural operator with learned deformations for pdes on general geometries. arXiv preprint arXiv:2207.05209, 2022a.
Zongyi Li, Miguel Liu-Schiaffini, Nikola Borislavov Kovachki, Burigede Liu, Kamyar Azizzade- nesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Learning dissipative dynamics in chaotic systems. In Advances in Neural Information Processing Systems, 2022b.
Zongyi Li, Nikola Borislavov Kovachki, Chris Choy, Boyi Li, Jean Kossaifi, Shourya Prakash Otta, Mohammad Amin Nabian, Maximilian Stadler, Christian Hundt, Kamyar Azizzadenesheli, et al. -eooeery-informeed neural operator for large-scale 3d pdes. arXiv preprint urXiv:2309.0058: 2023b.
Mario Lino, Stathi Fotiadis, Anil A Bharath, and Chris D Cantwell. Multi-scale rotation-equivariint graph neural networks for unsteady eulerian fluid dynamics. Physics of Fluids, 34(8), 2022.
Yaron Lipman, Ricky T. Q. Chen, Heli Ben-Hamu, Maximilian Nickel, and Matthew Le. Flow matching for generative modeling. In The Eleventh International Conference on Learning Repre- sentations, 2023. URL httpss//ppppooooooo/ee////uu///n
Phillip Lippe, Bastiaan S Veeling, Paris Perdikaris, Richard E Turner, and Johannes Brandstet- ter. Pde-refiner: Achieving accurate long rollouts with neural pde solvers. arXiv preprint arXiv:2308.05773: 2023.
Yuying Liu, J Nathan Kutz, and Steven L Brunton. Hicrarchical deep learning of multiscale dif- ferential equation time-steppers. Philosophical Transactions of the Royal Society A, 380(2229); 20210200, 2022a.
Zhuang Liu, Hanzi Mao, Chao- Yuan Wu, Christoph Feichtenhofer, Trevor Darrell, and Saining Xie. A convnet for the 2020s. In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pp. 11976-11988 2022b.
Ila Loshchilov and Frank Hutter. Sgdr: Stochastic gradient descent with warm restarts. In Interna- tional Conference on Learning Representations, 2017.
Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang, and George Em Karniadakis. Learning nonlinear operators via deeponet based on the universal approximation theorem of operators. Nature machine intelligence, 3(3):218-22) 2021.
Yiping Lu, Aoxiao Zhong, Quanzheng Li, and Bin Dong. Beyond finite layer neural networks: Bridging deep architectures and numerical differential equations. In International conference machine learning, pp. 3276-3285. PMLR, 2018
M Pino Martín, Ellen M Taylor, Minwei Wu, and V Gregory Weirs. A bandwidth-optimizzed weno scheme for the effective direct numerical simulation of compressible turbulence. Journal of Com- putational Physics, 220(1 (20(1):20 11:2002288 2006.
Nick McGreivy and Ammar Hakim. Weak baselines and reporting biases lead to overoptimism in machine learning for fluid-related partial differential equations. Nature machine intelligence, 1106 2024.
15

---

Published as a conference paper at ICLR 2026
Tung Nguyen, Johannes Brandstetter, Ashish Kapoor, Jayesh K Gupta, and Aditya Grover. Climax: A foundation model for weather and climate. In Proceedings of the 40th International Conference on Machine Learning. 2023.
Alexander Quinn Nichol and Prafulla Dhariwal. Improved denoising diffusion probabilistic models. In International conference machine learning, pp. 8162-8171. PMLR, 2021.
Elaine S Oran, Jay P Boris, and Jay P Boris. Numerical simulation of reactive flow, volume 2. Citeseer, 2001.
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al. Pytorch: An imperative style, high- performance deep learning library. Advances in neural information processing systems, 32, 2019.
Jaideep Pathak, Shashank Subramanian, Peter Harrington, Sanjeev Raja, Ashesh Chattopadhyay, Morteza Mardani, Thorsten Kurth, David Hall, Zongyi Li, Kamyar Azizzadenesheli, et al. Four- castnet: A global data-driven high-resolution weather model using adaptive fourier neural opera- tors. arXiv preprint arXiv:2202. 12 2022.
Saket Patkar, Mridul Aanjaneya, Wenlong Lu, Michael Lentine, and Ronald Fedkiw. Towards positivity preservation for monolithic two-way solid -luid coupling. Journal of Computational Physics, 312:82-114, 2016.
Tobias Pfaff, Meire Fortunato, Alvaro Sanchez-Gonzakee and Peter Battaglia. Learning mesh-based simulation with graph networks. In International Conference Learning Representations, 2021. XP
URL https://ponneevvvv..ee//u/uutlop
Michael Poli, Stefano Massaroli, Federico Berto, Jinkyoo Park, Tri Dao, Christopher Re, and Ste- fano Ermon. Transform once: Efficient operator learning in frequency domain. Alice H. Oh, Alekh Agarwal, Danielle Belgrave, and Kyunghyun Cho (eds.). Advances in Neural Information Processing Systems, 2022. URL
Jacob W Posey, Brayden Roque, Swagnik Guhathakurta, and Ryan W Houim. Mechanisms of prompt and delayed ignition and combustion of explosively dispersed aluminum powder. Physics of Fluids, 33(11), 2021.
Ilan Price, Alvaro Sanchez-Gonzale Ferran Alet, Timo Ewalds, Andrew El-Kadi, Jacklynn Stott, Shakir Mohamed, Peter Battaglia, Remi Lam, and Matthew Willson. Gencast: Diffusion-basee ensemble forecasting for medium-range weather. arXiv preprint arXiv:2312.15796 2023.
Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686-707, 2019.
Bogdan Raonic, Roberto Molinaro, Tim De Ryck, Tobias Rohner, Francesca Bartolucci. Rima Alai- fari, Siddhartha Mishra, and Emmanuel de Bézenac. Convolutional neural operators for robust and accurate learning of pdes. Advances in Neural Information Processing Systems, 36, 2024.
Junuthula Narasimha Reddy, NK Anand, and Pratanu Roy. Finite element and finite volume methods for heal transfer and fluid dynamics. Cambridge University Press, 2022.
Olaf Ronneberger, Philipp Fischer, and Thomas Brox. U-net: Convolutional networks for biomed- ical image segmentation. In Medical Image Computing and Computer- Intervention- MICCAI 2015: 18th International Conference, Munich, Germany, October 5-9, 2015, Proceed- ings, Part III 18, pp. 234-241. Springer, 2015.
Lars Ruthotto and Eldad Haber. Deep neural networks motivated by partial differential equations. Journal of Mathematical Imaging and Vision, 62(3):352-364, 2020.
Alvaro Sanchez-Gonzaled Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovee, and Peter Battaglia. Learning to simulate complex physics with graph networks. In International conference machine learning, pp. 8459-8468 PMLR, 2020.
16

---

Published as a conference paper at ICLR 2026
Jacob Seidman, Georgios Kissas, Paris Perdikaris, and George Pappas. Nomad: Nonlinear man- ifold decoders for operator learning. Advances in Neural Information Processing Systems, 35: 5601-5613, 2022.
Smruti Shah and NK Anand. Physics-inforee neural networks for periodic flows. Physics of Fluids, 36(7). 2024.
Noam Shazeer, * *Azalia Mirhoseini, *Krzysztof Maziarz, Andy Davis, Quoe Le, Geoffrey Hinton, and Jeff Dean. Outrageously large neural networks: The sparsely-gated mixture-of-experts layer. In International Conference on Learning Representations, 2017. URL. httpss//openneviewnnell forum
Zhijun Shen, Wei Yan, and Guangwei Yuan. A robust hilc-type riemann solver for strong shock.
Journal of Computational Physics, 309:185-2000 2016.
Gary A Sod. A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws. Journal of computational physics, 27(1):1-31, 1978.
Wenbin Song, Mingrui Zhang, Joseph G Wallwork, Junpeng Gao, Zheng Tian, Fanglei Sun, Matthew Piggott, Junqing Chen, Zuoqiang Shi, Xiang Chen, al. M2n: Mesh movement net- works for pde solvers. Advances in Neural Information Processing Systems, 35:7199-7210, 2022.
Yang Song, Jascha Sohl- Sohh-Dickstein, Diederik Kingma, Abhishek Kumar. Stefano Ermon, and Ben Poole. Score-based generative modeling through stochastic differential equations. In Interna- tional Conference on Learning Representations, 2021. URL https://opentrrii..ee//ouum iid=
PxTIGI2RRHS
Kimberly Stachenfeld, Drummond B Fielding, Dmitrii Kochkov, Miles Cranmer, Tobias Pfaff, Jonathan Godwin, Can Cui, Shirley Ho, Peter Battaglia, and Alvaro Sanchez-Gonzalez Learned coarse models for efficient turbulence simulation. In International Conference on Learning Rep- resentations, 2021.
Zhiqing Sun, Yiming Yang. and Shinjae Yoo. A neural pde solver with temporal stencil modeling. In International Conference on Machine Learning, pp. 33135-33115, PMLR, 2023.
Ben Thornber, Andrew Mosedale, and Dimitris Drikakis. On the implicit large eddy simulations of homogeneous decaying turbulence. Journal of Computational Physics, 226(2): 226(2):1902-192) 1902-1929, 2007.
Ben Thomber, Andrew Mosedale, Dimitris Drikakis, David Youngs, and Robin JR Williams. An improved reconstruction method for compressible flows with low mach number features. Journal of computational Physics, 227(10):483-4894 2008
Tijmen Tieleman. Lecture 6.5-rmsprop: Divide the gradient by a running average of its recent magnitude. COURSERA: Neural networks for machine learning, 4(2):26, 2012.
Eleuterio F Toro. Riemann solvers and numerical methods for fluid dynamics: practical introduc- fion. Springer Science Business Media, 2013.
Alasdair Tran, Alexander Mathews, Lexing Xie, and Cheng Soon Ong. Factorized fourier neural operators. The Eleventh International Conference on Learning Representations, 2023. URL
http:://peneeiie..eeffommmm hit.:rrrmmmii--mmmiMPPia.....
Mario Lino Valencia, Tobias Pfaff, and Nils Thuerey. Learning distributions of complex fluid sim- ulations with diffusion graph networks. In The Thirteenth International Conference on Learning Representations, 2025.
Bram Van Leer. Towards the ultimate conservative difference scheme. second-order sequel to godunov's method. Journal of computational Physics, 32(1):101-136, 1979.
Gege Wen, Zongyi Li, Kamyar Azizzadenesheli, Anima Anandkumar, and Sally M Benson. U- fno-an an enhanced fourier neural operator-basee based deep-lcarning model for multiphase flow. Ad- vances in Water Resources, 163:104180 2022.
17

---

Published as a conference paper at ICLR 2026
Gege Wen, Zongyi Li, Qirui Long, Kamyar Azizzadenesheli, Anima Anandkumar, and Sally M Benson. Real-time high-resolution co 2 geological storage prediction using nested fourier neural operators. Energy & Environmental Science, 16(4):1732-1741) 2023.
Haixu Wu, Huakun Luo, Haowen Wang, Jianmin Wang, and Mingsheng Long. Transolver: A fast transformer solver for PDEs on general geometries. In Ruslan Salakhutdinov, Zico Kolter, Katherine Heller, Adrian Weller. Nuria Oliver, Jonathan Scarlett, and Felix Berkenkamp (eds.), Proceedings of the 41st International Conference on Machine Learning, volume 235 of Pro- ceedings of Machine Learning Research, pp. 53681-53705. PMIR, 21-27 Jul 2024. URL
Tailin Wu, Takashi Maruyama, Qingqing Zhao, Gordon Wetzstein, and Jure Leskovec. Learn- ing controllable adaptive simulation for multi-scale physics. In NeurIPS 2022 Al for Science: Progress and Promises, 2022a. URL https://pppmeeeiiinenee//ouddcdp
Tailin Wu, Qinchen Wang, Yinan Zhang, Rex Ying, Kaidi Cao, Rok Sosic, Ridwan Jalali, Hassan Hamam, Marko Maucec, and Jure Leskovec. Learning large-scale subsurface simulations with hybrid graph network simulator. In Proceedings of the 28th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, pp. 4184-4194, 2022b.
Zhikai Wu, Shiyang Zhang, Sizhuang He, Sifan Wang, Min Zhu, Anran Jiao, Lu Lu, and David van Dijk. Tante: Time-adaptive operator learning via neural taylor expansion. arXiv preprint DrXivv2502.08877 2025.
Jiachen Yang, Tarik Dzanic, Brenden Petersen, Jun Kudo, Ketan Mittal, Vladimir Tomov, Jean- Sylvain Camier, Tuo Zhao, Hongyuan Zha, Tzanio Kolev, et al. Reinforcement learning for adaptive mesh refinement. In International conference on artificial intelligence and statistics, pp- 5997-6014. PMLR, 2023.
Mingrui Zhang, Chunyang Wang, Stephan C Kramer, Joseph G Wallwork, Siyi Li, Jiancheng Liu, Xiang Chen, and Matthew Piggott. Towards universal mesh movement networks. Advances in Neural Information Processing Systems, 37:14934-14961 2024a.
Weiqun Zhang, Ann Almgren, Vince Beckner, John Bell, Johannes Blaschke, Cy Chan, Marcus Day, Brian Friesen, Kevin Gott, Daniel Graves, Max Katz, Andrew Myers, Tan Nguyen, Andrew Nonaka, Michele Rosso, Samuel Williams, and Michael Zingale. AMReX: A Framework for Block-Structured Adaptive Mesh Refinement. Journal of Open Source Softwre, 4(37):1370, 2019. doi: 10.21105/jos..11370
Xuan Zhang, Limei Wang, Jacob Helwig, Youzhi Luo, Cong Fu, Yaochen Xie, .... and Shuiwang Ji. Artificial intelligence for science in quantum, atomistic, and continuum systems. arXiv preprint arXiv:2307.0842.. 2023.
Xuan Zhang, Jacob Helwig, Yuchao Lin, Yaochen Xie, Cong Fu, Stephan Wojtowytsch, and Shui- wang Ji. Sinenet: Learning temporal dynamics in dimendepn partial differential equa- tions. In The Twelfth International Conference on Learning Representations, 2024b. URL
Xuan Zhang, Haiyang Yu, Chengdong Wang, Jacob Helwig, Xiaofeng Qian, and Shuiwang Ji. Orbital transformers for predicting wavefunctions in time-dependent density functional the- ory. In The Fourteenth International Conference on Learning Representations, 2026. URL
18

---

Published as a conference paper at ICLR 2026
Appendix
CFD Problem Description and Numerical Models
19
A.1
Coal Dust Explosion
20
Circular Blast
21
Spherical Blast
22
A.4 4 Airfoil Shock
22
B Neural Solvers and Neural CFL Models
23
Training Details
26
Datasets
26
Dataset Statistics
................................
26
Training Pipeline
26
C.4
Parameter Counts, FLOPs and Peak GPU Memory
28
C.5
Solution Runtime Analysis
28
D Metries
31
Hypersonic Airfoil Shock
31
F Mass Conservation Analysis
31
G Long Circular Blast Analysis
35
H Neural CFL Ablation
35
Feature Importance Analysis
41
Solution Visualizations
42
K Extended Results
59
Limitations and Future Directions
59
Broader Impacts
64
A CFD PROBLEM DESCRIPTION AND NUMERICAL MODELS
We conduct the numerical simulations in this work using the HyBurn code, which implements high- order Godunov schemes to solve the governing equations. HyBurn is finite volume solver that employs parallelization and adaptive mesh refinement (AMR) via the AMReX library (Zhang al.,
2019). It uses a low- cow-dissipation, WENO-based high-order Godunov method (Houim & Kuo 2011
Houim & Oran 2016 Balsara Shu 2000 Martin et al. 2006,,Shen al. 2016 Harten al. 1983
and models turbulence using Implicit Large Eddy Simulation (Grinstein et al. 2007 Oran et al.
2001 Thormber al. 20088 Drikakis al. 2009 Thornber et al. 2007p. HyBurn solves compress-
ible Navier-Stokes equations for Cnllerian-Fulerian granular multiphase reactive flows, employing
advanced numerical techniques developed by Houim & Oran (2016). It uses an explicit third-order three-stage Kunge-Kutta time-stepping algorithm along with two independent CFL numbers one for the hyperbolic and the other for the parabolic terms of the compressible Wavier-Stoker equa- tions. The actual timestep size is determined by whichever is the smaller of the two. HyBurn has been extensively validated and verified against various compressible multiphase and reactive flow problems (Houim & Kuo, 2011 Houim & Oran 2016). More details about the algorithms and
19

---

Published as a conference paper at ICLR 2026
Spherical Blast: Density
t -3.117 x 10-5 [ 1.293 x 10-1
L 2.340 x 10
1 3.443 x 10-1
1 4.517 x 10-1
True
t= 3.122 x 10-5
{ = 1.294 X 10-1
t = 2.342 x 10-1
1 == 3.446 x 10-1
t t= 4.519 x 10-1
Predicted
Figure 9: Comparison of the ground truth (top) and predicted (bottom) density fields for the spherical blast data, We obtain predictions with autoregressive unrolling of ShockCast
Coal Dust Explosion: Grid Independence
Circular Blast: Grid Independence
Refinement Levels
Shock speed (ms=1)
Refinement Levels
Time (s x 10-3)
1
1, 1,102.52
1.0178
2
10033.99
2
1.0179
3
999.70
3
1.0196
(a) Coal dust explosion grid independence tests. We choose arbitrary point in the domain and measure the shock speed at this point for each of the configu- rations.
(b) Circular blast grid independence tests. We choose an arbitrary point in the domain and measure the time taken for the shock to reach this point for each of the three configurations.
Table I: Grid independence tests. In each setting, we re-ran a selected case for varying levels of mesh refinement.
the applications of Hy can be found in previous work (Guhathakurta Houim, 2023a
2021
Guhathakurta 2021 Guhathakurta & Houim 2023b Posey et al. 2021: Li & Houim 2024
2022
Hargis et al. 2024 Egeln et al. [023] Farrukh et al. 2025).
A.1
COAL DUST EXPLOSION
The first setting we consider is loosely based on the coal dust explosion simulations in Guhathakurtal & Houim (2023a) 202 Guhathakurta (2021); Guhathakurta & Houim (2023b). The major dif- ferences are that in this study we use a normal shock instead of a detonation to mimic the primary explosion in coal mine, and there we use is a 25 em by 5 cm two-dimensional rectangular channel containing air, with a thin uniform layer of coal dust settled on the bottom of the channel. We keep the left and right boundaries open, while the top and bottom boundaries are symmetry. This is a multiphase problem containing both gaseous air and granular coal particles. Near the left boundary, we initialize a normal shock using post-shock conditions to the left of the shock and quiescent air (300 K temperature, 1 atm pressure) ahead of it. We vary the initial strength of the shock between Mach 1.2 and 2.1 along with the parti- cle diameter between um and 150 um from case to case for total of 100 cases. We set the initial volume fraction of the dust layer to 47% and the particles to monodisperse. We model the coal particles being composed of inert ash only, with solid- -paaee density of 1330 kg m -3, We use two AMR levels to give an effective resolution of ~0.12 12 mm at the finest level. We use a fifth-order WENO interpolation scheme and the HLL (oro 2013) flux reconstruction method. We set the CFL numbers for both the hyperbolic and parabolic terms to 0.6.
chemical reactions involved. The computational geometry
Once the simulation starts, the normal shock travels the right shown in Figure 10 where it interacts with the dust layer. first compressing it and later generating instabilities at the gas-dust layer interface. These instabilities further grow with time into turbulent vortical structures, which
20

---

Published as conference paper at ICLR 2026
Coal Dust Explosion: Gas Velocity x-component, Shock Mach Number 1.85
<= 01,000
1,137
[= 2,188
== 3.128
== 4.048
4.967
5.889
= 6.814
6= 7.736
-TPP
10
Figure 10: Initial gas velocity x-component for selected coal dust explosion case. Times are in units of 10 seconds and the downsampling factor relative to the classical solver solution is x00x compared to 500x used for training ShockCast. The initial shock can be seen to be moving from left to right.
Circular Blast: Initial Densities
IPR 3.44
IPR 10.23
IPR 16.53
IPR 23.81
IPR 29.63
IPR 33.51
IPR 37.88
JPR 43.68
IPR 47.58
30
40
10
Figure 11:
Initial density field for the circular blast evaluation cases with varying Initial Pressure Ratios (IPR). Due to the ideal gas law, increasing the pressure also increases the density.
raise the dust in the channel and mix them with the air. The amount of mixing depends both the initial shock strength and the particle diameter. We cap the total simulation time at 3 ms for all cases.
We used a case with initial strength of the shock Mach and particle diameter 30 um to conduct a grid independence analysis. We used base grid of 520 104 and conducted the analysis by varying the AMReX maximum levels of refinement between I and 3. In Table Ta we compare the shock speed at arbitrary point in the domain for each of the three configurations. Because the shock speed differed least between 2 and 3 levels of refinement (3.4%), we used 2 levels of refinement as the grid configuration for data generation.
A.2 CIRCULAR BLAST
The second setting we consider is a two-dimensional circular blast case, which represents a two- dimensional version of the Sod's shock tube problem (Sod) 1978). We initialize circular region of high pressure such that the pressure inside the circle is substantially higher than its surroundings as shown in Figure II We vary the ratio of these initial pressures from 1.99 to 50 to produce set of 99 cases. To reduce the computational cost, we only simulate a quarter of the circle. We use symmetry boundary conditions for all boundaries to allow for the generated shocks and expansion waves to reflect when incident on them. We model a single gaseous phase throughout the com- putational domain. We set the initial temperatures 300 K everywhere, with the gas initially rest. We use two Adaptive Mesh Refinement (AMR) levels triggered by predefined pressure and density ratio thresholds between any two adjacent computational cells give an effective resolu- tion of ~0.98 mm at the finest level, which is sufficient to resolve the shocks. We use fifth-order MUSCL Van Leer 1979) interpolation scheme and the HLLC-LM (Fleischmann et al.) 2020) flux reconstruction method. We set the CFL numbers for both the hyperbolic and parabolic terms to 0.8.
Once the simulation starts, a circular shock travels radially outward, while an expansion wave travels in the opposite direction. This continues until the outward moving shock reflects from the boundaries and travels inwards toward the origin. The interaction of the reflected shocks with the post-shock gas generates instabilities which grow into turbulent structures. Once these reflected shocks reach the origin, they reflect again, thus propagating radially outward. This continues repeatedly, while with each reflection, the shocks lose strength. We cap the total simulation time 5 ms for all cases.
We used case with blast pressure of 50 atm to conduct the grid independence test. We used a base grid of 256 x 256 and performed the grid independence test by varying the AMREX maximum levels of refinement between 1 and 3. In Table lb, we compare the time taken by the shock to reach an arbitrary point in the domain for each of the three grid configurations. As the difference in the
21

---

Published as a conference paper at ICLR 2026
shock speeds minor across the three configurations, we chose the grid with 2 levels of maximum refinement for generating the cases.
A.3 SPHERICAL BLAST
The spherical blast setting extends the circular blast configuration to three spatial dimensions. As in the 2D case, we initialize region of high pressure embedded in quiescent ambient air. The ambient state is set to 300 K and atm, with the gas initially at rest throughout the domain. Inside the spherical region, we increase the pressure to obtain initial pressure ratios ranging from 2 to 34 across cases while keeping the temperature fixed at 300 K, leading to correspondingly higher initial densities through the ideal gas law. We model single gaseous phase governed by a constant-y ideal-gas equation of state with the same thermodynamic parameters as in the circular blast setting,
To reduce computational cost, we exploit symmetry and simulate only one octant of the spherical blast. The computational domain is a cubic box of side length 25 cm, and we place the blast center the corner of the cube so that the coordinate planes coincide with symmetry boundaries. We employ symmetry boundary conditions on all faces, allowing the autward-propggating spherical shock and inward-moving expansion waves to reflect in a physically consistent manner. We use base grid of 64 64 x 64 cells with two Adaptive Mesh Refinement (AMR) levels, triggered by pressure and density gradients as in the circular blast case. This yields an effective finest-level cell size of approximately 0.98 mm, which is sufficient to resolve the shock structure. We use the MUSCL (Van Leer 1979) reconstruction scheme with the HLLC-LM flux method (Fleischmann) et al. 2020p, and set the CFL numbers for both hyperbolic and parabolic terms to 0.8. We each simulation until 2 ms, capturing the initial blast, the formation of the outgoing spherical shock and inward rarefaction, and their early reflections from the symmetry planes across the full range of pressure ratios.
A.4 AIRFOIL SHOCK
The third setting we consider is a shock-airfoil interaction problem involving a NACA 0012 airfoil immersed in a compressible flow being exposed to a moving normal shock of varying strengths. This case evaluates ShockCast's s ability to model high-speed flows over curved solid boundaries using the immersed boundary method (IBM) implemented in HyBurn. The computational domain is a two- dimensional Cartesian rectangle of size 5 cm 1.25 cmn and the airfoil geometry is imported from an ASCII STL file and scaled by a factor of 0.02, yielding an effective chord length of approximately that the lower-left corner of its axis adis-aligned bounding box lies at (z, y) = (1 cm, 0.5 cm). The airfoil surface is treated as an adiabatic no-slip wall, with IBM guard cells and halo layers used to enforce the boundary conditions and maintain surface refinement during regridding.
2 cm. The geometry is then translated
We initialize normal shock inside the domain at x = 5 mm, with post-shock conditions to the left of the shock and quiescent ambient air (300 K, 1 atm) ahead of it. Across the dataset, we vary the initial shock strength between Mach numbers 1.2 and 2.1, and we vary the angle of attack of the airfoil between - and +8 by rotating the underlying STL geometry. This yields total of 100 cases. of which 90 are used for training and 10 for evaluation.
We employ Adaptive Mesh Refinement (AMR) with a base grid of 64 x 16 cells and up to four levels of refinement, with refinement ratio of 2x between successive levels. This gives an effective finest resolution of approximately 50 um to resolve the shock adequately. Refinement is triggered using shock sensors based on pressure rise and IBM surface proximity, ensuring that both the shock front and the airfoil geometry remain well-resolved throughout the simulation. We use fifth-order WENO reconstruction scheme with Rusanov flux and low-Mach correction to the Riemann solver.
Once the simulation begins, the normal shock propagates from left to right and impinges on the airfoil leading edge. The interaction generates strong compression on the windward side of the airfoil, followed by shock reflections and possible separation depending the angle of attack. The flow downstream of the airfoil develops complex vortical structures whose strength increases with both shock Mach number and angle of attack. We set the CFL numbers for the hyperbolic and parabolic terms to 0.9 and integrate the compressible Navier-Stokes equations using a shree-stage, third-order SSP Runge- Kutta method. All cases simulated until final time of 1.2 x 10-
22

---

Published as a conference paper at ICLR 2026
Initial
Latent Dimension
Model
Key Hyperparameters
Learning Rate
U-Net
2 x 10
64
downssapiing/uusaapling levels
downssapling/upsappligg levels, residual blocks per level
CNO
2 x 10-4
27
J-FNO
1 X 10-3
96
32 modes, 12 layers
Transolver
6 x 10-1
192
layers, 8 attention heads, 8 slices
ConvNeXt-T (See Section of Liu et al. (2022bb)
ConvNeXT
2 10-1
96
MGN
2 x 10-4
128
16 message passing layers
downsampinn/uusaapping levels, mes- sage passing blocks per level downsampling levels, 2 message passing blocks per level
DGN
2 x 10-4
96
DGN-CFL
x 10-5
128
U-Net-3D
2 x 10
32
dowasamplingg/upsaplligg levels
R-FNO-3D
x 10
96
16 modes, 12 layers
ConvNeXt-T (See Section 3 of Liu et al.
ConvNeXT-3D
2 x 10-4
96
(202bb))
Table 2: Hyperparameters for neural solver and neural CFL models.
which is sufficient for the shock to traverse the chord length and for the primary reflection patterns form.
B NEURAL SOLVERS AND NEURAL CFL MODELS
Various neural models for solving PDEs have emerged over the last decade, including Physics-
Informed Neural Networks (Raissi al. 2019 Biswas & Anand 2023
2024
Cho al. 2024
Shah & Anand (024)) and operator learning Lu et al. | Gupta et al. 2021
Li et al. 2022b a
20236 Poli et al.] 2022, Seidman al. 2022 Kovachki al.) 2023). Neural solvers have been
tailored to variety of applications of PDE modeling, including subsurface modeling (Deng et al.,
2022 Wu al. 20226). climate and weather modeling (Bi al.] [202]] Lam et al.
[022]] Pathak
et al.
2022 Price et al., 2023) Nguyen et al.] 2023). airfoil design (Bonnet et al. 2022 Helwig
al. 2024), and density functional theory (Zhang et al. (026). Neural solvers span a diverse array
of architectures, including convolutional models (Li ct al.
2021 Tran et al. 2023 Wen et al.
2022
Helwig et al., 2023 Wen et al. 2023 Bonev et al. 2023 Zhang et al. 2024b; Raonic et al.
2024
transformers Cao 2021 Li et al. 2023a Janny al. [203]] Hao et al. 2024 Alkin et al.
2024
and graph neural networks LLi et al. 2020a Brandstetter et al. 2022 Horie & Mitsume 2022).
Previous works have highlighted the importance of multiscale processing mechanisms in design- ing effective neural solvers Gupta & Brandstetter (2023). This at least in part due to the temporal coarsening approach taken by neural solvers. While differential operators defining PDEs are primar- ily local in nature, their effects become increasingly global as the timestep size is increased beyond that of the classical solver. It is therefore vital for neural solvers incorporate spatial processing mechanisms that enable modeling of phenomena on both local and global scales. Thus, we explore variety of neural solver backbones ShockCast spanning both hierarchical and parallel multi-scall processing mechanisms. Within the parallel framework. we consider both convolution- and attention-based based mechanisms.
U-Net. The U-Net (Ronneberger et al.) (2015) is one of the most prominent examples of hierar- chical mechanism. The U-Net composed of a downsampling path and an upsampling path. To achieve global receptive field, the input features sampled the original solution mesh are first sequentially downsampled using pooling operations or strided convolutions. At each resolution, the downsampled features convolved and point-wa wine activations applied. Following the down- sampling path, an inverse upsampling path is applied, consisting of upsampling operations which either transposed convolutions interpolation, with convolution and non-linearities again ap-
23

---

Published as a conference paper at ICLR 2026
Peak Train Memory (GiB)
Model
Parameters (M)
GigaFLOPS
U-Net
Affine
140.861
94.328
13.465
Euler
140.861
94.328
14.152
MoE
131.129
85.839
21.307
CNO
Affine
45.710
16.028
8.951
Euler
45.704
16.028
9.629
MoE
174.339
55.917
27.783
F-FNO
Affine
15.372
20.901
18.953
Euler
15.374
20.901
20.164
MoE
15.697
20.762
37.242
Transolver
Affine
11.139
191.112
41,895
Euler
10.336
191.110
41.881
MoE
9.911
182.924
62.432
ConvNeXT
27.822
1.736
11.117
MGN
4.719
382.129
39.383
DGN
4.443
124.484
17.445
DGN-CFL
0.795
15.577
4.305
Table 3: Model parameter counts, GigaFLOPs per forward pass, and peak GPU memory usage during training. For U-Net, CNO, F-FNO, Transolver, and ConvNeXT, FLOPs were computed with a batch size of 1 on the coal dust explosion dataset, while GPU memory was computed for a batch size of 32 single A100 GPU. For MGN, DGN, and DGN-CFL, FLOPs were computed with batch size on the airfoil shock dataset, while GPU memory was computed for a batch size on single RTX 6000 Ada GPU.
24

---

Published as a conference paper at ICLR 2026
plied at each resolution. To restore high-frequency details that were lost along the downsampling path, skip connections from respective resolutions along the downsampling path concatenate feature maps each stage of the upsampling path. The U-Net architecture we employ in our framework is the "modern U-Net" architecture from Gupta & Brandstetter 2023), which closely resembles architectures used by diffusion models (Ho et al. 2020).
CNO. The Convolutional Neural Operator (Raonic et al. 2024) adapts the U-Net into the neural
operator framework. Neural operators (Kovachki et al., 2023]] [Seimmnn al., (022) are class of
neural solver which aim to maintain the continuous nature of the underlying PDE solution despite the discretization of training data. In doing so, they enable trained neural solvers to be evaluated on
discretizations differing from the training data. CNO extends these properties to the U-Net architec- ture using anti-aliasing techniques from Karras et al. (2021).
F-FNO. While the previous architectures take a hierarchical approach to multi-scale processing, Fourier Neural Operators (Li et al. 202T) process information multiple scales in parallel using global Fourier convolutions. Convolution kernels are parameterized in the frequency domain such that the complex-valued weights to be learned represent the coefficients of the kernel function in the Fourier basis. Due to the convolution theorem, convolutions in the frequency domain are carried out via point-wise multiplication of frequency modes. Furthermore, because Fourier basis functions have global support, with high frequency modes describing fine details and low frequencies describ- ing the "background" of the function, information on multiple spatial scales processed in parallel. To execute these convolutions more efficiently, FNOs truncate the number of non-zero modes in each kernel to a threshold such that only the lowest frequencies are present in the Fourier expan- sion of the kernel function. Building on this efficiency, Factorized Fourier Neural Operators (Tran) ict al.) 2023) perform convolutions one spatial dimension at a time such that kernels are a function of only spatial dimension. This enables deeper F-FNOs, enhancing the expressive capacity of the architecture.
Transolver. As an alternative to parallel multi-scale processing with Fourier convolutions, atten- tion enables mesh points both distant and local to share information with one another. However, due to its quadratic complexity, adapting Transformers to PDE modeling tasks, where the number of mesh points can be on the order of thousands and above, presents computational challenges (Li al. 22233a Cao 2021 Hao al. (023). Transolver (Wu al. 2024) reduces this complexity by performing attention on a coarsened mesh. The coarsening is achieved by learnable soft pooling operation, where the soft assignments are not entirely based on clustering local points together.
ConyNeXT. For domains discretized on regular Cartesian grids, we employ a ConvNeXT neural CFL backbone (Liu et al. 2022b) to predict the timestep size. A convolutional stem with stride- 4 downsampling followed by a sequence of ConvNeXT stages progressively reduces the spatial resolution while expanding the channel dimension. Each ConvNeXT stage is composed of residual blocks with depthwise convolutions and pointwise MLPs interleaved with layer normalization. After the final stage, global pooling over the spatial dimensions yields a single feature vector, which is normalized and passed through a linear head to produce the scalar timestep prediction.
MeshGraphNets. The MeshGraphNets PPaff al. 2021) (MGN) architecture is GNN model capable of solving PDEs on arbitrary geometries with possibly irregular domains, such as the airfoil shock setting we consider. Each computational cell is treated as node whose features collect the local flow variables together, while edges connect nearby nodes and carry relative geometric information. In our airfoil shock setting, node features additionally include distance to the airfoil, while edge features are pairwise displacements and distances. We construct a eeneerest--nighhee graph with = using the mesh coordinates, and encode node and edge features with separate multilayer perceptrons. The encoded graph is then processed by a stack of message--pssinng blocks that alternate between edge updates and node updates, where edge blocks aggregate information from source and target nodes and update edge features, and node blocks aggregate incoming edge messages via a permutation pormutation-invariant reduction such as sunmation to update node states. Finally, a node decoder maps the processed node features back the predicted flow variables the mesh.
25

---

Published as a conference paper at ICLR 2026
Diffusion Graph Net. The Diffusion Graph Net (DGN) Valencia et al. 2025 combines a multi- scale graph neural network (Lino et al. 2022 Gao Ji 2019) with timestep conditioning to model flow dynamics arbitrary geometries. While the original architecture was used as denoising dif- fusion model (Valencia et al., 2025). we instead utilize it as a timestep-conditioned neural solver. We use the same input node features, edge features, and graph construction process described for MGN. Starting from the input graph, sequence of mesh coarsening operations constructs pro- gressively coarser graphs by clustering nearby nodes and pooling their incident edges. At each resolution, interaction-networf blocks update node and edge features by aggregating messages over local neighborhoods. Information is then propagated back to the finest scale via graph unpooling layers with skip connections, yielding a U-Ner-style encoder-decodee hierarchy in graph space. A final linear decoder maps the latent node features at the finest scale.
DGN Neural CFL. For domains with non-uniform geometries, we reuse the encoder (downsam- pling) half of the DGN backbone as a neural CFL model to predict the timestep size. We use the same graph construction and node/edge features as in the DGN solver. The resulting graph is passed through the multi-scale
pessagg-ppassing hierarchy: at each scale, interaction-network blocks up- date node and edge features, followed by mesh coarsening layers that cluster nearby nodes and pool their incident edges into progressively coarser sequence of graphs. On the final (coarsest) graph, we aggregate node features using global mean pooling over all nodes, and feed the resulting global embedding through a linear decoder to produce the predicted timestep size.
C TRAINING DETAILS
C.I DATASETS
We coarsen the coal dust explosion cases in time by saving every 100 steps, and applied further coarsening to the saved steps by a factor of 5 for overall coarsening factor of 500x relative to the CFD solver. The coarsest AMR level gave a spatial resolution of 104 x 520. We only use the left- most fifth of the domain along the horizontal axis such that the training resolution We train models on the velocity and temperature fields of the gas, as well as the volume fraction describing the percentage of coal comprising each computational cell.
104 x 104.
For the circular blast cases, we save solutions every 100 CFD steps. The coarsest AMR level for the circular blast setting gives spatial resolution of 256 x 256, which we coarsened further the training resolution of 128 x 128 using averaging. We train models on the velocity, temperature and density fields.
The initial shock strength and particle diameters for the Coal Dust Explosion are sampled ran- domly. Half of them are sampled using uniform sampling and the remaining half are sampled using Bridson's 's algorithm (Bridson, 2007) to fill the parameter space. The initial pressure ratios for the Spherical Blast dataset are uniformly spaced. All of the initial shock strength and angle of attacks for the airfoil shock cases are sampled uniformly. In all settings, the evaluation cases are selected so as to be as diverse as possible by first sampling points the parameter space using Bridson's algorithm and then finding the corresponding case closest to the sampled point.
C.2 DATASET STATISTICS
We visualize the distribution of dataset trajectory parameters in Figure 12, In Figuue 13], we visualize timestep sizes, while in Figure 14] we visualize the distribution of the number of steps per solution.
C.3 TRAINING PIPELINE
We implement our training pipeline in PyTorch (Paszke et al.) 2019) using PyTorch Lightnnigg (Fal- con 2019). Depending on model training memory requirements, we train models on between 1-2 80 GiB A100 GPUs or between 1-8 11 GiB RTX 2080 GPUs We optimize all models using the Adam optimizer (Kingma & Ba 2015) with a cosine learning rate scheduler (Loshchilov & Hutter 2017). We train all neural solver models for 400 epochs using batch size of 32, resulting in over
75K training updates for the coal dust explosion dataset and over 50K training updates for both the
circular blast dataset and the airfoil shock dataset. On the coal dust explosion and circular blast
26

---

Published as a conference paper at ICLR 2026
Solution parameters
Coal Dust Explosion
Circular Blast
Airfoil Shock
2.0
2.0
1.8
1.8
Mach
Mach
1.6
1.4
1.4
1.2
1.2
0.000000
0.00005 0.00010 0.00015
20
Diameter
Ratio
Attack Angle
Train
Train
Train
Evaluation
Evaluation
Evaluation
Figure 12: Distribution of solution parameters.
Timestep sizes for min/max timestep solutions
Coal Dust Explosion
x10
x 10 Circular Blast
110
Airfoil Shock
4.8
1.8
1.5
-7
1.2
3.2
www
0.9
2.4
0.000
0.001
0.002
0.003
0.000
0.002
0.004
0.00000
0.00005
0.00010
Diameter: 5.1 10-5,
Mach: 1.20
Ratio: 1.99
Attack Angle: 1.21, Mach: 1.21
Diameter: 2.4 10
Mach: 2.04
Ratio: 49.03
Attack Angle: 6.19, Mach: 2.04
Figure 13: Timestep size function of time for the cases with the least and most number of steps each setting.
datasets, we train neural CFL models for 800 epochs using batch size of 320 using training noise with level 0.01 (Sanchez-Gonzalez et al., (2020). On the coal dust explosion dataset, this results in over 12K training updates, while for the circular blast dataset, this results in over 6K training updates. For the airfoil shock dataset, we train the neural CFL model for 400 epochs using a batch
Timesteps per solution
Coal Dust Explosion
Circular Blast
Airfoil Shock
15.0
20-
30
12.5
15
Count
Coant
26
Count
7.5
5.0
10
2.5
0.0
Timestops
Timesteps
Timesteps
Figure 14: Number of timesteps per solution.
27

---

Published as a conference paper at ICLR 2026
Circular Blast: ShockCast vs Classical Solver Timing Comparison
10
5,423,84x
1,028.77
10a
Time (3)
107
81.36A
10"
1,35s
109
ShockCCost,
ShockCast,
Low Pidelity
High Fidelity
Figure 15:
Timing comparison for ShockCast on CPU and GPU versus low and high fidelity classical solver for circular blast case.
size of 32 with training noise with a level of 0.01, resulting in over 50K training updates. We present initial learning rates and other key hyperparameters used for all models in Table 2
C.4 PARAMETER COUNTS, FLOPS AND PEAK GPU MEMORY
In Table 3 we present parameter counts, forward FL.OPs, and peak training memory for the coal dust explosion setting and the airfoil shock setting. We compute FLOPS with batch size of one us- ing F1 FlopCount erMode from the torch. ut utils..ii........ flop-countez module, and report training memory observed on single A100 GPU using batch size of 32 in the coal dust explosion setting and on a single RTX 6000 Ada GPU with batch size of 4 in the airfoil shock setting.
offset increased computation when using the MoE timestep conditioning strategy by reducing the latent dimension of all models except CNO according to the square root of the number of experts used. We use four experts for all models except for Transolver, where we use two experts due to high memory consumption. For CNO, we ran into stability issues when attempting scale the number of parameters beyond those reported in (Herde al. 2024). As can be seen in Table 3 this led to reduced amount of computation for the CNO variants of ShockCast relative to the other neural solver backbones. However, when using MoE timestep conditioning, we found that we could stably train CNO with 4 experts with no reduction in the embedding dimension.
In experiments,
C.5 SOLUTION RUNTIME ANALYSIS
We compare the runtime of the classical solver ShockCast with the F-FNO neural solver backbone using affine timestep conditioning following the best practices described in McGreivy & Hakim (202). Specifically,
report both GPU and CPU runtimes for ShockCast, and additionally report both high-fidelity and low-fidelity runtimes for the classical solver. The high-fidelity classical solver refers the use of the original solver settings, while low-fidelity uses a setting that trades off degree of accuracy comparable to the neural solver error for increased efficiency. We compare using a circular blast case from the test set with initial pressure ratio 47.575. The low fidelity solver achieved by removing the finest level of the mesh. Times shown in Figure 15) while per-step errors relative to the high-fidelity solution are shown in Figure 16 The high fidelity, low fidelity, and neural solutions are visualized in Figure 17
In Table 4 we present the time to compute solution for ShockCast across all neural solvers back- bones and conditioning methods in both the Coal Dust Explosion and Circular Blast settings. timed on both GPU and CPU-only. We report runtimes for the high fidelity classical solver in the Coal Dust Explosion setting in Table
28

---

Published as a conference paper at ICLR 2026
Circular Blast: ShockCast vs Low Fidelity Classical Solver
0.200
ShockCast
Low Fidelity
0.175
0.150
0.125
Relative
0.100
0.075
0.050
0.025
0.000
0.000
0.001
0.002
0.003
0.004
0.005
Time (s)
Figure 16: Time versus Relative error for ShockCast solution and Low Fidelity classical solver solution for the circular blast setting.
Coal Dust Explosion: Runtime (s)
Circular Blast: Runtime (s)
Model
GPU
CPU
Model
GPU
CPU
CNO
CNO
Affine
2.15 (0.07)
42.51 (0.15)
Affine
1.67 (0.03)
40.58 (0.47)
Fuler
2.30 (0.03)
42.57 (0.08)
Fuler
1.94 (0.01)
41.86 (0.54)
MoE
7.14 (0.01)
108.51 (0.28)
MoE
6.12 (0.10)
117.45 (2.34)
F-FNO
F-FNO
Affine
1.41 (0.03)
58.05 (0.41)
Affine
1.15 (0.01)
56.32 (2.14)
Euler
(0.00)
56.38 (0.04)
Euler
1.18 (0.01)
56.45 (2.23)
MoE
3.75 (0.04)
82.49 (0.53)
MoE
(0.00)
79.91 (2.56)
Transolver
Transolver
Affine
(0.03)
159.82 (1.10)
Affine
2.34 (0.01)
194.56 (4.49)
Euler
2.21 (0.04)
167.54 (6.15)
Euler
2.33 (0.03)
198.53 (5.82)
MoE
3.39 (0.03)
175.19 (0.56)
MoE
3.76 (0.04)
226.63 (4.59)
U-Net
U-Net
Affine
1.66 (0.01)
85.44 (0.69)
Affine
1.61 (0.01)
94.76 (2.43)
Euler
1.66 (0.01)
87.36 (0.33)
Buler
1.64 (0.01)
95.49 (2.16)
MoE
2.78 (0.01)
79.29 (0.67)
MoE 2.51 (0.04)
84.96 (1.06)
Table 4: ShockCast runtime to compute a solution via autoregressive unrolling in both settings on CPU and GPU, presented mean (standard error).
Minimum
Mean
Maximum
15,592
67,441
128,675
Table 5: Classical solver runtime to compute solution on 16 CPU cores in seconds for the Coal Dust Explosion setting.
29

---

Published as a conference paper at ICLR 2026
Circular Blast, Low Fidelity Solver: Density
High Fidelity
Low Fidelity 10-1
Neural Solution
High Fidelity
Low Fidelity
Neural Solution
10-9
10-1
2.651
4.036 x 10-1
4.000 10
4.039 x 10
2,084 10-3
2 2,084 x 10
3.003 10-3
6.115 10
6,112 10
- 6.119 x 10
3.192 10
-- 3.192 10
3.212 10
9.303 10
9.304 10
8.308 10
3.501 10
3,501 10
3.519 10-3
1,530
1,530 10
1.500 x
== 4.117
== 4.117
4.144
1.641 10-3
1.841 10
1.84 x 10
4.373 10
1- 4.373 10
4.385 10
2.045 10
2,04 10
2.045 10-3
4.682 10-3
4.682 10
4.691 10-3
10
10'
Figure 17: Density for circular blast, low fidelity classical solver comparison.
30

---

Published as a conference paper at ICLR 2026
D METRICS
Correlation time proportion. We compute Pearson's correlation coefficient for each field and each timestep with the ground truth data. This requires the predicted fields to be on the same temporal grid as the ground truth data, for which we use linear interpolation. Specifically, for the function f : R + R sampled on the grid {t,} where t S fJ+1, we interpolate to the query point 1* [41, [n]
n-1
(IG))A+I - 17) + f((+1)(( -( 61)
f(t)) =
Σ L (1d)))))(t)))
t,+1 -
where 1l [a.b]] is the characteristic function given by
1l [a.00] (I)) [1 [a
O.w.
The correlation time (Kochkov et al., 2021 Lippe et al. 2023 Alkin et al. (024) for a given field defined as the last time i before the correlation sinks below a threshold, which we take to be 0.9 here. We then average this time across all fields and report it as a proportion of the full simulation time such that a perfect prediction would have correlation time proportion of 1.
Mean flow and Turbulence Kinetic Energy. The mean flow is defined as the flow states averaged over time, while the turbulence kinetic energy is the sum of the variances of the fluctuating part of the velocity field. These quantities are given by
1 //// u(t)dt
2T [/ (m(1) (u(1)- u) (v(1) d).
=
TKE
(5)
We approximate the integrals in Equation (5) Jussing the trapezoidal rule as
L. f(t)dt se
1 Es (f(t)+1) + f(t,)),
2
where D, =t+ tt.
E HYPERSONIC AIRFOIL SHOCK
To evaluate the ability of ShockCast to handle hypersonic (Mach 5) settings, extend the Airfoil Shock setting by varying the initial shock strength from Mach to Mach at a fixed angle of attack of (0. We found that despite having good performance the airfoil shock setting, moderately- sized MGN model (16 message passing layers, embedding dimension 128) did not have sufficient capacity for the hypersonic dynamics despite having large GPU footprint. In particular, we found that the y-velocity field was especially challenging to model and thus, we only model the z-velocity. density, and temperature fields in this setting. Additionally. we employed stronger U-Net model adapted to handle the non-uniform mesh. Specifically, we interpolated the non-uniform domain to fine, uniformly spaced 256 x 1024 grid. The U-Net projected input fields to latent space through a point-wise linear projection applied to the mesh points away from the airfoil surface, and learned embedding vector for mesh points on and inside of the airfoil. The U-Net then processed the data as usual. Following the final layer of the U-Net, we masked mesh points on and inside of the airfoil out of the loss. We visualize rollout predictions for held-out test case in Figures 18 to 20
F MASS CONSERVATION ANALYSIS
We analyze the physical consistency of predictions made by ShockCast by examining mass conser- vation in the circular blast setting. As the boundaries for the blast case do not allow any fluid to enter or exit, the total mass should remain constant across time up to numerical precision. That is, the mass function (t): IR + R defined as
M(t) 1( p(n, t)da
31

---

Published as a conference paper at ICLR 2026
Hypersonic Airfoil Shock (Initial Shock Strength=Mach 5.61): Gas Velocity z-component
True
Predicted
Residual Relative Error 0.114
t 1.971 x 10-6
f 2.166 x 10
1 4,029 x 10-6
8 4.324 x 10-
Relative Ervrr 0.116
6.106 10-
6.471 10-5
Relative Eiror 0.124
t = 9.378 x 10-6
1. 8.597 x 10-8
Relative Error 0.248
1= 1.156 x 10-5
1= 1.175 x 10-5
Relative Error 0.090
1.365 x 10-5
1.386 x 10-5
Relative Firor 0.090
x 10
1.503 x 10-5
Relative Error 0.U85
f= 1.779 x 10-5
f 1.800 x 10-5
Relative Error 0.U82
2,087 x 10-5
f 2.108 x 10-5
Relative Error 0.026
2.211 10
2313 10-5
Relative Error 0.027
2.496 x 50-5
J. 2.519 x 10-5
Relative Errrr 0.027
= 2.806 x 10-5
== 2.725 x 10-3
Relative Eror 0.028
3.013 x 10
1= 3.033 x 10
Relative Error 0.032
3.222 1
3.239 10
Relative 0.032
10
956
1922
1548
99
1747
Figure 18: Gas velocity z-component for hypersonic airfoil shock.
32

---

Published as a conference paper at ICLR 2026
Hypersonic Airfoil Shock (Initial Shock Strength=Mach 5.61): Gas Density
True
Predicted
Residual
t 1.971 x 10
{ 2,166 x 10
Relative Error 0.007
4.029 x 10
f 4.324 x 10""
Relative Error 0.104
6.106 10-6
6.471 10-
Relative Eroo 0.115
9.378 x 10-6
8.597 x 10-1
Relative Error 0.190
1.156 x 10-5
f 1,175 x 10-8
Relative Error
0.001
1.365 10-5
1,386 10-6
Relative Error 0.000
1,573 x 10-5
1,593 x 10-6
Relative Eroo 0.085
1= 1.779 x 10-5
t= 1,800 x 10-8
Relative Error = 0.083
2.087 x 10-5
2.108 10-5
Relative Froor 0.062
2,201 10
x 1
Relative Eiror 0.062
= 2.496 x 10-5
2,519 x 10-6
Relative Error 0.062
2,806 10-5
2,725 x 10-6
Relative Error 0.065
3.0333 5
3.033 5
Relative 0.006
3.222 x 10-
3,239 x 10-8
Relative Error 0.066
Figure 19: Gas density for hypersonic airfoil shock.
33

---

Published as a conference paper at ICLR 2026
Hypersonic Airfoil Shock (Initial Shock Strength=Mach 5.61): Gas Temperature
True
Predicted
Residual
t 1.071 x 10-6
: 2.166 x 10
Relative Error 0.096
t 4,029 x 10-6
8 4.324 x 10
Relative Ervor 0.097
6.106 10-
6,471 10
Relative Eiror 0.102
9.378 x 10-6
1. = 8.597 x 10-8
Relative Error 0.205
1= 1.156 x 10-5
1= 1.175 x 10-5
Relative Error 0.075
1.365 x 10-5
1.386 x 10-5
Relative Firor 0.076
x 10
8 1,503 x 10-5
Relative Error U.U72
=== 1.779 x 10-5
t 1,800 x 10-5
Relative Error 0.070
2.087 x 10-5
t 2.108 x 10-5
Relative Error 0.030
2.20 10
2:13 10-5
Relative Error 0.030
2.496 x 10-5
1. 2.519 x 10-
Relative Froor 0.030
= 2.806 x 10
= 2.725 x 10
Relative Eroor 0.031
3.01. x 10
1= 3.033 x 10
Relative Error 0.033
3.222 1
3.239 10
Relative 0.034
264
1781
3207
1934
1948
Figure 20: Gas temperature for hypersonic airfoil shock.
34

---

Published as a conference paper at ICLR 2026
Relative Mass Deviation from Initial Mass
Minimum Total Deviation
Maximum Total Deviation
Relative Deviation
0.002
Relative Deviation
0.001
0.000
0.000
0.001
0.002
0.003
0.004
0.005
0.000
0.001
0.002
0.003
0.004
0.005
True
Predicted
Figure 21: Relative mass deviation for the two cases with the minimum and maximum prediction total deviation, where the total computed by summing over time.
Long Circular Blast: Relative Error over Time
0.4
Relative Error
0.2
0.0
0.000
0.005
0.010
Time (s)
Figure 22: Visualization of per-step relative rollout error for the long circular blast cases. The error averaged over fields. The solid line is the mean over the test set, and the shaded region two times the standard error.
should remain constant, where Q is the spatial domain and p is the density. To analyze the evolution of the mass function over time, we define the relative mass deviation as
M M())
[M()
which positive for > the mass has changed since the simulation started. As can be seen in Figure 21, ShockCast's predicted total mass stays within 0.2% of M(0).
G LONG CIRCULAR BLAST ANALYSIS
In this section, we analyze performance of ShockCast in the long circular blast setting with a U-Net backbone. The rollout етто over time visualized Figure 22 Although the solution predicted by ShockCast has larger relative еттоr as more steps are taken, the rollout visualizations in Figures 23 show that predictions remain stable and in the distribution of solutions.
H NEURAL CFL ABLATION
In this section, we evaluate and ablate the neural CFL module in ShockCast by replacing it with two alternatives: an oracle model and a mean prediction model. The oracle conditions the neural solver on the ground-truth At for each timestep. As discussed in [Sectinn 3.1 this is not tractable approach, the ground truth At are determined by the full solution and therefore are not known before the solution is computed. However, comparing the performance of ShockCast to the oracle model is of interest, the neural CFL model is trained using the oracle A/ values. Therefore, the oracle represents an upper bound on achievable performance. The mean prediction model is more viable than the oracle and simpler alternative to the neural CFL model. At each timestep, the mean prediction model predicts the mean At computed from the training data, independent the currett flow state ut)..
35

---

Published as a conference paper at ICLR 2026
Long Circular Blast: Velocity z-component
True
Predicted
Residual
True
Predicted
Residual Relative Error 0.419
- 9.848 x 10-5
9,902 in
Rolative Eroor 0.003
5.728 10
t 5,730 x 10-3
Skng S0gg30
Shep
<- 7.116 x 10
7,130 10
Relarive Eroor 0.006
G.407 10
6.401 x 100
Relative Erao 0.22
Step
Step
Skep
1,300 x 10
1.339 inm
Relarive Ernor 0.000
7.042 100
2- 7,009 x 10
Relative Error 0.41
70
Seep 60
Soyp
2,007 x 10
2,008 10
Relarive Frnor 0.015
7.694 10
1- 7.644 x 10
Relative Emro 0.4%8
Stop 500-30
Skep 75
Soyp
2,002 10
2.063 10
Relarive Erior 0.000
8.245 10
2- 8.270 x 10
Relative Eroo 0.61
Slop 30
Skpp XL
Skpp
< 3.267 10
3. 169 10
Relative Eroor 0.508
8.883 10
8,920 x 10
Relative Emoo 0.74
Step 31
Sup 39
Sup
C- 3,950 x 10
1- 3.010 10
Relarive Proor 0.177
9.557 10
9:481 x 10
Relative Emao 0.702
Stip SM
500 37
Stpp 92
== 4.322 30-3
1 4.309 10
Relative Eroor 0.162
1 1.017 10
1.018 x 10-2
Relative Error 0.564
Step 13
34
1917 9
Step
= 5.100 3
5.117 5
Relarive Froor 0.250
1.070 10
1.075 x 10
Relative Fmor 0,00
Skep 105
Seep
11P
-00
107
10
-107
107
107
107
Figure 23: Velocity z-component for long circular blast.
36

---

Published as a conference paper at ICLR 2026
Long Circular Blast: Velocity y-component
True
Predicted
Residual
True
Predicted
Residual Relative Error 0.417
- 9.848 x 10-5
9,902 un
Rolative Eroor 0.003
5,728
10
t 6,730 x 10-3
Skgg
Shep
<- 7.115 x 10
7,130 10
Relarive Erior 0.006
G.417 10"
- 6.401 x 100
Relative Error 0.54
Step S4
Skep 61
Step
1.300 x 10
t- 1.339 inm
Relarive Ernor 0.000
7.042 10"
2- 7.009 x 10
Relative Error 0.40l
Seep
Soyy12
Sey
2,007 10-3
2,008 10
Relarive Ernor 0.015
7.694 10
2- 7.644 x 10
Relative Fmro 01.77
S09 19
Skep 75
Shop
2,002 x 10-6
2,063 10
Relarive Erior 0.03
8.245 x 10
2- 8.270 x 10-
Relative Error 0.609
Step 36
Sepp X0
3:267 100
169 10
Relative Eroor 0.501
8.883 10
8,929 x 10
Relative Erro 0. 0.770
Step 31
Suy 99
Step
3,800 10
3:100 10
Relarive Eroor 0.191
9.007 10
1- 9:481 x 10
Relative Error 0.750
Step SM
Step 92
<= 4.302 30
4.309 10
Relative Error 1816
1.017 10
1.018 x 10-2
Relative Error 0.576
Step 13
30pp24 34
Step 90
-- 1.100 10
5.117 10
Relarive Frror 0.206
7.070 10
1.075 x 10-2
Relative Firo 0.021
S00 S
Seep
100
107
107
107
107
107
10
Figure 24: Velocity y-component for long circular blast.
37

---

Published as a conference paper at ICLR 2026
Long Circular Blast: Density
True
Predicted
Residual
True
Predicted
Residual Relative Error 0.362
- 9.848 x 10-5
9,902 un
Rolative Error 0.002
5.728 10--
6,730 x 10-3
Shep
t- 7.115 10
7,130 10
Relarive Eroor 0.003
G.407 10"
- 6.401 x 10
Relative Error 0.170
Step T
M4
Skep 61
1.330 x 10
1.339 in
Relarive Ernor 0.006
2 7.042 10"
2- 7.009 x 10
Relative Error 0.13
So9p 70
Step
Seeg
-- 2,007 x 10
2,008 10
Relarive Frnor 0.000
7.694 10
1- 7.644 x 10
Relative Faror 0.187
Stop 30
S09 19
Shep 75
2,002 10
2,063 10
Relarive Erior 0.001
8.245 10
2- 8.270 x 10
Relative Erao 0.20
Stop 20
Supp Ki
< 3.267 10
169 10
Relative Eroor 0.215
8.883 10
8,929 x 10
Relative Error 0.25
Step 31
34
Sep
t- 3,900 10
1- 3.100 10
Relarive Erior 0.148
9.557 10
- 9.481 x 10
Relative Emro 0,236
Stip SM
50p.33 33
Step 92
Sapp
== 4.322 3
4.309 10
Relative Eroor 0.163
1.017 10
1.018 x 10-2
Relative Error 018
Step 13
Step
= 1,10 3
== 3.117 5
Relarive Frnor 0.150
1= 1.079 10
1.075 x 10-2
Relative Firo 0,29g
Sa0 S
31
SHep 100
19
107
10'
10'
Figure 25: Density long circular blast.
38

---

Published as a conference paper at ICLR 2026
Long Circular Blast: Temperature
True
Predicted
Residual
True
Predicted
Residual Relative Eruo 0
- 9.848 x 10-5
9,902 un
Rolative Eroor 0.001
5.728 10
5,730 x 10-3
Stap
Siep
<- 7.115 10
7.130 10
Relarive Eroor 0.00G
6.407 100
2- 6.401 x 10
Relative Eraor 0.22
S4
Skep 61
Step
1.300 x 10
1.339 i0
Relarive Ernor 0.00
2 7.042 100
1- 7.009 x 100
Relative Error 0.154
Step
1007
Segg
1- 2,007 10
1- 2,008 10
Relarive Frnor 0.34
2- 7.694 10
7.644 x 10
Relative Earor 0. 0,150
Soy 19
Swa 73
Stop
2,002 10
2.003 10
Relarive Erior 0.100
8.245 10
2- 8.270 x 10
Relative Erao 0 220
Shap X2
Stop
3.267 107-
150 10
Relative Eroor 0.191
8.883 10
8,929 x 10
Relative Emro 0.24
Step 31
2A7 99
Step
t- 3,900 10
1- 3.10 10
Relarive Eroor 0.171
9.057 10
7- 9-481 x 10
Relative Emro 0.24
Step 3M
Stap K4
Step 92
4.322 30
4.309 10
Relative Error 156
1.017 10
1.018 x 10-2
Relative Error 0.24
Step 13
Step
= 5.100 1
== 3.117 3
Relarive Froor 0.206
1.077 10
1.075 x 10-2
Relative Fior 01944
Sepp' 53
31
1AA FOR
Step 105
107
107
Figure 26: Temperature for long circular blast.
39

---

Published as a conference paper at ICLR 2026
Coal Dust Explosion: Correlation Proportion @ 90
Gas Temperature
Volume Fraction
0.8
Correlation Proportion
Correlation Proportion
0.6
0.6
0.4
0.4
0.2
0.2
0.0
0.0
GT ShockCast Mean At
GT At
ShockCast Mean At
Gas Velocity z-component
Gas Velocity y-component
Correlation Proputiin
0.8
Correlation Propurtion
0.3
0.6
0.2
0.4
0.1
0.2
0.0
GT M
ShockCast Mean AM
0.0
GT AM
ShockCast
Mean A
Figure 27: Ablation on the neural CFL model in the coal dust explosion setting.
Circular Blast: Correlation Proportion @ 90
Density
Temperature
1.00
1.00
Correlation Proportion
Correlation Proportion
0.75
0.75
0.50
0.50
0.25
0.25
0.00
0.00
GT At
ShockCast Mean At
GT At
ShockCast Mcan At
Velocity x-component
Velocity y-component
1.00
1.00
Correlation Proportion
Correlation Proportion
0.75
0.75
0.50
0.50
0.25
0.25
0.00
0.00
GT AM
ShockCast Mean At
GT At
ShockCast
Mean At
Figure 28: Ablation on the neural CFL model in the circular blast setting.
40

---

Published as a conference paper at ICLR 2026
Coal Dust Explosion: Correlation Proportion Difference x10 2 ())
Mean Correlation Proportion
Model
Mean Correlation Proportion
Difference
CNO
Affine
60.03 (0.53)
-9.24 (0.11)
Euler
61.04 (0.92)
-12.00 (2.08)
MoE
61.05 (0.85)
-1..58 (2.60)
F-FNO
Affine
61.05 (0.16)
14.770 (0.51)
Euler
60.95 (0.13)
- 15..76 (0.56)
MoE
59.88 (0.94)
17..42 (0.94)
Transolver
Affine
60.55 (0.76)
-1.61 (2.27)
Fuler
59.89 (1.10)
10.47 (0.39)
MoE
60.22 (0.17)
-..88 (0.74)
U-Net
Affine
61.25 (0.12)
9.84 (0.72)
Euler
60.89 (0.26)
9.66 (0.52)
MoE
60.96 (0.29)
10.57 (0.64)
Table 6: Correlation proportion averaged over fields for ShockCast and difference between the ShockCast correlation proportion and the correlation proportion for the neural solver component
using the ground truth At in the circular blast setting.
Coal Dust Explosion: Neural CFL Feature Importance Study
Mean Abvahhte Normalized Hiror
0.35
M((t)
Q()(()
0.30
(gu(()
0.25
Sound Speed
|vel,
[vely]
Speed
0.10
0.0
0.2
0.4
0.6
0.8
1.0
Dropout Probability
Figure 29: Feature importance analysis for the Neural CFL model in the coal dust explosion setting.
We present results for the oracle (GT At) and mean prediction (Mean At) models compared to ShockCast using U-Net backbone with time-conditioned layer norm in Figures 27 and 28 In both settings, we observe that the mean model performs substantially worse than ShockCast. This is likely because the predicted At for the mean prediction model is independent of the current flow state ull), despite the dependence between these variables observed in Figure 13) This results in mean timestep size AL and current flow state u(l) being an dustr--distiiuution input pair for the neural solver. We additionally observe that ShockCast achieves performance close to its upper bound in the oracle model. We provide comprehensive comparisons of ShockCast to the oracle model with various neural solver backbones and conditioning methods in Tables 6lad
I FEATURE IMPORTANCE ANALYSIS
To assess the degree to which each of the input features to the neural CFL model influence predic- tions, we conduct feature importance analysis by progressively corrupting each feature in the coal dust explosion setting. We analyze the effect this perturbation has by measuring the Mean Absolute
41

---

Published as a conference paper at ICLR 2026
Circular Blast: Correlation Proportion Difference x x10 (4)
Mean Correlation Proportion
Model
Mean Correlation Proportion
Difference
CNO
Affine
96.61 (0.15)
0.16 (0.11)
Euler
96.56 (0.08)
0.07 (0.10)
MoE
96.46 (0.02)
0.05 (0.03)
F-FNO
Affine
97.84 (0.14)
0.10 (0.08)
Euler
97.81 (0.05)
0.01 (0.02)
MoE
97.62 (0.03)
0.12 (0.06)
Transolver
Affine
96.92 (0.24)
0.18 (0.09)
Fuler
96.87 (0.17)
0.17 (0.09)
MoE
96.82 (0.10)
0.30 (0.05)
U-Net
Affine
98.34 (0.26)
0.04 (0.04)
Euler
98.09 (0.16)
0.22 (0.09)
MoE
97.82 (0.05)
0.08 (0.02)
Table 7: Correlation proportion averaged over fields for ShockCast and difference between the ShockCast correlation proportion and the correlation proportion for the neural solver component using the ground truth At in the circular blast setting.
Normalized Error, defined as
A
Å aA
where A and Å are the ground truth and predicted timesteps, and as is the standard deviation of the timestep sizes over the training data. We corrupt inputs by applying dropout of increasing levels one field at a time. Results are shown in Figure 29 As can be seen, the flow state u(t) has the largest influence on predictions, followed by the partial derivative with respect to y. This is consistent with the physics of the coal dust explosion setting, in which the shock repeatedly reflects between the upper and lower channel boundaries, generating strong vertical gradients that require small timesteps to resolve as the shock approaches the upper or lower boundary.
J SOLUTION VISUALIZATIONS
In this section, we visualize true and predicted fields for ShockCast on selected solutions from the evaluation datasets. For the coal dust explosion and circular blast setting, ShockCast predictions are made using the F-FNO with Fuler conditioning neural solver backbone, while for the airfoil shock setting, ShockCast predictions are made with the DGN neural solver. The selected solution has shock Mach number of 1.85 in the coal dust explosion setting, max Mach number of 2.68 in the circular blast setting, and shock Mach number of 1.51 in the airfoil shock setting. We visualize TKE fields in Figures 30 to 32 and mean flow fields in Figures 33 to 35 When visualizing the instantaneous fields, subsample the true solution in time to reduce the number of snapshots. For each of the snapshots uI) in the subsampled solution, we find the snapshot from the solution â()) autoregressively unrolled by ShockCast with the closest predicted time f to t. We then visualize the true snapshots u(t) alongside the corresponding closest-in-time predicted snapshots ú(f), well as the residual between each pair (u(t) ú(f). We present these visualizations for each of the fields in Figures 36 to 39 for the coal dust explosion case, Figures 40) to 4 for the circular blast case, and Figures 44 to 47 for the airfoil shock case.
42

---

Published as a conference paper at ICLR 2026
Coal Dust Explosion: TKE Predicted
True
Residual
1000
500
500
1000
Relative Error 0.035
1000
2000
3000
4000
5000
Figure 30: TKE for coal dust explosion.
Circular Blast: TKE Predicted
True
Residual
2000
2000
Relative Error 0.021
10000
20000
30000
40000
50000
G0000
Figure 31: TKE for circular blast.
Airfoil Shock: TKE
True
Predicted
Residual (Rel. Err. 0.080)
10
103
10
10
-10
102
10
Figure 32: TKE for airfoil shock.
43

---

Published as a conference paper at ICLR 2026
Coal Dust Explosion: Mean Flow
True
Predicted
Residual
Gas Velocity z-component
0
-
Relative Error 0.003
200
300
400
Gas Velocity y-component
2.5
0.0
-25
-5.0
Relative Error 0.218
10
0.10
Volume Fraction
0.05
0..00
0.05
0..00
Relative Error = 0.071
0.00
0.05
0.10
0.15
0.20
0.25
0.30
0.35
Gas Temperature
Relative Error = 0.002
375
400
425
450
475
500
525
Figure 33: Mean flow for coal dust explosion.
44

---

Published as a conference paper at ICLR 2026
Circular Blast: Mean Flow Predicted
True
Residual
Velocity x-component
10
Relative Error = 0.025
50
100
150
200
Velocity y-component
10
10
Relative Error 0.027
50
100
150
200
Density
-2
Relative Error 0.021
10
12
14
16
50
Temperatrr
-50
Relative Error 0.024
200
300
400
500
600
Figure 34: Mean flow for circular blast.
45

---

Published as a conference paper at ICLR 2026
Airfoil Shock: Mean Flow
True
Predicted Velocity z-component
Residual Relative Error 0.038
100
200
300
400
100
Velocity y-component
Relative Error 0.105
-50
100
150
200
25
25
Density
Relative Error = 0.017
0.0
0.5
1.0
1.5
2.0
2.5
0.2
0.0
Temperature
Relative Error 0.006
100
200
300
400
10
10
Figure 35: Mean flow for airfoil shock.
46

---

Published as a conference paper at ICLR 2026
Coal Dust Explosion: Gas Velocity z-component
True == 1,625
Predicted
Residual Rolative Error 0.00
True
Predicted
Residual Relative Erroe 0.065
= 1.621 x 10-3
1.586 10--
1,882 10
t 2.009 10
2 2.37 x 10-1
Relative Error 0.003
1,786 10-1
1.782 10
Relative Ernoe 0.004
3,857 10m
3,852 10
Relarive Error (.02)
1,995 10
1,001 10
Relative Erooe 0.006
5.854
10
Relarive Frror 0.01
2.087 10
2.081
= 7,303 10
A= 7.323 x 10
Relaive Eroor 0.022
2.290 x 10
= 2,281 10-3
Relative Error 0.065
8.855 10
8.309 10
Rolatiee Erar 0.08
2.444 10-2
2.437 10-3
Relative Eroor 0,008
1.055
10
Relativ Error 0.088
2.500
2.6-2
Relative 0,082
= 1.287 10
6= 1,232 x 10
Rolaive Error 0.088
2.807 x 10
< 2.796 10
Relative Erroe 0.009
1.406 10
1.431 10
Relative Error 0.006
2.964 10
2.951 10:
Relative Erooe 0.008
10
10
10
Figure 36: Gas velocity z-component for coal dust explosion.
47

---

Published as a conference paper at ICLR 2026
Coal Dust Explosion: Gas Velocity y-component
True == 1,225 10m1
Predicted
Residual Relatve Error 0.05
True
Predicted
Residual Relative Erroe 0.107
= 1.621 x 10-3
1.586 10--
1,882 10
< 2,009 10
2 2.357 x 10
Relative Error 0.25
1,786 10-
< 1.782 1013
Relative Ernoe 0.156
3,857 10ml
3,852 10-1
Relarive Eraor 01261
1,595 10
1,001 10
Relative Erooe 0.215
5.854
3.818 10-4
Relarive Frror 0.197
2.087
2.081
Relative 0,1 126
= 7,363 10
4= 7.323 x 10
Relaive Er 0.263
2,290 10-3
= 2,281 10-3
Relative Error 0,0
8.855 10
8.309 10
Rolative Erao 0.320
2.444 1
<= 2.437 10-3
Relative Eroo 0,173
1,055
10
Relativ Error 0.207
2.500
= 2.622
Relative 0.528
=: 1.287 10
4= 1,232 x 10
Rolaive Error 0.2
2,807 x 10
< 2.796 10
Relative Erroe 0,208
1.406 16
1.431 10
Relarive Error 0.301
2.964 10
2.951 10:
Relative Erooe
0 0,229
107
10
Figure 37: Gas velocity y-component for coal dust explosion.
48

---

Published as a conference paper at ICLR 2026
Coal Dust Explosion: Volume Fraction
True == 1,025 10-1
Predicted
Residual Relative Ervor 0,001
True
Predicted
Residual Relaive Ervor 0.1
= 1,021 x L0-s
1,596 10
1,582 10-1
2,309 10
2.357 Lo-+
Relative Ertur 0.12)
1,786 10-
1,782 10-7
Relative Ervor 0,01
t 5.851 10
= 5.819 x 10-4
Relative Ertor 0.088
2.087
= 2,081 10
Residual Norm 0.206
C- 7,063 10
<- 7.123 x 10
Relative Eror 0.122
t- 2,200 10
t- 2,284 10-2
Residual Norm 0.123
8,855 10-1
8.509 -0
Relative Ertor 0.161
2.464 10-
2,457 10-
Residual Norm 0.068
Residual
( 1,377 10
1.252 x 10
Relative Ercur 0.17L
= 2.907
= 2,790 10
Residual Norm 0.000
<- 1.43 10
<- 1.431 x 3
Relative Enar 0.130
<- 2.964 10
-- 2,901 10
Residual Norm 0.000
10"2
10-1
-00 10F4
10
10-1
-ir 16
Figure 38: Volume fraction for coal dust explosion.
49

---

Published as a conference paper at ICLR 2026
Coal Dust Explosion: Gas Temperature
True == 1,625 1016
Predicted
Residual Relatvv Error 0.08
True
Predicted 1,882 1003
Residual Relative Erroe 0.065
1.621 x 10-3
1.588 10-3
- 2.009 10
2 2.357 x 10-1
Relative Error 0.044
1,786 10-1
1.782 101
Relative Ernor 0.065
3,857 10m
3,852 10
Relarive Error 0.055
1,395 10
1,001 10-
Relative Eroor 0.005
5.854
3.818 10-4
Relarive Frror 0.03
2.087 10
2:081
= 7,303 10
A= 7.323 x 10
Relaiv Eroo 0.000
2.290 10
< 2,281 10-3
Relative Error 0.083
8,855 10
8.309 10
Rolatie Erar 0.001
2.444 10-2
<= 2.437 10-3
Relative Eroor 0.003
1.055
10
Relativ Error 0.088
2.500 10
2.6-2
Relative 0,027
== 1,287 10
4= 1.232 x 10
Rolaive Ereor 0.007
2,807 x 10
< 2.796 10
Relative Erroe 0,00g
1.406 16
1.431 10
Relarive Error 0.006
2.964 10
2.951 10
Relative Erooe 0.004
3xx
4
107
Figure 39: Gas temperature for coal dust explosion.
50

---

Published as a conference paper at ICLR 2026
Circular Blast: Velocity z-component
True t= 1,00 1011
Predicted
Residual Relative Error 0.00
True
Predicted
Residual Relative Error 0.081
1.000 x 10
2.013 10--
2.01 10-3
- 4.089 10
2 1.069 10
Relative Error 0.001
2.947 x 10-1
< 2.968 10-3
Relative Ernor 0.062
6.196
10-1
Relarive Error 0.00
3,164 10
3.165 10-
Relative Eroor 0.086
9-426
Relarive Frror 0.00
x
3-478
< 1,263 10-3
2= 1,23 x 10
Relaive Eroor 0.088
<= 3,779 x 10-3
<= 3.781 10-3
Relative Error 0,136
1.551 10-3
1.381 10
Rolative Eraor 0.001
4,081 -
4.093 10-3
Relative Eroo 0, L39
= 1.705
1.765 10
Relative Error 0.00
4,345 10
4.352
Relative 0,126
= 2.079 10
- 2.061 x 10
Rolatie Ereor 0.011
1,651 x 10
< 4,581 10
Relative Erooe 0,261
2.309 16
2.309
Relarive Error 0.0.7
4.901 10
4.906 10
Relative Erooe 150
16
Figure 40: Velocity z-component for circular blast.
51

---

Published as a conference paper at ICLR 2026
Circular Blast: Velocity y-component
True t= 1,00 101+
Predicted
Residual Relative Error 0.00
True
Predicted 2.01 10-3
Residual Relative Error 0.080
1.000 x 10
2.013 10--
- 4.089 10
2 1.069 10
Relative Error 0.001
2,947 x 10-
< 2.948 1013
Relative Ernor 0.069
6.196 10ml
10-1
Relarive Error 0.00
3.164 10
3.16% 10
Relative Eroor 0.088
9:426
Relarive Frror 0.002
3,471
3-478
Relative 0.135
< 1,203 10-3
A= 1,233 x 10
Relatve Er 0.088
3,779 x 10
= 3.781 10
Relative Error 0.1
1.551 10
1.381 10
Rolative Eraor 0.001
4.084 -0
4.093 10
Relative Eroo 0.111
L7755
1,765 10
Relative Error 0.00
4,345
4.352
Relative 0,1 0.127
= 2.079 10
= 2.061 x 10
Rolaive Ereor 0.011
4,654 x 10
< 4,581 10
Relative Erroe 0,255
2.309 10
2:309
Relarive Error 0.0.8
4.901 10
4.908 10
Relative Erroe 156
10
Figure 41: Velocity y-component for circular blast.
52

---

Published as a conference paper at ICLR 2026
Circular Blast: Density
True x- 1,00 10-1
Predicted 1,006 10
Residual Relative Error 0.001
True
Predicted
Residual Relarive Error 0.060
2,633 x 10
2.6.4 x 10
1,089 10
4,089 10
Relative Ervor 0.001
2,917 10-1
2,948 10-1
Relative Eroo 0.08
9.427 x 10
9.425 10
Relative Ervor 0.002
3,471 x 10-1
3.470 10
Rolative Eroo 0.128
1,063
1,03 10
Relative Eraor 0.00g
3.779 x 10
-- 3,784 10
Relarive Error 0
1,551
1.551 10
Relative Exvor 0.007
4,081 10
4.063 10
Relative Eror 0.117
10
2,079
== 2.08L 10-3
Relative Ervor 0.0L8
4.651 x 10-1
4.584 10
Rolative Eroo 0 0.2
2,899 x 10
== 2,59 10-3
Relative Er 0.011
= 1.901 x 10
4= 4.906 10
Relative Eroor 0.156
10P
10F
10
Figure 42: Density for circular blast.
53

---

Published as a conference paper at ICLR 2026
Circular Blast: Temperature
True t= 1.00 10
Predicted 1.000 x 10
Residual Relative Error 0.000
True
Predicted == 2.11 10
Residual Relative Erroe 0.065
2.013 10--
4.099 10
2- 1.069 10
Relative Error 0.001
2,947 x 10
< 2.968 10
Relative Ernor 0.088
6.196
10
Relarive Error 0.00
3.164 10
3.16% 10
Relative Erooe 0.12
9.426
Relarive Frror 0.00
3.478
= 1.203 10-3
A= 1 2233 x 10
Relaive Er 0.088
<= 3,779 x 10-
<= 8.781 10
Relative Error 065
1.551 10
1.381 10
Rolative Eraor 0.088
4.081 -
4.093 10--
Relative Eroo 082
1.703 1079
1.765 10
Relative Error 0.01
4,345 10
4.352
Relative Eroor 0.10
= 2.099 10
= 2.061 x 10
Rolaive Ereor 0.099
1,64 10
<= 4,581 10
Relative Erooe 0,22
2.309 10
2:309
Relarive Eraor 0.943
4.901 10
4.906 10:
Relative Erooe
0.21
20
10
Figure 43: Temperature for circular blast.
54

---

Published as a conference paper at ICLR 2026
Airfoil Shock: Velocity z-component
True
Predicted
Residual Relative Error 0.106
6,335 10-9
7.404 10-*
1.419 x 10-5
1,470 10-5
Relative Error 0.086
2,138 x 10-3
1- 2.186 x 10
Relative Error 0.072
t= 2,849 x 10-3
== 2.891 x 10-5
Relative Error
0.059
t -.7777 x 10-3
3,807 10-6
Relative Error
0.057
4.636 x 10
4.630 10-5
Relative Frmor 0.059
5,231 10-5
5.252 10-
Relative Frmor 0.057
x 10-5
6,168 10-5
Relative Error 0.0G7
1= 6.744 x 10
1= 6.781 x 10-3
Relative Error 0.058
t= 7.678 x 10-5
t= 7.702 x 10-5
Relative Error 0.055
t- 8:277 x 10-3
L- 8.314 x 10
Relative Error 0.066
9.198 x 10-5
9.234 10-5
Relative Error 0.009
9,814 1
9.853 10
Relative Error 0.071
10-4
1.077 10-4
Relative Eiror
0.070
XM
648
-100
100
0
(0)
700
Figure 44: Velocity z-component for airfoil shock.
55

---

Published as a conference paper at ICLR 2026
Airfoil Shock: Velocity y y-component
True
Predicted
Residual Relative Error 0.201
t 6.035 x 10
f 7.404 x 10
1.419 x 10"2
f 1,470 x 10
Relative Error 0.242
2.138 10-5
2,186 10-8
Relative Eror 0.191
2.819 x 10-5
2,891 x 10-8
Relative Error = 0.
8- 3.777 x 10
f 3.807 x 10-5
Relative Error 0.183
4,636 10-5
4.659 10-6
Relative Error 0.152
5,231 x 10-5
5,252 x 10-6
Relative Erro = 063
1 6. 131 x 10-5
6.168 x 10-8
Relative Error = 0.171
6.744 x 10-5
6.781 x 10-5
Relative Frror 0818
7,678 x 10
x 10
Relative Eiror 0.212
8,277 x 10-5
8,314 x 10-6
Relative Error 0..44
9.198 x 10-5
9.234 x 10-6
Relative Error 0.295
9.844 1
(= 9,853 8
Relative 0.354
1.073 x 10
1,077 x 10-4
Relative Error 0.453
Figure 45: Velocity y-component for airfoil shock.
56

---

Published as a conference paper at ICLR 2026
Airfoil Shock: Density
True
Predicted
Residual Relative Error 0.039
t 6.035 x 10
f 7.404 x 10
1.419 x 10"2
f 1,470 x 10"
Relative Error 0.034
2.138 10-5
2.186 10-8
Relative Eror 0.031
2.819 x 10-5
x 10-8
Relative Error 0.028
3.777 x 10-5
3.807 x 10-5
Relative Error
0.028
4,636 10-5
4.659 10-6
Relative Error 0,029
= 5.231 x 10-5
5,252 x 10-6
Relative Erro 0.028
= 6. 131 x 10-5
6.168 x 10-8
Relative Error = 0.029
6.744 x 10-5
6.781 x 10-5
Relative Frror 0.029
10
x 10
Relative Eiror 0,002
8,277 x 10-5
8,314 x 10-6
Relative Error 0.022
9.198 x 10-5
f 9.234 x 10-6
Relative Error 0.001
9.814 5
(= 9,853 8
Relative 0.032
1.073 x 10-4
1,077 x 10-4
Relative Error 0.034
Figure 46: Density for airfoil shock.
57

---

Published as a conference paper at ICLR 2026
Airfoil Shock: Temperature
True
Predicted
Residual Relative Error 0.017
t 6.035 x 10
f 7,404 x 10
1.419 x 10"2
f 1,470 x 10"5
Relative Error 0.016
2.138 10-5
2.186 10-8
Relative Emor 0.014
2.819 x 10-5
x 10-8
Relative Error 0.012
3.777 x 10-5
3.807 x 10-5
Relative Error
0.013
4,636 10-5
4.659 10-6
Relative Error
0.013
= 5,231 x 10-5
5,252 x 10-6
Relative Erro = 0.012
= 6. 131 x 10-5
6.168 x 10-5
Relative Error = 0.012
6.744 x 10-5
6.781 x 10-5
Relative Frror 0.013
7,678 x 10
x 10
Relative Eiror 0.011
8,277 x 10-5
8,314 x 10-6
Relative Error 0.011
9.198 x 10-5
f 9.234 x 10-6
Relative Error 0.013
9.814 -
(= 9,853 8
Relative 0.013
1.073 x 10-4
1,077 x 10-4
Relative Error 0.013
Figure 47: Temperature for airfoil shock.
58

---

Published as a conference paper at ICLR 2026
Coal Dust Explosion: One-Step Prediction Relative Error x x10 (J)
Gas Velocity x-component
Gas Velocity y-component
Volume Fraction
Gas
Model
Mean
Temperature
CNO
Affine
0.94 (0.01)
9.94 (0.02)
3.03 (0.01)
0.38 (0.00)
3.57 (0.00)
Euler
0.97 (0.01)
10.27 (0.08)
3.15 (0.04)
0.38 (0.00)
3.69 (0.02)
MoE
0.92 (0.00)
9.93 (0.04)
2.86 (0.02)
0.36 (0.00)
3.52 (0.01)
F-FNO
Affine
0.93 (0.00)
10.27 (0.06)
3.04 (0.01)
0.36 (0.00)
3.65 (0.02)
Euler
0.93 (0.00)
10.20 (0.01)
3.01 (0.03)
0.36 (0.00)
3.62 (0.01)
MoE
0.93 (0.00)
10.33 (0.03)
3.05 (0.01)
0.36 (0.00)
3.67 (0.01)
Transolver
Affine
1.22 (0.02)
12.98 (0.19)
2.62 (0.01)
0.44 (0.01)
4.32 (0.05)
Fuler
1.18 (0.02)
12.83 (0.22)
2.60 (0.04)
0.44 (0.01)
4.26 (0.07)
MoE
1.20 (0.01)
12.83 (0.07)
2.53 (0.02)
0.43 (0.00)
4.25 (0.02)
U-Net
Affine
0.92 (0.01)
10.32 (0.06)
2.80 (0.02)
0.35 (0.00)
3.59 (0.02)
Euler
0.91 (0.01)
10.27 (0.05)
2.82 (0.02)
0.35 (0.00)
3.59 (0.02)
MoE
0.93 (0.01)
10.36 (0.07)
2.88 (0.01)
0.35 (0.00)
3.63 (0.02)
Table 8: Relative error for one-step predictions on evaluation split of coal dust explosion cases.
K
EXTENDED RESULTS
In this section, we present the numerical values of average evaluation errors and their corresponding standard errors as mean (standard error). In Tables 8 and [ we present one-step errors for Shock- Cast. We note that the timestep predicted by the neural CFL model will not perfectly match the ground truth timestep such that the prediction from the neural solver model will be for time which differs from the ground truth. To compute the unrolled errors in Tables 10 and [T and correlation time proportions shown in Tables 12 and 13] we linearly interpolate ShockCast predictions in time to be sampled the volume fraction field at later timesteps can be sparse, and so we clamp the norm of the ground truth field in the denominator of the relative error to have minimum value of 1. For the mean flow results, which we show in Tables 14 and [5] and TKE results that we present in Table 16] the target quantities involve integrating the instantaneous fields with respect time, and thus, no interpolation is required.
the same temporal grid as the ground truth data. As can be seen in FFiuure 38
L LIMITATIONS AND FUTURE DIRECTIONS
As previously discussed, neural solvers can benefit from time-adaptive schemes, as varying the timestep size according to the rate of change can lead more balanced one-step objectives across flow states with varying gradient sharpness. Here, we have supervised our neural CFL model using timesteps resulting from coarsening temporal mesh computed using the CFL condition. How- ever, approaches that learn to adapt timestep sizes based on a policy that balances solution accuracy with computational cost, is done by Wu al. (2022a) for spatial remeshing, may lead further improvements. While the settings here contain dynamics comprising some of the most prevalent phenomena in high-speed flows, including shocks, blasts, and a fluid-solid interaction, future works should look to study other phenomena such as detonations and boundary layers. As discussed [Sec tion 3.1. the use of adaptive time-stepping results in more balanced training objectives, which often results in improved generalization due to variance reduction (Duchi & Namkoong) 2019). Neverthe- less, neural solvers in general do not include the same theoretical convergence guarantees enjoyed by classical methods. Future works should look to extend these results from the classical setting to machine learning methods. As discussed in Section 3.4, the MoE iimestep-conditions strategy does not leverage the efficiency advantages of sparse MoE architectures, which apply conditional computation over experts. While this design allows for a simpler implementation without load- balancing losses (Shazeer al. 2017 Fedus et al. 2022), future work should explore incorporating
59

---

Published as a conference paper at ICLR 2026
Circular Blast: One-Step Prediction Relative Error x10 (4)
Velocity x-component
Velocity
Model
Density
Temperature
Mean
y-component
CNO
Affine
2.08 (0.01)
2.08 (0.01)
1.82 (0.01)
1.73 (0.01)
1.93 (0.01)
Fuler
2.05 (0.01)
2.08 (0.02)
1.79 (0.01)
1.69 (0.01)
1.90 (0.01)
MoE
1.75 (0.03)
1.75 (0.04)
1.49 (0.01)
1.42 (0.01)
1.60 (0.01)
R-FNO
Affine
1.87 (0.00)
1.87 (0.00)
1.93 (0.00)
2.05 (0.00)
1.93 (0.00)
Euler
1.85 (0.01)
1.84 (0.01)
1.93 (0.00)
2.04 (0.00)
1.92 (0.01)
MoE
1.87 (0.00)
1.86 (0.00)
2.00 (0.01)
2.14 (0.00)
1.97 (0.00)
Transolver
Affine
1.21 (0.01)
1.21 (0.01)
0.95 (0.01)
0.93 (0.01)
1.07 (0.01)
Fuler
1.27 (0.02)
1.27 (0.02)
1.01 (0.01)
0.99 (0.01)
1.14 (0.01)
MoE
1.28 (0.01)
1.29 (0.01)
1.02 (0.01)
0.99 (0.01)
1.15 (0.01)
U-Net
Affine
1.52 (0.00)
1.51 (0.01)
1.34 (0.00)
1.37 (0.00)
1.44 (0.00)
Euler
1.52 (0.00)
1.52 (0.00)
1.34 (0.00)
1.38 (0.01)
1.44 (0.00)
MoE
1.71 (0.01)
1.70 (0.01)
1.49 (0.01)
1.52 (0.01)
1.61 (0.01)
Table 9: Relative error for one-step predictions on evaluation split of circular blast cases.
Coal Dust Explosion: Unrolled Prediction Relative Error x10-² (L)
Gas Velocity x-component
Gas Velocity y-component
Volume Fraction
Gas
Model
Mean
Temperature
CNO
Affine
3.09 (0.08)
45.30 (1.25)
18.02 (1.00)
1.22 (0.05)
16.91 (0.59)
Euler
3.07 (0.09)
45.71 (1.15)
17.53 (0.34)
1.19 (0.03)
16.88 (0.35)
MoE
3.04 (0.02)
45.53 (0.17)
17.58 (0.70)
1.18 (0.01)
16.83 (0.15)
F-FNO
Affine
2.87 (0.08)
42.70 (0.38)
16.72 (0.09)
1.10 (0.01)
15.85 (0.13)
Fuler
2.80 (0.00)
42.51 (0.66)
16.98 (0.27)
1.09 (0.01)
15.84 (0.20)
MoF
2.89 (0.06)
43.27 (0.73)
17.16 (0.32)
1.13 (0.01)
16.11 (0.14)
Transolver
Affine
3.33 (0.07)
42.59 (1.20)
19.59 (0.21)
1.20 (0.02)
16.68 (0.26)
Euler
3.33 (0.02)
43.95 (0.65)
19.26 (0.70)
1.22 (0.03)
16.94 (0.30)
MoE
3.21 (0.11)
42.61 (0.84)
19.56 (0.60)
1.22 (0.03)
16.65 (0.38)
U-Net
Affine
3.03 (0.05)
44.66 (0.49)
16.26 (0.23)
1.12 (0.01)
16.27 (0.16)
Fuler
2.93 (0.03)
44.82 (0.40)
16.86 (0.25)
1.12 (0.00)
16.43 (0.08)
MoF
3.00 (0.04)
45.02 (0.47)
17.12 (0.10)
1.14 (0.01)
16.57 (0.15)
Table 10: Relative error for unrolled predictions on evaluation split of coal dust explosion cases.
60

---

Published as a conference paper at ICLR 2026
Circular Blast: Unrolled Prediction Relative Error x10 - (L)
Velocity x-component
Velocity
Model
Density
Temperature
Mean
y-component
CNO
Affine
6.93 (0.08)
6.99 (0.12)
7.63 (0.02)
7.96 (0.03)
7.37 (0.05)
Euler
6.98 (0.13)
7.06 (0.12)
7.58 (0.06)
7.92 (0.06)
7.38 (0.08)
MoE
6.74 (0.06)
6.77 (0.05)
7.48 (0.07)
7.85 (0.05)
7.21 (0.01)
F-FNO
Affine
5.89 (0.03)
5.86 (0.01)
5.59 (0.01)
5.89 (0.01)
5.81 (0.01)
Euler
5.70 (0.04)
5.73 (0.06)
5.56 (0.03)
5.87 (0.03)
5.71 (0.04)
MoE
5.91 (0.08)
5.95 (0.10)
5.68 (0.01)
6,00 (0.03)
5.89 (0.04)
Transolver
Affine
7.35 (0.31)
7.35 (0.32)
7.59 (0.05)
8.07 (0.02)
7.59 (0.17)
Euler
7.54 (0.16)
7.31 (0.11)
7.63 (0.09)
8.11 (0.04)
7.65 (0.09)
MoE
7.27 (0.17)
7.23 (0.23)
7.62 (0.10)
8.08 (0.06)
7.55 (0.14)
U-Net
Affine
5.45 (0.07)
5.47 (0.10)
5.16 (0.03)
5.44 (0.04)
5.38 (0.05)
Euler
5.50 (0.01)
5.44 (0.02)
5.24 (0.06)
5.55 (0.06)
5.43 (0.03)
MoF
5.46 (0.02)
5.50 (0.02)
5.30 (0.04)
5.56 (0.03)
5.45 (0.01)
Table 11: Relative error for unrolled predictions on evaluation split of circular blast cases.
Coal Dust Explosion: Correlation Time Proportion x10-² (1)
Gas Velocity a-component
Gas Velocity y-component
Volume Fraction
Gas
Model
Mean
Temperature
CNO
Affine
80.00 (0.00)
19.05 (0.17)
80.04 (1.10)
61.05 (1.06)
60.03 (0.53)
Fuler
80.00 (0.00)
21.40 (2.96)
78.53 (1.59)
64.21 (0.92)
61.04 (0.92)
MoF
80.00 (0.00)
19.48 (0.32)
79.27 (1.07)
65.45 (2.44)
61.05 (0.85)
F-FNO
Affine
80.00 (0.00)
19.76 (0.28)
78.74 (0.32)
65.71 (0.30)
61.05 (0.16)
Euler
80.00 (0.00)
20.72 (1.46)
78.21 (0.23)
64.86 (0.94)
60.95 (0.13)
MoE
80.00 (0.00)
16.17 (2.99)
77.64 (0.68)
65.69 (0.25)
59.88 (0.94)
Transolver
Affine
80.00 (0.00)
21.76 (2.71)
75.25 (0.28)
65.19 (0.48)
60.55 (0.76)
Buler
80.00 (0.00)
19.23 (4.55)
75.80 (0.69)
64.48 (0.48)
59.88 (1.10)
MoE
80.00 (0.00)
19.36 (0.34)
76.46 (0.74)
65.06 (0.05)
60.22 (0.17)
U-Net
Affine
80.00 (0.00)
21.90 (1.68)
79.92 (1.13)
63.17 (0.10)
61.25 (0.12)
Euler
80.00 (0.00)
20.02 (0.15)
79.27 (1.01)
64.27 (0.56)
60.89 (0.26)
MoE
80.00 (0.00)
22.08 (1.50)
78.65 (0.96)
63.12 (0.11)
60.96 (0.29)
Table 12: Correlation time proportion for unrolled predictions on evaluation split of coal dust explo-
sion cases.
61

---

Published as a conference paper at ICLR 2026
Circular Blast:
Correlation Time Proportion x10 (t)
Velocity I-component
Velocity
Model
Density
Temperature
Mean
y-component
CNO
Affine
100.00 (0.00)
100.00 (0.00)
93.38 (0.39)
93.07 (0.23)
96.61 (0.15)
Euler
100.00 (0.00)
100.00 (0.00)
93.35 (0.12)
92.89 (0.22)
96.56 (0.08)
MoE
100.00 (0.00)
100.00 (0.00)
93.04 (0.09)
92.79 (0.06)
96.46 (0.02)
F-FNO
Affine
100.00 (0.00)
100.00 (0.00)
96.07 (0.42)
95.30 (0.15)
97.84 (0.14)
Euler
100.00 (0.00)
100.00 (0.00)
95.77 (0.06)
95.49 (0.14)
97.81 (0.05)
MoE
100.00 (0.00)
100.00 (0.00)
95.52 (0.06)
94.96 (0.05)
97.62 (0.03)
Transolver
Affine
100.00 (0.00)
100.00 (0.00)
94.11 (0.50)
93.58 (0.47)
96.92 (0.24)
Euler
100.00 (0.00)
100.00 (0.00)
94.11 (0.39)
93.37 (0.39)
96.87 (0.17)
MoE
100.00 (0.00)
100.00 (0.00)
94.18 (0.42)
93.09 (0.01)
96.82 (0.10)
U-Net
Affine
100.00 (0.00)
100.00 (0.00)
97. (0.84)
96.13 (0.21)
98.34 (0.26)
Euler
100.00 (0.00)
100.00 (0.00)
96.25 (0.42)
96.10 (0.23)
98.09 (0.16)
MoF
100.00 (0.00)
100.00 (0.00)
95.70 (0.07)
95.57 (0.13)
97.82 (0.05)
Table 13: Correlation time proportion for unrolled predictions on evaluation split of circular blast cases.
Coal Dust Explosion: Mean Flow Relative Error x10-2 2 (L)
Gas Velocity x-component
Gas Velocity y-component
Volume Fraction
Gas
Model
Mean
Temperature
CNO
Affine
1.00 (0.02)
23.03 (0.39)
12.42 (0.17)
0.42 (0.01)
9.22 (0.10)
Fuler
0.97 (0.05)
22.68 (0.86)
12.71 (0.30)
0.39 (0.01)
9.19 (0.27)
MoE
0.91 (0.04)
22.79 (0.25)
12.25 (0.28)
0.37 (0.01)
9.08 (0.10)
F-FNO
Affine
0.82 (0.10)
22.59 (0.12)
11.80 (0.12)
0.35 (0.01)
8.89 (0.08)
Buler
0.79 (0.04)
22.29 (0.25)
11.99 (0.31)
0.32 (0.01)
8.85 (0.13)
MoE
0.80 (0.02)
22.44 (0.13))
12.55 (0.24)
0.34 (0.00)
9.03 (0.10)
Transolver
Affine
1.44 (0.09)
23.07 (0.60)
14.32 (0.17)
0.39 (0.03)
9.81 (0.10)
Euler
1.37 (0.10)
23.47 (0.70)
13.46 (0.13)
0.39 (0.01)
9.67 (0.14)
MoE
1.25 (0.16)
22.78 (0.11)
13.20 (0.15)
0.40 (0.02)
9.41 (0.09)
U-Net
Affine
1.06 (0.07)
22.50 (0.24)
11.17 (0.30)
0.36 (0.02)
8.77 (0.10)
Euler
1.03 (0.05)
23.33 (0.30)
11.74 (0.41)
0.37 (0.02)
9.12 (0.18)
MoE
1.06 (0.02)
22.53 (0.48)
11.62 (0.25)
0.35 (0.01)
8.89 (0.17)
Table 14: Relative error for mean flow evaluation split of coal dust explosion cases.
62

---

Published as a conference paper at ICLR 2026
Circular Blast: Mean Flow Relative Error x x10 (1)
Velocity z-component
Velocity y g.component
Model
Density
Temperature
Mean
CNO
Affine
6,25 (0.12)
6.20 (0.24)
2.35 (0.06)
2.33 (0.07)
4.28 (0.11)
Euler
6.24 (0.23)
6.25 (0.21)
2.30 (0.05)
2.28 (0.06)
4.27 (0.13)
MoE
7.94 (2.10)
7.98 (2.14)
2.30 (0.06)
2.24 (0.05)
5.12 (1.08)
F-FNO
Affine
5.77 (0.22)
5.70 5.70 (0.15)
1.73 (0.01)
1.68 (0.02)
3.72 (0.10)
Euler
5.77 (0.19)
5.81 (0.12)
1.72 (0.02)
1.66 (0.03)
3.74 (0.08)
MoE
5.12 (0.13)
5.12 (0.07)
1.71 (0.01)
1.69 (0.02)
3.41 (0.03)
Transolver
Affine
6.21 (0.13)
6.04 6.04 (0.23)
2.27 (0.01)
2.37 (0.01)
4.22 (0.08)
Euler
11.59 (1.40)
10.82 (1.95)
2.43 (0.01)
2.49 (0.01)
6.83 (0.83)
MoF
6.52 (0.28)
6.30 (0.42)
2.34 (0.05)
2.44 (0.02)
4.40 (0.18)
U-Net
Affine
7.94 (2.63)
7.90 7.90 (2.45)
1.62 (0.05)
1.53 (0.03)
4.75 (1.29)
Euler
7.97 (2.34)
7.92 (2.32)
1.67 (0.05)
1.58 (0.02)
4.79 (1.18)
MoE
12.89 (0.02)
12.95 (0.12)
1.76 (0.03)
1.63 (0.02)
7.31 (0.05)
Table 15: Relative error for mean flow on evaluation split of circular blast cases.
Coal Dust Explosion: TKE
Circular Blast: TKE
Relative Error x x10-2 (4)
Relative Error x x10-2 (L)
Model
TKE
Model
TKE
CNO
CNO
Affine
11.55 (0.23)
Affine
2.60 (0.02)
Euler
10.91 (0.16)
Euler
2.69 (0.07)
MoE
10.94 (0.11)
MoE
2.64 (0.05)
F-FNO
F-FNO
Affine
11.01 (0.36)
Affine
2.29 (0.03)
Euler
10.28 (0.11)
Euler
2.23 (0.06)
MoF
10.93 (0.04)
MoE
2.16 (0.03)
Transolver
Transolver
Affine
12.72 (0.50)
Affine
2.68 (0.08)
Euler
12.44 (0.32)
Euler
2.95 (0.04)
MoE
11.79 (0.88)
MoE
2.77 (0.06)
U-Net
U-Net
Affine
9.85 (0.24)
Affine
2.27 (0.12)
Euler
9.21 (0.12)
Buler
2.30 (0.06)
MoE
8.91 (0.16)
MoE
2.36 (0.03)
Table 16: Relative error for TKE on evaluation splits.
63

---

Published as a conference paper at ICLR 2026
sparsity into MoE-based timestep conditioning. Such approach holds strong potential for scaling time-adaptive neural solvers while balancing computational and memory cost.
M BROADER IMPACTS
Neural PDE solvers have been applied accelerating dynamics simulations across a variety of real-world applications of PDE modeling, including weather and climate forecasting, acrodynamics modeling, and subsurface modeling. As neural solvers often do not include guarantees on gener- alization stability over long time-integration periods, it is vital to perform rigorous validation before relying on predictions in applications. Here, we have explored the potential of neural solvers accelerate modeling of high-speed flows. High-speed flows play an important role in the design of a variety of applications with potential for societal impact, including spacecraft, missiles, and atmospheric reentry vehicles. It is therefore important to closely monitor the development of works along this direction.
64


Reviews:
[{"review_id": "f97a8N0Qmb", "reviewer": "Reviewer_Ssu2", "summary": "The paper proposes ShockCast, a two-phase deep learning framework for adaptive time-stepping in high-speed flow modeling. The first phase employs a Neural CFL model to predict the time-step size ($\\Delta t$) based on the current flow state, while the second phase uses a time-conditioned neural solver to evolve the flow field by the predicted $\\Delta t$. The authors generate two new datasets (spherical blast and coal dust explosion) and demonstrate that ShockCast effectively handles the sharp gradients and varying time scales in supersonic/hypersonic flows.", "questions": "1.How were the initial conditions (e.g., Mach numbers, pressure ratios) for the datasets chosen? Are they representative of real-world scenarios?\n\n2.Can you visualize/analyze which features (e.g., velocity gradients, sound speed) most influence the predicted $\\Delta t$?\n\n3.Small errors in $\\Delta t$ prediction may compound during autoregressive rollout. How does ShockCast handle stability over long simulations?\n\n4.Does ShockCast generalize to unstructured meshes or three-dimensional (3D) flows?\n\n5.The paper lacks an analysis of error accumulation beyond 100 autoregressive steps.\n\n6.The requirement for high-fidelity solver data for training may be prohibitive for certain users.\n\n7.Please specify the computational costs of the proposed methods and the actual speedup achieved compared to numerical simulations.", "limitations": "", "rating": 2, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "First work to address adaptive time-stepping in neural solvers for high-speed flows, filling a critical gap in ML-based CFD. Potential to accelerate simulations in aerodynamics, aerospace, and explosion modeling.", "weaknesses": "1.Evaluated only on two synthetic datasets (spherical blast & coal dust explosion). Training neural CFL + solver may still be expensive compared to classical adaptive methods.\n\n2.Be overly dependent on data.\n\n3.The proposed component is not sufficiently validated.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "GElBmFxqsd", "reviewer": "Reviewer_1474", "summary": "This paper, ShockCast, proposes a novel two-phase deep learning framework for adaptive time-stepping in high-speed flow modeling (e.g., supersonic and hypersonic regimes). High-speed flows exhibit transient sharp gradients (shocks) requiring dynamic adjustment of the timestep size (Δt) using the Courant-Friedrichs-Lewy (CFL) condition, which is computationally expensive for classical solvers.\n\nShockCast addresses this by decomposing the task:\n\nNeural CFL Phase ($\\psi$): A ConvNeXt backbone predicts the optimal, large timestep $\\Delta t$ based on the current flow state. This module is trained to emulate the $\\Delta t$ choices from the classical solver used for data generation, circumventing issues caused by coarse computational meshes.\n\nNeural Solver Phase ($\\phi$): The flow state is evolved by the predicted $\\hat{\\Delta}t$. The authors introduce three novel timestep conditioning strategies for various neural solver backbones (F-FNO, U-Net, CNO, Transolver): Euler Residuals, Mixture of Experts (MoE), and Affine/Spatial-Spectral Conditioning.\n\nThe framework is evaluated on two new supersonic flow datasets—Coal Dust Explosion (multiphase) and Circular Blast (single-phase)—and achieves strong performance in accurately predicting both instantaneous fields and integrated physical quantities (TKE and Mean Flow).", "questions": "Missing Speedup Quantification: The core justification is computational efficiency, yet the paper fails to state the final speedup factor (e.g., $1000\\times$) of ShockCast relative to the original classical solver runtime (which can take tens of thousands of seconds). This critical number must be explicitly provided\n\n\nrade-Off for MoE Complexity: The Mixture of Experts (MoE) component, while showing excellent results, dramatically increases memory consumption (e.g., up to $\\sim 62$ GiB)14. The efficiency/performance trade-off needs deeper analysis, as the simpler Euler Residuals or Affine methods often perform comparably", "limitations": "", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "Improved Training Objective for Transient Dynamics: The core motivation is strong: adaptive time-stepping naturally balances the training objective by inversely scaling $\\Delta t$ according to the rate of change. This more evenly distributes the learning difficulty across states with smooth and sharp gradients (i.e., shocks), a highly pertinent consideration for high-speed flows.\n\n\nPhysically-Informed Neural CFL Model: The Neural CFL phase successfully emulates the true adaptive time mesh, as shown by the close match between predicted and true $\\Delta t$ during autoregressive rollout. Furthermore, incorporating physically-motivated inputs like spatial gradients ($\\nabla u$) and CFL features substantially improves the $\\Delta t$ prediction accuracy for the complex multiphase Coal Dust Explosion scenario\n\n\nNovel and Effective Conditioning Strategies: The introduction of Euler Residuals and Mixture of Experts (MoE) as timestep conditioning strategies is technically insightful. The results demonstrate that these methods are highly competitive, achieving the lowest TKE error for the Circular Blast (F-FNO backbone with MoE/Euler) and best Mean Flow/TKE performance for the Coal Dust Explosion (U-Net backbone with MoE/Euler)", "weaknesses": "1. Missing Quantification of Speedup: The paper's primary motivation is the immense computational cost of classical high-speed flow solvers. However, the results section fails to quantify the final speedup achieved by the full ShockCast pipeline (inference runtime) relative to the original classical solver. Without this figure, the practical utility of the entire framework remains unproven. (Only the classical solver runtime is given in Table 5, min: $\\sim 15$K seconds, mean: $\\sim 67$K seconds)\n\n2. Complexity of MoE Implementation: The MoE conditioning significantly increases peak training memory (e.g., F-FNO: $18.9$ GiB (Affine) vs. $37.2$ GiB (MoE); Transolver: $41.8$ GiB (Affine) vs. $62.4$ GiB (MoE))9. While the complexity is offset by reducing the latent dimension for some models, the substantial jump in memory requirement suggests that the MoE approach is challenging to implement and scale, requiring clarification on the trade-off.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "vwZPJyNtGA", "reviewer": "Reviewer_WzPk", "summary": "The paper proposes ShockCast, a two-phase deep learning framework for modeling high-speed flows with adaptive time-stepping.\n\n- Phase one uses a “Neural CFL” model to predict the timestep ∆t;\n\n- Phase two employs a timestep-conditioned neural solver to evolve the flow field.\n\nThe authors evaluate their framework on two internally generated supersonic datasets (circular blast and coal dust explosion), comparing multiple backbone architectures and conditioning strategies.", "questions": "- If the ∆t distribution is known in training data, why we need the two-phases modeling?\n\n- If the two-phases modeling is usefull, the predicted ∆t from Phase one is inherently impossible for it to perfectly match the true ∆t  (based on a continuous value space) . In this case, what is the significance of designing the Phase one?\n\n- Can you demonstrate that ShockCast preserves key physical invariants or avoids divergence during long rollouts?\n\n- How much computational speedup or advantage (vs. one-phase modeling with time variable inputs) is achieved in practice?\n\n- Have you tested the model on higher Mach (>5) or 3D flow cases to assess scalability?\n\n- Have you tested the model on irregular domain cases to assess scalability?\n\nShould you be able to satisfactorily address the points I've raised above, I will accordingly provide a positive rating.", "limitations": "", "rating": 4, "confidence": 5, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "- Clear motivation: aiming to addresses PDEs under varying time intervals.\n\n- Clear architecture modularity: Two phases modeling.\n\n- Inclusion of physically inspired features (∇u, wave speed, CFL terms) is a positive step toward physics-awareness.\n\n- Comprehensive validations across several backbones and conditioning methods.\n\n- Figures are clear and the manuscript is overall readable.", "weaknesses": "- Conceptual shallowness: The so-called “Neural CFL” model merely regresses ∆t from local features; it does not derive from or guarantee compliance with the true CFL stability condition. There is no theoretical guarantee that predicted timesteps are stable or physically valid.\n\n- Lack of physical consistency: The paper never checks conservation of mass, momentum, or energy, which is essential for high-Mach or compressible flows.\n\n- Experimental validation is weak – Both datasets are self-generated 2D toy problems. Under the same physical conditions, which configuration for each benchmark is the most effective?\n\n- No efficiency or stability analysis – The core motivation of adaptive stepping (e.g., computational savings or others) is never quantified; no wall-clock or rollout stability plots of flow evolution.\n\n- Limited generalization – Only Mach < 3 cases are tested; unclear if the method scales to realistic hypersonic or turbulent regimes.\n\n- Over-claimed novelty – Prior works (e.g., continuous-time neural solvers, time-conditioned FNOs, and physics-aware operator networks) already include similar temporal adaptivity ideas; Also many works for supersonic/hypersonic flow modeling (not the first machine learning framework).\n\n- Method complexity vs. benefit – The two-phase design and multiple conditioning mechanisms add heavy machinery without showing clear improvement over a simpler time-conditioned baseline. Reporting a clear and full coparison (vs. graph-based methods, transformer-based methods, FNO-based methods, and etc. ) on one new table.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "QxEJdpR9AG", "reviewer": "Reviewer_b48r", "summary": "This paper proposes using an ML model to predict time steps for fluid simulations that use adaptive time steps.  It then uses those timesteps in an autoregressive fashion to solve compressible flow problems with a neural solver.", "questions": "- It's probably worth making the connection in the into that adaptive time stepping is not some CFD-specific thing - it is core to ML as well, e.g. \"adaptive learning rates\" with SGD/Adam/etc.\n- Probably also worth noting in the intro that the problem with adaptive time stepping is more extreme than what you describe - if you don't have good timesteps, not only will you miss fine-scale flow phenomena, but your simulation will just crash due to predicting negative pressures etc. - before the solution ever has a chance to diverge - see e.g. Patkar et al., \"Towards positivity preservation for monolithic two-way solid–fluid coupling\" (in particular the issue of positivity preservation that is highlighted - even for relatively low-Mach flows)\n- In 2.3, since you are interested in CFD/engineering, it may be worth noting that methods like PINNs are not convergent (they are, in their vanilla versions at least, one-shot predictors of the solution of a PDE, which can't be convergent).  Recent hybrid solver approaches like that of Kaneda et al., \"A deep conjugate direction method for iteratively solving linear systems\" use a neural network for preconditioning but use a classical solver loop to ensure convergence (that paper happens to be for incompressible flow, but same ideas for compressible flow systems).", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "- The paper makes an interesting observation that when using a neural PDE solver, you can't necessarily rely on the CFL number computed on the fine computational mesh used for the simulation, because the neural model will effectively downsample the mesh inside the network, and so the adaptive time step chosen should be on that coarser scale.\n- The numerical experiments seem convincing to me of the utility of the proposed method.", "weaknesses": "- I think \"Our work represents the first steps towards developing machine learning models for high-speed flows\" is a bit of an inflated statement that could be toned down.\n- Although it's a fair point that a neural solver will have a different internal grid resolution than what is used for the PDE outside the network, presumably one knows what is the architecture they're dealing with and what that coarse grid resolution is.  Why not just compute an adaptive time step - classically - using that coarser scale (which should be known from looking at the network architecture you choose) and eschewing the proposed neural CFL predictor?\n- If the answer to the above is \"you could do either, but neural CFL predictor is faster while still just as accurate,\" can we see a study on that?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)