Paper deep dive
Green-Routed Neural Operators:\\Physics Determines Where the Network Reads
Chenhao Si, Ming Yan
Intelligence
Status: not_run | Model: - | Prompt: - | Confidence: 0%
Entities (0)
Relation Signals (0)
No relation signals yet.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:We identify a mismatch between the physical role of transport fields in many PDEs and their usual role in neural operators: PDEs use them to select read coordinates, whereas neural operators typically treat them only as input values. We address this mismatch with the Green-Routed Neural Operator (GRNO), which uses the governing equation to determine where latent features are sampled. A parameter-free equation adapter evaluates the diagnostic relation and constructs a departure map whose values are the read coordinates. A multiscale encoder-decoder combines centered and routed reads of latent features to learn the complete finite-time update. Across five two- and three-dimensional PDE systems, GRNO achieves the lowest mean final relative $L^2$ error on four under 40-step autoregressive evaluation and remains competitive on Keller-Segel. Fixed-weight route interventions reveal strong dependence on direction and spatial alignment in four systems, with weak dependence in Keller-Segel. In independently trained ablations, GRNO achieves lower mean errors than variants that supply the transport field only as an input feature, substitute a learned displacement for the equation-specified route, or apply the route with a spatial misalignment, across all five systems. It also substantially outperforms directly advecting the physical state and learning the remaining update, indicating that equation-specified read coordinates provide an effective structural prior for long-horizon PDE forecasting.
Tags
Links
- Source: https://arxiv.org/abs/2610.05337v1
- Canonical: https://arxiv.org/abs/2610.05337v1
Trouble viewing inline? Open PDF directly →
Full Text
81,291 characters extracted from source content.
Expand or collapse full text
Green-Routed Neural Operators: Physics Determines Where the Network Reads Chenhao Si Affiliation: School of Data Science Affiliation: The Chinese University of Hong Kong, Shenzhen Affiliation: Shenzhen, China Email: chenhao.si@link.cuhk.edu.cn Ming Yan11 1 Corresponding author. Affiliation: School of Data Science Affiliation: The Chinese University of Hong Kong, Shenzhen Affiliation: Shenzhen, China Email: yanming@cuhk.edu.cn Abstract We identify a mismatch between the physical role of transport fields in many PDEs and their usual role in neural operators: PDEs use them to select read coordinates, whereas neural operators typically treat them only as input values. We address this mismatch with the Green-Routed Neural Operator (GRNO), which uses the governing equation to determine where latent features are sampled. A parameter-free equation adapter evaluates the diagnostic relation and constructs a departure map whose values are the read coordinates. A multiscale encoder–decoder combines centered and routed reads of latent features to learn the complete finite-time update. Across five two- and three-dimensional PDE systems, GRNO achieves the lowest mean final relative L2L^2 error on four under 40-step autoregressive evaluation and remains competitive on Keller–Segel. Fixed-weight route interventions reveal strong dependence on direction and spatial alignment in four systems, with weak dependence in Keller–Segel. In independently trained ablations, GRNO achieves lower mean errors than variants that supply the transport field only as an input feature, substitute a learned displacement for the equation-specified route, or apply the route with a spatial misalignment, across all five systems. It also substantially outperforms directly advecting the physical state and learning the remaining update, indicating that equation-specified read coordinates provide an effective structural prior for long-horizon PDE forecasting. October 4, 2026 1 Introduction Neural operators learn maps between functions and provide fast surrogates for partial differential equations (PDEs) [1, 2, 3]. During autoregressive forecasting, predictions become inputs, allowing errors to accumulate across the rollout [4, 5, 6]. Existing architectures model spatial interactions through kernels, multiscale convolutions, attention, and coordinate-dependent queries; we examine how known transport fields should determine their read coordinates. In many evolution PDEs, the equation already specifies more than which field values are relevant: it specifies the coordinates from which the evolving state should be read. Let q denote the state to be advanced, with μ specifying the physical parameters, geometry, and boundary conditions. The current state may first determine an auxiliary field ψμ[q] _μ[q] through an instantaneous diagnostic relation, after which the equation converts ψμ[q] _μ[q] into a velocity, drift, or phase-space vector μ[q] v_μ[q]. We call ψμ[q] _μ[q] the diagnostic field and μ[q] v_μ[q] the diagnostic carrier. The carrier determines a departure map Δt,μ[q] _ t,μ[q], whose values are the read coordinates. Sampling the previous state at these coordinates gives the transport contribution q∘Δt,μ[q]q _ t,μ[q]. Using μ[q] v_μ[q] as an input channel exposes its values at fixed grid indices, whereas using it to construct the departure map changes the coordinates from which features are read. Endpoint supervision constrains the map qn↦qn+1q^n q^n+1 while leaving the carrier-to-coordinate computation unspecified. Models such as Flowers make sampling locations explicit but learn their displacements from forecast data [7]. Our mismatch lies in the computational role of the known carrier: the PDE uses it to determine read coordinates, while a generic endpoint learner leaves this role implicit. We address this mismatch with the Green-Routed Neural Operator (GRNO). A parameter-free equation adapter computes the carrier and departure map, whose values are the read coordinates used to sample latent features in a multiscale encoder–decoder. Centered and routed reads of latent features jointly determine the complete finite-time update, and the diagnostic field, diagnostic carrier, and departure map are recomputed from every predicted state. Thus, the equation determines where the network reads, while the network learns how the state evolves (Figure 1). Across five systems, GRNO achieves the lowest mean final error on four under 40-step autoregressive evaluation and remains competitive on Keller–Segel. Holding trained weights fixed and varying only the departure displacement reveals strong route dependence in four systems. Independently trained controls compare GRNO with variants that supply the transport field only as an input feature, substitute a learned displacement for the equation-specified route, apply the route with a spatial misalignment, or directly advect the physical state and learn the remaining update. This work identifies a simple mismatch—a PDE can use a known field to choose coordinates while a neural operator receives it only as values—and provides an architecture and matched experiments for correcting and testing it. 2 Related Work Neural-operator interactions. Graph and multipole kernels model nonlocal interactions [8, 9]; Fourier and multiwavelet operators learn spectral interactions [2, 10, 11]. CNO and U-NO use multiscale architectures [12, 13], while transformer operators learn interactions through attention and coordinate queries [14, 15, 16, 17]. Representations and coordinates. G-FNO, SFNO, and ReViT encode spatial symmetries [18, 19, 20], RieszNO represents directional structure [21], and RD-JEPA uses predictive latent pretraining for few-trajectory transfer across related PDEs [22]. Spatial transformers and deformable convolutions learn sampling transformations [23, 24]; Geo-FNO and GINO connect physical domains to latent grids [25, 26]. Flowers and PARADIS learn displacements or latent advection [7, 27]. GRNO constructs read coordinates from the diagnostic carrier. Hybrid numerical models. Hybrid solvers learn discretizations [28, 29], corrections within numerical updates [30], or inference-time corrections with numerical simulation [31]. GRNO uses the physical transport computation to select latent features, while learning the complete finite-time update. Green structure. Prior work learns Green kernels, constructs Green-inspired elliptic operators, or embeds discrete Green superposition [32, 33, 34, 35, 36, 37]. GRNO uses the known Green inverse of a diagnostic relation to construct read coordinates for latent features. 3 Motivation: The PDE determines where information is read Problem setup. Given fully observed training trajectories of a PDE with known coefficients and boundary conditions, we learn a map ℳθM_θ from the current field qn=q(t0+nΔt)q^n=q(t_0+n t) to the next field. Forecasting starts from one exact frame and proceeds recursively: q^0=q0 q^0=q^0 and q^n+1=ℳθ(q^n) q^n+1=M_θ( q^n), with prescribed static inputs omitted from the notation. Models are trained separately for each system. We use Navier–Stokes to examine how known diagnostic carriers should enter this prediction map. 3.1 A concrete example: diagnostic routing in Navier–Stokes Let ωn=ω(tn)ω^n=ω(t_n) denote the vorticity at time tnt_n, with Δt=tn+1−tn t=t_n+1-t_n. On the periodic domain Ω=2 =T^2, forced two-dimensional Navier–Stokes can be written as ∂tω+[ω]⋅∇ω=νΔω−αω+f, _tω+ u[ω]\!·\!∇ω=ν ω-αω+f, (1) Δψ=ω−ω¯,[ω]=∇⟂ψ. ψ=ω- ω, u[ω]=∇ ψ. (2) Here ν is viscosity, α is linear drag, f is a prescribed forcing, and ω¯ ω is the spatial mean. The vorticity ω measures local fluid rotation and is the state being advanced. At each time, solving the Poisson equation in Eq. (2) gives the diagnostic field ψ from the current ω. Its rotated gradient ∇⟂ψ∇ ψ gives the transporting velocity [ω] u[ω], which is the diagnostic carrier. Solving the Poisson equation and applying ∇⟂∇ gives [ω]()=∫Ω(,)(ω()−ω¯),(,)=∇⟂G(,), u[ω]( x)= _ K( x, y) (ω( y)- ω )\,d y, K( x, y)= _ x G( x, y), (3) where G is the mean-zero periodic Green function and K is the associated Biot–Savart kernel. This kernel weights vorticity contributions across the domain to obtain the velocity at each target point [38]. Now fix a target coordinate x at time tn+1t_n+1. Backtracking the velocity field defines a departure map Δtn _ t^n, where Δtn() _ t^n( x) is the earlier coordinate of the material arriving at x. Along this characteristic, ωn+1() ω^n+1( x) =ωn(Δtn())+Δtn(), =ω^n\! ( _ t^n( x) )+C_ t^n( x), (4) Δtn() _ t^n( x) =−Δtn()+(Δt2),n:=[ωn]. = x- t\, u^n( x)+O( t^2), u^n:= u[ω^n]. The pullback ωn∘Δtnω^n _ t^n reads the old vorticity at the departure coordinate, while ΔtnC_ t^n collects viscosity, forcing, drag, and the remaining finite-time effects [39]. For intuition, pure uniform rightward transport by one grid spacing makes each target read its left neighbor. The departure map extends this lookup to spatially varying velocities. The update therefore contains the explicit sequence ωn→Biot–Savartn→backtraceΔtn→sampleωn∘Δtn→+Δtnωn+1.ω^n Biot--Savart u^n backtrace _ t^n sampleω^n _ t^n +\,C_ t^nω^n+1. (5) Appending n u^n as an input channel exposes its values at fixed indices. Eq. (4) uses the same velocity differently: it moves the read coordinates before vorticity is read. Figure 1A illustrates this value-versus-coordinate distinction. Figure 1: Equation-specified routing and GRNO. (A) The current state qnq^n determines a diagnostic field ψ[qn]ψ[q^n] and a diagnostic carrier [qn] v[q^n]. Backtracing the carrier gives the departure map [qn] [q^n], which selects the routed physical sample qn∘[qn]q^n\! \! [q^n]; [qn]C[q^n] is the remainder containing all other finite-time effects. GRNO explicitly realizes the carrier-to-coordinate-to-sample chain by using the carrier to construct read coordinates. (B) The parameter-free equation adapter evaluates the known diagnostic relation, supplies the carrier as an encoder feature, and uses a backward trace to construct the displacement n d^n and departure map n ^n. (C) A local multiscale encoder produces features sℓs_ . At each scale, GRNO combines the centered read sℓ()s_ ( x) with the routed read rℓ()=ℐℓ[sℓ](+ℓ())r_ ( x)=I_ [s_ ]( x+ d_ ( x)). The decoder fuses the two reads, predicts the complete state increment, and recomputes the diagnostic field, diagnostic carrier, and departure map from every predicted state during autoregressive rollout. 3.2 Where the mismatch arises: values versus read coordinates Changing the current state affects both the sampled values and the coordinates used to read them. For (q):=q∘[q]T(q):=q [q], let D[q]eDT[q]e denote the first-order response to a state perturbation e. The chain rule gives D[q]e=e∘[q]⏟state changes at a fixed map+(∇q∘[q])D[q]e⏟the departure map also changes.DT[q]e= e [q]_state changes at a fixed map+ (∇ q [q])D [q]e_the departure map also changes. (6) The second term captures the response created when the carrier moves the read coordinates. In GRNO, this response enters through the departure map [q] [q]. When D[q]D [q] contains a Green inverse, a remote state perturbation can move a departure point and thereby change a local sample. Endpoint supervision leaves the route underdetermined. For an endpoint map P and any alternative ~[q] [q], define ~[q]:=(q)−q∘~[q]⟹(q)=q∘~[q]+~[q]. C[q]:=P(q)-q [q] (q)=q [q]+ C[q]. (7) With an unrestricted remainder, a change in the departure map can be offset by a change in the remainder, producing the same next state. Eq. (7) therefore admits multiple internal decompositions for identical training pairs. The governing equation specifies which departure map GRNO should use. The mismatch is the absence of an explicit link from the known carrier to the read coordinates: q→diagnostic relation[q]→backtrace[q]→samplingq∘[q].q diagnostic relation v[q] backtrace [q] samplingq [q]. (8) Routing is most useful when the equation-specified read coordinates materially affect the finite-time read. When centered and routed reads are similar, their predictive contributions are correspondingly similar. Section 5 measures this dependence across the five systems. 3.3 General form: diagnostically routed PDEs The Navier–Stokes computation extends to systems with different states and diagnostic relations. Let q:Ω→ℝcq: ^c denote the evolving state, and let μ collect physical parameters, geometry, and boundary conditions. The current state first determines a diagnostic field through an instantaneous diagnostic relation, μψμ[q]=μ(q),A_μ _μ[q]=S_μ(q), (9) where μ(q)S_μ(q) forms the source from the current state and μ−1A_μ^-1 solves the diagnostic relation with the prescribed boundary and uniqueness conditions. The equation then converts the diagnostic field into a carrier, μ[q]=μ(q):=μloc(q)+ℬμψμ[q], v_μ[q]=V_μ(q):= v_μ^loc(q)+B_μ _μ[q], (10) where μloc v_μ^loc is an optional known local contribution and ℬμB_μ converts the diagnostic field into a carrier. For Navier–Stokes, μ=ΔA_μ= , μ(ω)=ω−ω¯S_μ(ω)=ω- ω, ℬμ=∇⟂B_μ=∇ , and μloc=0 v_μ^loc=0. The carrier enters the evolution equation as ∂tq+μ[q]⋅∇q=ℛμ(q,ψμ[q]), _tq+ v_μ[q]\!·\!∇ q=R_μ(q, _μ[q]), (11) with ℛμR_μ collecting diffusion, reaction, forcing, and other dynamics not represented by the transport term. Conservative and phase-space forms are given in Appendix C. When the diagnostic inverse has a Green representation, the state-to-carrier map can be written as μ[q] v_μ[q] =μloc(q)+ℬμμ−1μ(q), = v_μ^loc(q)+B_μA_μ^-1S_μ(q), (12) μ[q]() v_μ[q]( x) =μloc(q)()+∫Ωμ(,)μ(q)()d,μ=ℬμ,Gμ. = v_μ^loc(q)( x)+ _ K_μ( x, y)S_μ(q)( y)\,d y, K_μ=B_μ, xG_μ. Here GμG_μ is the Green function of the diagnostic relation. The carrier then determines the short-time departure map Δt,μ[q]()=−Δtμ(q)()+(Δt2). _ t,μ[q]( x)= x- t\,V_μ(q)( x)+O( t^2). (13) Let Δt,μP_ t,μ denote the exact physical map over one step. Relative to the equation-specified departure map, it can be decomposed as Δt,μ(q)=q∘Δt,μ[q]+Δt,μ[q],Δt,μ[q]:=Δt,μ(q)−q∘Δt,μ[q]. P_ t,μ(q)=q _ t,μ[q]+C_ t,μ[q], _ t,μ[q]:=P_ t,μ(q)-q _ t,μ[q]. (14) The first term reads the state at the equation-specified read coordinates; the remainder contains all other finite-time effects. Eq. (14) defines the remainder relative to the equation-specified route; choosing a different route induces a different remainder. We call a PDE diagnostically routed when an instantaneous diagnostic relation determines a carrier that, in turn, determines the read coordinates for the evolving state. Purely local reaction–diffusion equations lie outside this class, as do nonlocal fields that enter only as additive sources. 3.4 How route perturbations accumulate over long rollouts We next examine how a coordinate error changes the sampled state and affects subsequent predictions. Consider a controlled perturbation of Eq. (14): qn+1 q^n+1 =qn∘n[qn]+n[qn], =q^n _n[q^n]+C_n[q^n], (15) q~n+1 q^\,n+1 =q~n∘(n[q~n]+n)+n[q~n]+rn. = q^\,n ( _n[ q^\,n]+ η_n )+C_n[ q^\,n]+r_n. Here n η_n perturbs the read coordinates, while rnr_n perturbs the remainder. Assume locally that composition is bounded by BnB_n, the route is Lipschitz with constant LΞ,nL_ ,n, and the remainder is Lipschitz with constant L,nL_C,n. With en=q~n−qne_n= q^\,n-q^n and Mn=‖∇qn‖L∞M_n=\|∇ q^n\|_L^∞, ‖en+1‖2≤An‖en‖2+Mn‖n‖2+‖rn‖2,An:=Bn+MnLΞ,n+L,n.\|e_n+1\|_2≤ A_n\|e_n\|_2+M_n\| η_n\|_2+\|r_n\|_2, A_n:=B_n+M_nL_ ,n+L_C,n. (16) The term Mn‖n‖2M_n\| η_n\|_2 bounds the sampling error caused by a coordinate displacement; MnM_n measures how sharply the state varies in space. The term MnLΞ,nM_nL_ ,n captures feedback from an existing state error to the next departure map. Iterating the recursion gives ‖eN‖2≤(∏k=0N−1Ak)‖e0‖2+∑n=0N−1(∏k=n+1N−1Ak)(Mn‖n‖2+‖rn‖2).\|e_N\|_2≤ ( _k=0^N-1A_k )\|e_0\|_2+ _n=0^N-1 ( _k=n+1^N-1A_k ) (M_n\| η_n\|_2+\|r_n\|_2 ). (17) A route defect can propagate through later rollout steps because each predicted state determines the next diagnostic solve and departure map. Eq. (17) separates the propagation of existing state error from newly introduced route and update errors. The assumptions, proof, and volume-preserving special case are provided in Appendix A. These observations motivate an architecture that preserves the centered computation while exposing the equation-specified route as a separate latent read. 4 Green-Routed Neural Operator Figure 1B–C illustrates the resulting design: the equation determines the read coordinates, while the neural core learns the complete finite-time update. 4.1 Equation-specified routes For state qnq^n, the parameter-free equation adapter evaluates n=μ(qn)=μloc(qn)+ℬμμ−1μ(qn), v^n=V_μ(q^n)= v_μ^loc(q^n)+B_μA_μ^-1S_μ(q^n), (18) and, for spatial transport, constructs a midpoint departure displacement, n+12 y^n+ 12 =μ(−Δt2n()), = b_μ\! ( x- t2 v^n( x) ), n+12 v^n+ 12 =ℐ[n](n+12), =I[ v^n]( y^n+ 12), (19) n d^n =−Δtn+12, =- t\, v^n+ 12, n ^n =μ(+n). = b_μ( x+ d^n). Here ℐI interpolates grid values at off-grid coordinates, and μ b_μ applies the coordinate boundary rule. The field n+12 v^n+ 12 is the current n v^n sampled at n+12 y^n+ 12. Vlasov–Poisson instead uses a symmetric phase-space backtrace consisting of a half-step trace in position, a full-step trace in velocity, and a second half-step trace in position; the exact construction is given in Appendix D. Periodic indices are wrapped, bounded indices are clamped, and outflow interpolation contributions are masked as specified there. The normalized neural input concatenates the state, carrier, static context, and normalized time step, xμn=cat[q,μ(qn),v,μ(n),cμ,Δt/Δtμ].x_μ^n=cat[N_q,μ(q^n),N_v,μ( v^n),c_μ, t/ t_μ]. (20) Here Δtμ t_μ is the reference frame spacing for system μ, and cμc_μ denotes the prescribed static context. 4.2 Multiscale mixing of centered and routed reads Our four-scale encoder–decoder uses skip connections and local depthwise blocks [40, 41]. Its encoder produces features sℓs_ . The displacement is average-pooled to each scale to sample the same feature tensor, rℓ()=ℐℓ[sℓ](μ,ℓ(+ℓ())),ℓ=0,…,3.r_ ( x)=I_ [s_ ]\! ( b_μ, ( x+ d_ ( x)) ), =0,…,3. (21) At each target grid point, sℓ()s_ ( x) provides the local feature vector and rℓ()r_ ( x) provides the vector interpolated from the departure location. The output grid stays fixed. At the coarsest scale, the decoder receives the two reads directly; at each finer scale, it combines them with the upsampled decoder feature from the next coarser scale. Both reads share one encoder; layer widths and block details are given in Appendix E. 4.3 State increment, rollout, and training The finest decoder feature z0z_0 is mapped by the learned output head HθH_θ to predict the complete state increment q^n+1−q^n q^\,n+1- q^\,n, including transport and the remaining dynamics: q^n+1=q^n+Πμ(σΔq,μHθ(z0)), q^\,n+1= q^\,n+ _μ\! ( _ q,μH_θ(z_0) ), (22) where σΔq,μ _ q,μ is the fixed increment scale computed from training data and Πμ _μ enforces a known integral constraint when applicable. At the next step, Eqs. (18)–(19) are recomputed from q^n+1 q^\,n+1. Training backpropagates through a K=4K=4 rollout using ℒ(θ)=ZK−1∑k=1Kρk−1(ℓstate(k)+λ∇ℓ∇(k)+λspecℓspec(k)),L(θ)=Z_K^-1 _k=1^Kρ^k-1 ( _state^(k)+ _∇ _∇^(k)+ _spec _spec^(k) ), (23) where ZK=∑k=1Kρk−1Z_K= _k=1^Kρ^k-1. Loss definitions and hyperparameters are given in Appendix F. The diagnostic solve costs (NlogN)O(N N) for our spectral adapters; the local neural core and routed reads are linear in the number of grid points, and the equation adapter adds no trainable parameters. 5 Experiments We evaluate long-horizon predictive accuracy and isolate the contribution of equation-specified routing across five PDE systems. 5.1 Benchmarks and protocol We evaluate Kolmogorov flow, dissipative SQG, Vlasov–Poisson, 3D Keller–Segel, and 3D porous convection. These systems span periodic and bounded domains, physical space and phase space, and five distinct diagnostic relations for computing the carrier (Appendix C). The baselines cover four spatial interaction mechanisms. ReViT [20] uses local canonical frames for rotationally equivariant attention. RieszNO [21] combines directional Riesz features with spectral mixing. NIPS-AR [17], our autoregressive adaptation of NIPS, combines linear attention with channel-independent Fourier kernels. Flowers [7] learns multihead displacements for feature resampling, making it the closest coordinate-warp comparison to GRNO’s equation-specified read coordinates. Baseline comparisons use matched splits, normalization, training windows, loss, optimizer updates, and checkpoint selection. Models train through K=4K=4 steps and are tested by free H=40H=40 rollouts from prescribed exact frames, with each prediction reused as the next input. We average errors over 12 prescribed test-trajectory windows within each seed and report mean± deviation over seeds 0,1,2\0,1,2\ for the baseline comparison. Route interventions use seeds 0,1,2\0,1,2\, whereas independently trained ablations use seed 0. Full protocols and dimensional or autoregressive adaptations are specified in Appendices F–G. 5.2 Predictive utility GRNO achieves the lowest mean final error on four of five systems (Table 1), reducing error relative to the strongest baseline by approximately 46%46\% on Kolmogorov flow, 28%28\% on SQG, 47%47\% on Vlasov–Poisson, and 80%80\% on porous convection. The strongest comparator is ReViT on SQG and Flowers on the other three tasks. On Keller–Segel, Flowers, RieszNO, and GRNO obtain errors of 0.320.32, 0.340.34, and 0.370.37 in the table’s units. Figure 2a shows lower errors across substantial portions of these four rollouts, extending beyond the four supervised steps to the 40-step test horizon. On Kolmogorov flow and porous convection, GRNO also has smaller standard deviations across training seeds than Flowers: 0.560.56 versus 5.395.39 and 0.160.16 versus 2.712.71, respectively, in the units of Table 1. Figure 2b shows that GRNO recovers the alternating vorticity regions and filaments in Kolmogorov flow and the concentration patterns in porous convection, with localized discrepancies. Appendix J.1 compares both models’ spatial errors across ten rollout horizons. The following experiments isolate the contribution of equation-specified routing. Table 1: Long-horizon predictive accuracy. Final relative L2L^2 error (×10−2× 10^-2) at H=40H=40. Entries are mean± deviation over three training seeds. Lower is better; bold marks the lowest mean. Model Kolmogorov SQG Vlasov–Poisson Keller–Segel Porous convection ReViT [20] 27.84±1.1827.84± 1.18 3.25±0.033.25± 0.03 4.86±0.904.86± 0.90 1.23±0.021.23± 0.02 19.35±1.2719.35± 1.27 Flowers [7] 14.54±5.3914.54± 5.39 13.12±0.5013.12± 0.50 0.30±0.010.30± 0.01 0.32±0.090.32± 0.09 16.51±2.7116.51± 2.71 NIPS-AR [17] 134.25±12.96134.25± 12.96 51.65±0.0651.65± 0.06 19.33±6.1019.33± 6.10 1.46±0.391.46± 0.39 69.33±32.3969.33± 32.39 RieszNO [21] 39.01±0.4839.01± 0.48 8.07±0.638.07± 0.63 0.49±0.050.49± 0.05 0.34±0.010.34± 0.01 34.84±1.7734.84± 1.77 GRNO (ours) 7.80±0.567.80± 0.56 2.34±0.242.34± 0.24 0.16±0.010.16± 0.01 0.37±0.020.37± 0.02 3.34±0.163.34± 0.16 Figure 2: Long-horizon GRNO accuracy and system-dependent route use. a, Relative L2L^2 error over 40-step free rollouts. GRNO gives the lowest error on four systems and remains competitive on Keller–Segel. b, Reference fields, GRNO predictions and absolute errors at h=40h=40 for Kolmogorov flow and porous convection; pointwise errors remain localized around the recovered structures. c, Ratios of final error under the zero, reverse, and shift interventions relative to the correct route, with network weights fixed. In four systems, intervention errors are 5.0–124 times those obtained with the correct route; in Keller–Segel, the ratios remain near one. d, The error differences appear early and persist throughout the rollout. e, Space–time concentration along the selected porous-convection transect for the reference, Carrier feature only, and GRNO. In this selected window, Carrier feature only develops temporal banding, whereas GRNO follows the advancing front. 5.3 Does the fitted model use the route? The accuracy gains above raise a direct mechanistic question: do the fitted GRNOs actually exploit the equation-specified read coordinates, or could comparable results be achieved while ignoring the route? We answer this by keeping each checkpoint and all other operations unchanged, and modifying only the departure displacement in Eq. (19): zero sets n=0 d^n=0, reverse uses −n- d^n, and shift spatially misaligns the displacement with the current state. Table 2: Fixed-weight route interventions. Final relative L2L^2 error (×10−2× 10^-2) at H=40H=40. Entries are mean± standard deviation across the three seed means. System Correct Zero Reverse Shift Kolmogorov 7.80±0.567.80± 0.56 210.18±0.17210.18± 0.17 196.34±0.21196.34± 0.21 130.67±0.63130.67± 0.63 SQG 2.34±0.242.34± 0.24 11.77±0.3611.77± 0.36 24.72±0.2724.72± 0.27 27.58±0.3127.58± 0.31 Vlasov–Poisson 0.16±0.010.16± 0.01 6.71±0.096.71± 0.09 19.75±0.2119.75± 0.21 20.31±0.0620.31± 0.06 Keller–Segel 0.37±0.020.37± 0.02 0.38±0.020.38± 0.02 0.39±0.010.39± 0.01 0.38±0.010.38± 0.01 Porous convection 3.34±0.163.34± 0.16 41.08±0.9141.08± 0.91 47.76±1.1747.76± 1.17 44.42±1.2144.42± 1.21 Each intervention increases mean final error by at least 16.7×16.7× on Kolmogorov flow, 5.0×5.0× on SQG, 41.0×41.0× on Vlasov–Poisson, and 13.0×13.0× on porous convection, while Keller–Segel changes little (Table 2). Figure 2c–d shows degradation factors and rollout profiles; Appendix J.2 shows the corresponding fields and pointwise errors. The fitted predictors depend strongly on route direction and alignment in four systems, and weakly in Keller–Segel. 5.4 Does equation-specified routing improve prediction? Route reliance in a fixed checkpoint does not by itself establish that equation-specified routing improves the learned predictor. We therefore retrain all variants independently to test whether equation-specified read coordinates reduce error. A 2×22× 2 ablation independently enables or disables carrier features and equation-specified latent routing in the shared backbone. Centered only uses neither component, Carrier feature only and Route only use one component each, and GRNO uses both. This separates the value of supplying the carrier from the value of using it to select read coordinates. Table 3: Independently trained ablations and route controls. Mean final relative L2L^2 error (×10−2× 10^-2) at H=40H=40, for seed 0 and 12 prescribed test windows per system. The first four rows form the 2×22× 2 factorial. Bold marks the lowest mean among evaluated variants. Variant Kolmogorov SQG Vlasov–Poisson Keller–Segel Porous convection Centered only 8.1659 4.2917 0.3158 0.3563 6.5352 Carrier feature only 9.3106 3.7392 0.2849 0.3788 6.2905 Route only 9.0460 3.3361 0.2256 0.3710 4.8331 Full (GRNO) 7.8023 2.3354 0.1635 0.3718 3.3426 Learned route 8.7553 2.8387 0.1782 0.3892 3.4775 Shifted equation-specified route 9.5599 3.3667 0.1913 0.3791 4.0076 State advection + residual 11.3067 6.2341 0.6775 0.5907 9.3214 With carrier features present, adding equation-specified routing reduces mean final error by 16.2%16.2\%, 37.5%37.5\%, 42.6%42.6\%, 1.8%1.8\%, and 46.9%46.9\% on Kolmogorov flow, SQG, Vlasov–Poisson, Keller–Segel, and porous convection, respectively (Table 3). GRNO achieves the lowest mean on four systems, whereas Centered only leads on Keller–Segel. On SQG and Vlasov–Poisson, either Carrier feature only or Route only improves on Centered only, and GRNO improves further. On Kolmogorov flow, the improvement over Centered only emerges when carrier features and routing are combined. Figure 2e compares Carrier feature only with GRNO along the same porous-convection transect. We next test how the read coordinates are constructed. The Learned route control replaces the equation-specified departure displacement with a displacement head trained using the forecasting loss. The Shifted equation-specified route control spatially misaligns the equation-specified departure displacement with the current state. Both retain carrier features and are independently trained using their respective read coordinates. GRNO has lower mean errors than both controls on all five systems. Finally, we compare routing the physical state with routing latent features. The State advection + residual control uses the same equation-specified read coordinates to resample the physical state and then applies a centered network to learn the remaining update, without latent routing. GRNO reduces final error relative to this control by 31.0%31.0\%, 62.5%62.5\%, 75.9%75.9\%, 37.1%37.1\%, and 64.1%64.1\% on Kolmogorov flow, SQG, Vlasov–Poisson, Keller–Segel, and porous convection, respectively. State advection + residual has the highest mean error among the evaluated variants on all five systems. These ablations indicate that the gain depends on using the equation-specified route for latent feature reads while learning the complete finite-time update. 6 Discussion GRNO uses the diagnostic carrier both as an input feature and to determine where latent features are read. The diagnostic relation determines the diagnostic field, from which the carrier and departure map are computed. The values of the departure map are the read coordinates used by the neural operator. The equation therefore specifies the read coordinates, while the network learns the complete finite-time update. This interface applies across two-dimensional flow, phase-space dynamics, and three-dimensional transport. Fixed-weight route interventions show that four fitted predictors depend strongly on route direction and alignment. In the independently trained controls, GRNO has lower mean error than Learned route and Shifted equation-specified route on all five systems. It also has lower mean error than State advection + residual on every benchmark. These results indicate that performance depends on both supplying the carrier and using it to determine latent read coordinates. Dependence on routing varies by system. Keller–Segel shows little sensitivity to fixed-weight route interventions. In the ablations, GRNO improves only slightly on Carrier feature only, and Centered only has the lowest error. In the other four systems, the trained models are sensitive to route interventions, and GRNO has the lowest mean error in the ablations. These results support equation-specified latent routing when the route materially affects finite-time prediction. AI use statement The authors used generative AI tools to help polish the text, check the mathematical arguments, give feedback on the experimental design and controls, implement and debug the model code, and write the plotting code based on the model code from their trained models and the requirements they gave. The authors conducted all numerical experiments and reviewed the AI-assisted code, text, and figures. All reported numerical results were obtained from these experiments. The authors take responsibility for the final content of this work. Ethics statement This work uses synthetic numerical PDE data and does not involve human subjects, personal data, or sensitive datasets. Reproducibility Statement All equations, datasets, adapters, boundary rules, model configurations, losses, training budgets, and evaluation metrics are specified in the appendices. The code will be shared after the paper is accepted. References [1] L. Lu, P. Jin, G. Pang, Z. Zhang, and G. E. Karniadakis (2021) Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators. Nature Machine Intelligence 3 (3), p. 218–229. External Links: Document Cited by: §1. [2] Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar (2021) Fourier neural operator for parametric partial differential equations. In International Conference on Learning Representations, External Links: Link Cited by: §1, §2. [3] N. Kovachki, Z. Li, B. Liu, K. Azizzadenesheli, K. Bhattacharya, A. Stuart, and A. Anandkumar (2023) Neural operator: learning maps between function spaces with applications to PDEs. Journal of Machine Learning Research 24 (89), p. 1–97. External Links: Link Cited by: §1. [4] J. Brandstetter, D. Worrall, and M. Welling (2022) Message passing neural PDE solvers. In International Conference on Learning Representations, External Links: Link Cited by: §1. [5] P. Lippe, B. Veeling, P. Perdikaris, R. Turner, and J. Brandstetter (2023) PDE-Refiner: achieving accurate long rollouts with neural PDE solvers. In Advances in Neural Information Processing Systems, Vol. 36. External Links: Link Cited by: §1. [6] B. List, L. Chen, K. Bali, and N. Thuerey (2025) Differentiability in unrolled training of neural physics simulators on transient dynamics. Computer Methods in Applied Mechanics and Engineering 433, p. 117441. External Links: Document Cited by: §1. [7] T. Muser, A. Spitzer, M. Lassas, M. V. de Hoop, and I. Dokmanić (2026) Flowers: a warp drive for neural PDE solvers. In Forty-third International Conference on Machine Learning, External Links: Link Cited by: §1, §2, §5.1, Table 1. [8] Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar (2020) Neural operator: graph kernel network for partial differential equations. External Links: 2003.03485, Link Cited by: §2. [9] Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, A. Stuart, K. Bhattacharya, and A. Anandkumar (2020) Multipole graph neural operator for parametric partial differential equations. In Advances in Neural Information Processing Systems, Vol. 33, p. 6755–6766. External Links: Link Cited by: §2. [10] J. Guibas, M. Mardani, Z. Li, A. Tao, A. Anandkumar, and B. Catanzaro (2022) Efficient token mixing for transformers via adaptive Fourier neural operators. In International Conference on Learning Representations, External Links: Link Cited by: §2. [11] G. Gupta, X. Xiao, and P. Bogdan (2021) Multiwavelet-based operator learning for differential equations. In Advances in Neural Information Processing Systems, Vol. 34. External Links: Link Cited by: §2. [12] B. Raonić, R. Molinaro, T. De Ryck, T. Rohner, F. Bartolucci, R. Alaifari, S. Mishra, and E. de Bézenac (2023) Convolutional neural operators for robust and accurate learning of PDEs. In Advances in Neural Information Processing Systems, Vol. 36. External Links: Link Cited by: §2. [13] M. A. Rahman, Z. E. Ross, and K. Azizzadenesheli (2023) U-NO: U-shaped neural operators. Transactions on Machine Learning Research. External Links: Link Cited by: §2. [14] S. Cao (2021) Choose a transformer: Fourier or Galerkin. In Advances in Neural Information Processing Systems, Vol. 34. External Links: Link Cited by: §2. [15] H. Wu, H. Luo, H. Wang, J. Wang, and M. Long (2024) Transolver: a fast transformer solver for PDEs on general geometries. In Proceedings of the 41st International Conference on Machine Learning, Proceedings of Machine Learning Research, Vol. 235, p. 53681–53705. External Links: Link Cited by: §2. [16] S. Wang, J. H. Seidman, S. Sankaran, H. Wang, G. J. Pappas, and P. Perdikaris (2025) CViT: continuous vision transformer for operator learning. In The Thirteenth International Conference on Learning Representations, External Links: Link Cited by: §2. [17] N. Liu and Y. Yu (2025) Neural interpretable PDEs: harmonizing Fourier insights with attention for scalable and interpretable physics discovery. In Proceedings of the 42nd International Conference on Machine Learning, Proceedings of Machine Learning Research, Vol. 267, p. 40085–40099. External Links: Link Cited by: §2, §5.1, Table 1. [18] J. Helwig, X. Zhang, C. Fu, J. Kurtin, S. Wojtowytsch, and S. Ji (2023) Group equivariant Fourier neural operators for partial differential equations. In Proceedings of the 40th International Conference on Machine Learning, Proceedings of Machine Learning Research, Vol. 202, p. 12907–12930. External Links: Link Cited by: §2. [19] B. Bonev, T. Kurth, C. Hundt, J. Pathak, M. Baust, K. Kashinath, and A. Anandkumar (2023) Spherical Fourier neural operators: learning stable dynamics on the sphere. In Proceedings of the 40th International Conference on Machine Learning, Proceedings of Machine Learning Research, Vol. 202, p. 2806–2823. External Links: Link Cited by: §2. [20] H. Wei, B. List, and N. Thuerey (2026) ReViT: rotational-equivariant vision transformers for neural PDE solvers. In Forty-Third International Conference on Machine Learning, Note: Oral presentation External Links: Link Cited by: §2, §5.1, Table 1. [21] S. Liu, X. Yang, and Y. Chen (2026) Riesz neural operator for solving partial differential equations. In The Fourteenth International Conference on Learning Representations, External Links: Link Cited by: §2, §5.1, Table 1. [22] C. Si and M. Yan (2026) RD-JEPA: predictive latent pretraining for few-trajectory transfer across Reaction-Diffusion equations. External Links: 2609.29403, Link Cited by: §2. [23] M. Jaderberg, K. Simonyan, A. Zisserman, and K. Kavukcuoglu (2015) Spatial transformer networks. In Advances in Neural Information Processing Systems, Vol. 28. External Links: Link Cited by: §2. [24] J. Dai, H. Qi, Y. Xiong, Y. Li, G. Zhang, H. Hu, and Y. Wei (2017) Deformable convolutional networks. In Proceedings of the IEEE International Conference on Computer Vision, p. 764–773. External Links: Link Cited by: §2. [25] Z. Li, D. Z. Huang, B. Liu, and A. Anandkumar (2023) Fourier neural operator with learned deformations for PDEs on general geometries. Journal of Machine Learning Research 24 (388), p. 1–26. External Links: Link Cited by: §2. [26] Z. Li, N. Kovachki, C. Choy, B. Li, J. Kossaifi, S. Otta, M. A. Nabian, M. Stadler, C. Hundt, K. Azizzadenesheli, and A. Anandkumar (2023) Geometry-informed neural operator for large-scale 3D PDEs. In Advances in Neural Information Processing Systems, Vol. 36. External Links: Link Cited by: §2. [27] C. A. Pereira, S. Gaudreault, V. Dallerit, C. Subich, S. Panday, S. Wei, S. Zhang, S. Rout, E. Haber, R. J. Spiteri, et al. (2026) Learning to advect: a neural semi-Lagrangian architecture for weather forecasting. arXiv preprint arXiv:2601.21151. Cited by: §2. [28] Y. Bar-Sinai, S. Hoyer, J. Hickey, and M. P. Brenner (2019) Learning data-driven discretizations for partial differential equations. Proceedings of the National Academy of Sciences 116 (31), p. 15344–15349. External Links: Document Cited by: §2. [29] D. Kochkov, J. A. Smith, A. Alieva, Q. Wang, M. P. Brenner, and S. Hoyer (2021) Machine learning–accelerated computational fluid dynamics. Proceedings of the National Academy of Sciences 118 (21), p. e2101784118. External Links: Document Cited by: §2. [30] K. Um, R. Brand, Y. Fei, P. Holl, and N. Thuerey (2020) Solver-in-the-loop: learning from differentiable physics to interact with iterative PDE-solvers. In Advances in Neural Information Processing Systems, Vol. 33. External Links: Link Cited by: §2. [31] Z. Fan, Y. Sun, S. Yang, and Y. Lu (2026) Physics-informed inference time scaling for solving high-dimensional partial differential equations. In The Fourteenth International Conference on Learning Representations, External Links: Link Cited by: §2. [32] N. Boullé, C. J. Earls, and A. Townsend (2022) Data-driven discovery of green’s functions with human-understandable deep learning. Scientific Reports 12, p. 4824. External Links: Document Cited by: §2. [33] N. Boullé, S. Kim, T. Shi, and A. Townsend (2022) Learning green’s functions associated with time-dependent partial differential equations. Journal of Machine Learning Research 23 (218), p. 1–34. External Links: Link Cited by: §2. [34] H. Melchers and J. Prins (2026) Neural green’s operators for parametric partial differential equations. Computer Methods in Applied Mechanics and Engineering 455, p. 118893. Cited by: §2. [35] H. Wang, J. Li, A. Dwivedi, K. Hara, and T. Wu (2024) BENO: boundary-embedded neural operators for elliptic PDEs. In The Twelfth International Conference on Learning Representations, External Links: Link Cited by: §2. [36] H. Tang, H. Chen, and C. Li (2026) Generalized spherical neural operators: green’s function formulation. In The Fourteenth International Conference on Learning Representations, External Links: Link Cited by: §2. [37] Y. Tan, Q. Yao, and Y. Wang (2026) DGNet: discrete green networks for data-efficient learning of spatiotemporal PDEs. In The Fourteenth International Conference on Learning Representations, External Links: 2603.01762, Link Cited by: §2. [38] A. J. Majda and A. L. Bertozzi (2001) Vorticity and incompressible flow. Cambridge Texts in Applied Mathematics, Cambridge University Press. External Links: Document Cited by: §3.1. [39] A. Staniforth and J. Côté (1991) Semi-lagrangian integration schemes for atmospheric models—a review. Monthly Weather Review 119 (9), p. 2206–2223. External Links: Document Cited by: §3.1. [40] O. Ronneberger, P. Fischer, and T. Brox (2015) U-Net: convolutional networks for biomedical image segmentation. In Medical Image Computing and Computer-Assisted Intervention – MICCAI 2015, Lecture Notes in Computer Science, Vol. 9351, p. 234–241. External Links: Document Cited by: §4.2. [41] Z. Liu, H. Mao, C. Wu, C. Feichtenhofer, T. Darrell, and S. Xie (2022) A ConvNet for the 2020s. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, p. 11976–11986. External Links: Document Cited by: §4.2. [42] S. M. Cox and P. C. Matthews (2002) Exponential time differencing for stiff systems. Journal of Computational Physics 176 (2), p. 430–455. External Links: Document Cited by: §C.2. [43] G. Strang (1968) On the construction and comparison of difference schemes. SIAM Journal on Numerical Analysis 5 (3), p. 506–517. External Links: Document Cited by: §C.2. [44] C. Z. Cheng and G. Knorr (1976) The integration of the Vlasov equation in configuration space. Journal of Computational Physics 22 (3), p. 330–351. External Links: Document Cited by: §C.2. [45] I. Loshchilov and F. Hutter (2019) Decoupled weight decay regularization. In International Conference on Learning Representations, External Links: Link Cited by: §F.1. Appendix A Proofs for the route-sensitivity bounds Assume that, in a neighborhood of the reference trajectory, ‖f∘(n[p]+)‖2 \|f ( _n[p]+ η)\|_2 ≤Bn‖f‖2, ≤ B_n\|f\|_2, (24) ‖n[p]−n[q]‖2 \| _n[p]- _n[q]\|_2 ≤LΞ,n‖p−q‖2, ≤ L_ ,n\|p-q\|_2, (25) ‖n[p]−n[q]‖2 \|C_n[p]-C_n[q]\|_2 ≤L,n‖p−q‖2. ≤ L_C,n\|p-q\|_2. (26) Let en=q~n−qne_n= q^\,n-q^n and MnM_n bound ‖∇qn‖L∞\|∇ q^n\|_L^∞ along the segment joining the two departure maps. Add and subtract qn∘(n[q~n]+n)q^n ( _n[ q^\,n]+ η_n) in Eq. (15). The first resulting term is bounded by Bn‖en‖2B_n\|e_n\|_2; the mean-value theorem and the second inequality bound the coordinate change by Mn(LΞ,n‖en‖2+‖n‖2)M_n(L_ ,n\|e_n\|_2+\| η_n\|_2); the final inequality bounds the remainder difference. Summing the three bounds gives the one-step recursion in Eq. (16). Repeated substitution proves Eq. (17). For a measure-preserving continuum characteristic, change of variables gives Bn=1B_n=1; bounded-domain handling and discrete interpolation are covered by retaining BnB_n. Appendix B Supplementary guide The supplement is organized around the claims made in the main paper. Appendix A proves the route-sensitivity bounds. Appendix C instantiates the diagnostically routed PDE class and specifies every dataset. Appendix D gives the parameter-free equation adapter, traces, and boundary rules. Appendix E records the shared trainable architecture, and Appendix F fixes the training and evaluation contract. Appendix G documents baseline provenance and adaptations. Appendices H–I define the reported metrics and describe both the fixed-weight interventions and independently trained controls. Appendix J provides forecasting and route-intervention visualizations. Table 4: Claim–evidence map. This table separates evidence for predictive utility from evidence for the proposed routing mechanism. Question Control Evidence Does the complete model predict well? Independently trained models share data, loss, update budget, checkpoint rule, and evaluator. Main Table 1 Does the trained model use the route? Keep weights and carrier features fixed; change only the departure displacement. Figure 2 and Appendix I Does equation-specified routing improve prediction after training? Independently train models with different carrier features, coordinate constructions, and routing targets. Main Table 3 and Appendix I Appendix C PDE instantiations and datasets C.1 Variants of the common template Eq. (9)–(11) use an advective notation for readability. When the evolved quantity is a density, the corresponding conservative form is ∂tq+∇⋅(μ[q]q)=ℛμ(q,ψμ[q]). _tq+∇\!·\! ( v_μ[q]q )=R_μ(q, _μ[q]). (27) For a phase-space state f(,,t)f( x, v,t), the same idea applies on the product domain: ∂tf+x[f]⋅∇xf+v[f]⋅∇vf=ℛ(f), _tf+ a_x[f]\!·\! _xf+ a_v[f]\!·\! _vf=R(f), (28) where the diagnostic solve may determine either component of the phase-space carrier (x,v)( a_x, a_v); Table 5 gives the exact instantiation for every task. Table 5: Diagnostic carriers used in the experiments. Each row lists the evolved state, diagnostic relation, diagnostic carrier that determines the read coordinates, and boundary rule. System Evolved state q Diagnostic relation Diagnostic carrier Boundary rule Kolmogorov ω Δψ=ω−ω¯ ψ=ω- ω =∇⟂ψ u=∇ ψ periodic x,yx,y SQG θ ψ=−(−Δ)−1/2θψ=-(- )^-1/2θ =∇⟂ψ u=∇ ψ periodic x,yx,y Vlasov–Poisson f(x,v)f(x,v) ρ=∫fvρ= f\,dv, ∂xE=1−ρ _xE=1-ρ (ax,av)=(v,−E)(a_x,a_v)=(v,-E) periodic x; zero outflow in v Keller–Segel ρ (−Δ+λ2)c=ρ−ρ¯(- +λ^2)c=ρ- ρ =χ∇c v=χ∇ c periodic x,y,zx,y,z Porous convection C =−∇p+Cz u=-∇ p+C e_z, ∇⋅=0∇\!· u=0 Darcy velocity u periodic x,yx,y; bounded z C.2 Equations and coefficients Kolmogorov flow. On [0,2π)2[0,2π)^2, the reference trajectories solve ∂tω+⋅∇ω=νΔω−αω−4cos(4x2),Δψ=ω−ω¯,=∇⟂ψ, _tω+ u·∇ω=ν ω-αω-4 (4x_2), ψ=ω- ω, u=∇ ψ, (29) with ν=10−2ν=10^-2, α=0.1α=0.1, and frame spacing Δt=0.1 t=0.1. The pseudo-spectral reference solver uses 2/32/3 de-aliasing, ETDRK2 [42], and 20 internal steps per stored frame. Initial vorticity is obtained by Fourier filtering independent standard Gaussian grid values to ‖∞≤8\| k\|_∞≤ 8, removing the zero mode, and normalizing the root-mean-square amplitude to one. The field is advanced for 500 stored-frame intervals before collection. Dissipative surface quasi-geostrophy. On the same periodic domain, ∂tθ+⋅∇θ=κΔθ−γ(−Δ)−1θ+f,ψ=−(−Δ)−1/2θ,=∇⟂ψ, _tθ+ u·∇θ=κ θ-γ(- )^-1θ+f, ψ=-(- )^-1/2θ, u=∇ ψ, (30) with κ=0.008κ=0.008, γ=0.15γ=0.15, and Δt=0.1 t=0.1. Each trajectory has a frozen random forcing supported on 3≤||≤43≤| k|≤ 4 with RMS amplitude 0.100.10. The de-aliased ETDRK2 solver uses 40 internal steps per stored frame. Initial conditions use the same Gaussian filtering procedure as Kolmogorov flow, with ‖∞≤8\| k\|_∞≤ 8, root-mean-square amplitude 0.50.5, and a 200-frame warm-up. Vlasov–Poisson. The phase-space domain is x∈[0,4π)x∈[0,4π) and v∈[−6,6)v∈[-6,6), with ∂tf+v∂xf−E∂vf=0,ρ(x)=∫−66f(x,v)v,∂xE=1−ρ, _tf+v\, _xf-E\, _vf=0, ρ(x)= _-6^6f(x,v)\,dv, _xE=1-ρ, (31) and Δt=0.1 t=0.1. A spectral Strang solver [43, 44] applies a half x translation, a full v translation, and a second half x translation, using five internal steps per stored frame. The reference generator implements these translations by FFT on the truncated velocity interval; GRNO instead assigns zero to routed reads that leave [−6,6)[-6,6). The sampled initial distributions have negligible support at that boundary. Each trajectory independently selects a Landau or two-stream initial condition with probability 1/21/2. Writing Mσ(v)=exp[−v2/(2σ2)]/(2πσ)M_σ(v)= [-v^2/(2σ^2)]/( 2πσ), the unnormalized fields are fL(x,v) f_L(x,v) =Mσ(v)[1+αcos(x/2+ϕ)+14αcos(x−2ϕ)], =M_σ(v) [1+α (x/2+φ)+ 14α (x-2φ) ], (32) fT(x,v) f_T(x,v) =12[Mσ(v−v0)+Mσ(v+v0)][1+αcos(x/2+ϕ)+15αcos(3x/2−ϕ)]. = 12 [M_σ(v-v_0)+M_σ(v+v_0) ] [1+α (x/2+φ)+ 15α (3x/2-φ) ]. (33) Parameters are sampled independently and uniformly: ϕ∈[0,2π)φ∈[0,2π); σ∈[0.8,1.2]σ∈[0.8,1.2] and α∈[0.05,0.18]α∈[0.05,0.18] for Landau; v0∈[1.6,2.4]v_0∈[1.6,2.4], σ∈[0.35,0.60]σ∈[0.35,0.60], and α∈[0.03,0.10]α∈[0.03,0.10] for two-stream. Each field is divided by its discrete mean density, (ΔxΔv/Lx)∑i,jf(xi,vj)( x v/L_x) _i,jf(x_i,v_j), so that the spatially averaged density is one. Screened Keller–Segel. On [0,2π)3[0,2π)^3, ∂tρ=DΔρ−∇⋅(ρ),=χ∇c,(−Δ+λ2)c=ρ−ρ¯, _tρ=D ρ-∇\!·(ρ v), v=χ∇ c, (- +λ^2)c=ρ- ρ, (34) with D=0.10D=0.10, χ=0.30χ=0.30, λ=1λ=1, and Δt=0.25 t=0.25. The 2/32/3-de-aliased Fourier ETD2 solver uses eight internal steps per stored frame. Initial densities have the form ρ=exp(0.2ξ)ρ= (0.2ξ), where ξ is a unit-variance random field supported on 1≤||≤31≤| k|≤ 3; each sample is rescaled to mean one. Buoyancy-driven porous convection. On the unit cube, with periodic horizontal coordinates and bounded vertical coordinate z, ∂tC+∇⋅(C)=RaD−1ΔC,=−∇p+Cz,∇⋅=0, _tC+∇\!·( uC)=Ra_D^-1 C, u=-∇ p+C e_z, ∇\!· u=0, (35) with RaD=1000Ra_D=1000 and Δt=0.1 t=0.1. We impose C=1C=1 at the top, ∂zC=0 _zC=0 at the bottom, and impermeability uz=0u_z=0 at both vertical boundaries. Pressure is fixed to zero mean. The finite-volume reference solver combines MUSCL advection, SSP/RK2 time stepping, an exact mixed-boundary diffusion semigroup, and Strang splitting; its internal step is at most 0.020.02 and also satisfies an advective CFL bound of 0.350.35. Trajectories are generated on 1283128^3 and volume-averaged to 64364^3. The vertical coordinate increases downward from the top boundary, z=0z=0. Before the warm-up, concentration is initialized as C0(x,y,z)=clip[0,1][b(z)+0.08b(z)(1−b(z))ξ(x,y)],b(z)=erfc(z2RaD−1t0).C_0(x,y,z)=clip_[0,1] [b(z)+0.08b(z)(1-b(z))ξ(x,y) ], b(z)=erfc\! ( z2 Ra_D^-1t_0 ). (36) Here t0t_0 is sampled uniformly from [0.9,1.1][0.9,1.1]. The horizontal field ξ is a sum of sine and cosine modes with integer wavevectors 1≤||≤81≤| k|≤ 8, using one representative of each pair ,−\ k,- k\. Its coefficients are independent standard Gaussian variables, and the field is normalized to unit root-mean-square amplitude. The initialized field is advanced for 0.20.2 before storage. C.3 Dataset contract Table 6: Dataset contract. Shapes list trajectories × stored frames; every frame also contains the grid shown in the second column. Splits are fixed before training. System Grid Train Validation Test Stored Δt t Kolmogorov 1282128^2 64×5164× 51 8×1018× 101 12×10112× 101 0.100.10 SQG 1282128^2 64×5164× 51 8×1018× 101 12×10112× 101 0.100.10 Vlasov–Poisson 128x×128v128_x× 128_v 64×10164× 101 8×1218× 121 12×16112× 161 0.100.10 Keller–Segel 64364^3 64×5164× 51 8×1018× 101 12×10112× 101 0.250.25 Porous convection 64364^3 64×10164× 101 8×1018× 101 12×10112× 101 0.100.10 Each training run records a trajectory/start-time manifest; a result is included in the matched table only after its manifest hash agrees with the task reference. Manifests log array shapes, coefficients, split membership, solver version, and file hashes; any run with a mismatched contract is excluded. Appendix D Equation adapter, traces, and boundaries The adapter is deterministic, has no trainable parameters, and is recomputed from the current predicted state at every rollout step. D.1 Diagnostic solves Periodic active scalars. For Kolmogorov flow and SQG, the zero Fourier mode is fixed to zero and ψ^()=−||−rq^(),^()=(ik2,−ik1)ψ^(), ψ( k)=-| k|^-r q( k), v( k)=(ik_2,-ik_1) ψ( k), (37) with r=2r=2 for vorticity and r=1r=1 for SQG. This is the discrete Green multiplier for the respective periodic diagnostic relation. Vlasov–Poisson. The adapter integrates f over velocity, removes the neutral Fourier mode, and solves Eq. (31) spectrally for E(x)E(x). It returns the phase-space carrier (v,−E(x))(v,-E(x)) rather than treating E as a second evolved state. Keller–Segel. The screened inverse is diagonal in Fourier space: c^()=ρ^()||2+λ2(≠0),=χ∇c. c( k)= ρ( k)| k|^2+λ^2 ( k≠ 0), v=χ∇ c. (38) The mean source is removed as in Eq. (34). Porous convection. The adapter first forms the buoyancy on vertical cell faces: interior faces use the average of their two adjacent concentrations, and the top and bottom face values are zero. The discrete divergence of this face buoyancy is the pressure right-hand side. The pressure equation uses the finite-volume Laplacian, diagonalized by Fourier transforms in the periodic horizontal directions and a cosine transform vertically; the zero mode is set to zero. Face velocities are computed from this buoyancy and the corresponding discrete pressure gradient, with zero normal velocity at the top and bottom. Averaging adjacent face velocities gives the cell-centered carrier supplied to the network. D.2 Backward traces For spatial transport, GRNO uses the midpoint rule from Eq. (19): the adapter first samples the carrier at −12Δt() x- 12 t\, v( x) and then places the final read at −Δt1/2() x- t\, v_1/2( x). Vlasov–Poisson uses the symmetric phase-space trace matching the order of the characteristics: x1/2 x_1/2 =x−12Δtv, =x- 12 t\,v, (39) E1/2 E_1/2 =ℐ[E](x1/2), =I[E](x_1/2), vd v_d =v+ΔtE1/2, =v+ t\,E_1/2, xd x_d =x1/2−12Δtvd. =x_1/2- 12 t\,v_d. The routed read is taken at (xd,vd)(x_d,v_d). The sign follows x˙=v x=v and v˙=−E v=-E in Eq. (31) while integrating backward from the target point. D.3 Interpolation and boundary rules Table 7: Boundary behavior of routed interpolation. The rule is applied identically at every latent scale. Coordinate type Coordinate handling Sample outside domain Periodic wrap modulo domain length wrapped value Bounded wall clamp to the nearest grid index value at the nearest boundary grid point Phase-space outflow mask each out-of-range index zero contribution Bilinear interpolation is used in two dimensions and trilinear interpolation in three dimensions. At scale ℓ , a departure displacement dℓ,jd_ ,j shifts the fractional grid index by dℓ,j/Δxℓ,jd_ ,j/ x_ ,j, where Δxℓ,j=Lj/Nℓ,j x_ ,j=L_j/N_ ,j. The full-resolution displacement is average-pooled to each scale and remains in physical units. Each neighboring index in the interpolation stencil is wrapped, clamped, or assigned zero interpolation weight according to Table 7. Appendix E Shared GRNO architecture The trainable core uses a common 2D architecture and its direct 3D counterpart. Within a dimension, the parameter names, tensor shapes, stage depths, and channel widths are identical across PDEs. Only the equation adapter, boundary metadata, context channel, and deterministic output constraint vary by task. Table 8: Dimension-specific GRNO configuration. “Axial” in 3D denotes the average of three depthwise kernels, one per spatial axis. Quantity 2D / phase space 3D Input channels 5 6 Base width 48 24 Stage widths (48,96,144,192)(48,96,144,192) (24,48,72,96)(24,48,72,96) Encoder depths (2,2,3,4)(2,2,3,4) (2,2,3,4)(2,2,3,4) Decoder blocks per scale 1 1 Local depthwise kernel 5×55× 5 axial 55-tap Number of scales 4 4 Trainable parameters 1,534,9451,534,945 391,561391,561 Trainable adapter parameters 0 0 Inputs and normalization. Each task advances one scalar channel. The neural input has shape B×5×Ny×NxB× 5× N_y× N_x in two dimensions and B×6×Nz×Ny×NxB× 6× N_z× N_y× N_x in three dimensions. For Vlasov–Poisson the two grid axes are (v,x)(v,x). The channel order is state, carrier components in (x,y)(x,y) or (x,y,z)(x,y,z) order, one context channel, and Δt/Δtμ t/ t_μ. Using fixed scales computed from training data, q,μ(q)=q/sq,μ,2D and phase space,(q−mq,μ)/sq,μ,3D,[v,μ()]j=vj/sv,j,μ.N_q,μ(q)= casesq/s_q,μ,&2D and phase space,\\ (q-m_q,μ)/s_q,μ,&3D, cases [N_v,μ( v)]_j=v_j/s_v,j,μ. (40) Input scale denominators are floored at 10−710^-7. For Kolmogorov flow and SQG, state, velocity, forcing and increment scales are read from the training normalization file; forcing is divided by its scale, and the context channel is the normalized forcing pattern. For Vlasov–Poisson, state and increment scales are root-mean-square values estimated from up to 512 training frame pairs using sampling seed 17. Its carrier scales are vmaxv_ and the root-mean-square electric field, and its context is zero. For the 3D tasks, state mean, state scale and increment scale are read from the saved training statistics. Each carrier component is scaled by its root-mean-square value over 24 evenly selected training states. Porous convection uses cμ(z)=2z/Lz−1c_μ(z)=2z/L_z-1 as context; Keller–Segel uses zero. Multiscale fusion. The grids have side lengths (128,64,32,16)(128,64,32,16) in 2D and (64,32,16,8)(64,32,16,8) in 3D. Every encoder scale produces one feature tensor sℓs_ and its routed copy rℓr_ . The coarsest decoder input is [s3,r3][s_3,r_3]; at the remaining scales, zℓ=Decℓ(Wℓ[Upℓ(zℓ+1),sℓ,rℓ]),ℓ=2,1,0.z_ =Dec_ \! (W_ [Up_ (z_ +1),s_ ,r_ ] ), =2,1,0. (41) Brackets denote channel concatenation, and WℓW_ is a learned pointwise projection. Downsampling uses average pooling and a pointwise projection. Upsampling uses a pointwise projection followed by nearest-neighbor interpolation. A local residual block applies u⟼u+W2GELU(W1LN(DW(u))).u u+W_2GELU (W_1LN(DW(u)) ). (42) Here DWDW is a depthwise spatial convolution, with the three axial outputs averaged in 3D. Layer normalization acts across channels independently at each grid point, with learned scale and shift and numerical floor 10−610^-6. The pointwise maps expand channels by a factor of two and project them back. The increment head HθH_θ is a learned pointwise projection from the finest decoder feature z0z_0 to the normalized state increment and is initialized to zero. For tasks with an exactly conserved spatial integral, the predicted increment is projected to zero mean. Porous convection uses the identity projection because mass enters through its Dirichlet top boundary. These projections are fixed and contain no learned coefficients. Appendix F Training and evaluation contract F.1 Optimization and exposure Table 9: Matched training and evaluation protocol. Four 3D micro-batches are accumulated to obtain an effective batch size of four. Setting 2D / phase space 3D Epochs 100 100 Windows per epoch 2,048 512 Micro-batch / accumulation 4/14/1 1/41/4 Effective batch 4 4 Optimizer updates per epoch 512 128 Training rollout K 4 4 Validation / test horizon 10/4010/40 10/4010/40 Optimizer AdamW AdamW Initial learning rate 3×10−43× 10^-4 3×10−43× 10^-4 Weight decay 10−410^-4 10−410^-4 Schedule cosine to 1.5×10−51.5× 10^-5 cosine to 1.5×10−51.5× 10^-5 Gradient clipping 1.0 1.0 Precision bfloat16 AMP bfloat16 AMP Seeds: baseline comparison 0,1,2\0,1,2\ 0,1,2\0,1,2\ Seeds: fixed-weight route interventions 0,1,2\0,1,2\ 0,1,2\0,1,2\ Seed: independently trained ablations 00 00 For a fixed task, all models receive the same window manifest. An epoch is therefore defined by data exposure and optimizer updates, not by how many micro-batches fit in memory. Three-dimensional models use gradient accumulation only to meet the A100 memory budget. Training uses AdamW [45]; the learning-rate scheduler advances once per epoch. For validation, all admissible (trajectory, start-frame) pairs are ordered by trajectory and then start frame. The evaluator selects 64 evenly spaced pairs in 2D and 32 in 3D, or all pairs if fewer are available. Checkpoint selection minimizes relative L2L^2 averaged over these windows and all ten validation horizons. F.2 Loss Let eb(k)=q^bn+k−qbn+ke_b^(k)= q_b^\,n+k-q_b^n+k for micro-batch sample b. Let B be the micro-batch size, N the number of grid points, and sq,μs_q,μ the fixed state scale defined in Appendix E. All reported tasks have one physical output channel. The state loss is ℓstate(k)=1BN∑b=1B∑|eb(k)()/sq,μ|2. _state^(k)= 1BN _b=1^B _ i |e_b^(k)( i)/s_q,μ |^2. (43) Define ⟨a⟩=(BN)−1∑b,ab() a =(BN)^-1 _b, ia_b( i). The gradient loss is ℓ∇(k)=∑j⟨|δje(k)|2⟩max∑j⟨|δjqn+k|2⟩,ε. _∇^(k)= _j | _je^(k)|^2 \ _j | _jq^n+k|^2 , \. (44) The difference δj _j uses grid spacing Lj/NjL_j/N_j: centered differences with wraparound on periodic axes, centered interior differences and one-sided endpoint differences otherwise. Let ℱF_P transform only the periodic axes. The spectral loss averages per-sample ratios: ℓspec(k)=1B∑b=1B‖|ℱq^bn+k|−|ℱqbn+k|‖22max‖ℱqbn+k‖22,ε. _spec^(k)= 1B _b=1^B \||F_P q_b^\,n+k|-|F_Pq_b^n+k| \|_2^2 \\|F_Pq_b^n+k\|_2^2, \. (45) The spectral norm sums over every transformed coefficient and any remaining grid axis. The gradient and spectral terms use physical field values. The implementation uses ε=10−8 =10^-8, λ∇=0.05 _∇=0.05, λspec=0.02 _spec=0.02, ρ=0.8ρ=0.8, and ZK=∑k=1Kρk−1Z_K= _k=1^Kρ^k-1 in Eq. (23). In 3D, gradients from four micro-batches are accumulated for each optimizer update. Gradients pass through all four rollout steps, including the recomputed diagnostic solves, traces, and interpolations. F.3 Free autoregressive evaluation For test trajectory i∈0,…,11i∈\0,…,11\, let sis_i be the prescribed start frame. Evaluation uses q^i0=qisi,q^ih+1=ℳθ(q^ih),h=0,…,39, q_i^0=q_i^s_i, q_i^h+1=M_θ( q_i^h), h=0,…,39, (46) and compares q^ih q_i^h with qisi+hq_i^s_i+h. The three 2D systems use si=0s_i=0. Porous convection uses si=20(imod4)s_i=20(i 4). Keller–Segel uses the integer start indices (0,5,10,16,21,27,32,38,43,49,54,60)(0,5,10,16,21,27,32,38,43,49,54,60), obtained by evenly spacing twelve indices over the admissible interval [0,60][0,60]. Errors are averaged over the twelve windows within each seed. Tables 1 and 2 report the mean and sample standard deviation across the three seed means. Table 3 reports seed 0. Appendix G Baseline provenance and matched adaptations Table 10: Baseline roles and adaptations. For dimensional or autoregressive adaptations, we retain each published model’s interaction mechanism and modify only the input and output interfaces required by the common one-frame protocol. The table summarizes the evaluated implementations. Model Retained interaction mechanism Protocol adaptation ReViT canonical local frames and rotationally equivariant transformer blocks Physical-feature adapter constructs the vector/scalar layout from each task; one-frame residual decoder. Flowers learned multihead coordinate warps and latent resampling The same state, context, boundary sampler, and autoregressive evaluator are used; read coordinates are learned from forecast data. NIPS-AR channel-independent Fourier kernel and linear attention Released spatial block placed in a recurrent one-frame wrapper; reported as an autoregressive adaptation. RieszNO / RNO-3D Riesz directional representation plus spectral mixing Released 2D construction; direct axial 3D extension is explicitly marked as a dimensional adaptation. For every baseline reported in the matched comparison, the implementation uses the model-specific adaptation listed in Table 10. Its checkpoint and manifest are verified against the common array files, sampled windows, training-only statistics, K=4K=4 loss, optimizer schedule, effective batch, update count, validation selection, and H=40H=40 evaluator. Micro-batch size is adjusted only for memory. Each experiment manifest records the source repository revision, architecture configuration, parameter count, peak allocated memory, and forward time for that run. Appendix H Metrics and result records For each test window, errors are evaluated in physical units. At rollout step h, the relative L2L^2 error is EL2(h)=‖q^h−qh‖2max‖qh‖2,ε,ε=10−8.E_L^2(h)= \| q_h-q_h\|_2 \\|q_h\|_2, \, =10^-8. (47) The rollout mean is H−1∑h=1HEL2(h)H^-1 _h=1^HE_L^2(h). The reported spectral error is Espec(h)=[‖|ℱq^h|−|ℱqh|‖22max‖ℱqh‖22,ε]1/2.E_spec(h)= [ \||F_P q_h|-|F_Pq_h| \|_2^2 \\|F_Pq_h\|_2^2, \ ]^1/2. (48) Thus the per-window spectral metric is the square root of the spectral ratio used in training. Relative H1H^1 adds the discrete derivative energies to the squared field norms in numerator and denominator before taking the square root. Spatial correlation is computed after subtracting each field’s spatial mean. Persistence uses q^h=q0 q_h=q_0, and its skill score is Skill(h)=1−EL2model(h)maxEL2persistence(h),ε.Skill(h)=1- E_L^2^model(h) \E_L^2^persistence(h), \. (49) Table 11: Task-specific physical metrics recorded by the common evaluator. These metrics are reported in addition to relative L2L^2 error. Task Additional records Kolmogorov / SQG correlation, spectral error, variance error, and kinetic-energy or scalar-variance error Vlasov–Poisson total-mass drift, negative-density fraction, charge-density error, and electric-energy error Keller–Segel mass drift, negative-mass fraction, peak and variance error, and screened free-energy error Porous convection anomaly error, plume-mask IoU, top dissolution-flux error, plume-depth error, scalar-dissipation error, variance error, and fractions below zero or above one Errors are computed for each trajectory and rollout step. Within each seed, results are averaged over the 12 prescribed test windows. Tables 1 and 2 report the mean and sample standard deviation across the three seed means, whereas Table 3 reports seed 0. Appendix I Expanded route-mechanism study I.1 What each intervention changes Table 12: Fixed-weight route interventions. Network weights, construction of carrier features, centered branches, and prescribed test windows are held fixed. Mode Operation on displacement d Question tested Correct d from the equation adapter reference computation Zero d←0d← 0 Is the carrier used as a coordinate rather than only as an input feature? Reverse d←−d←-d Does route orientation matter? Shift deterministic spatial roll of d Does alignment between route and current state matter? At each rollout step, the carrier and departure displacement are recomputed from that run’s current predicted state. The chosen intervention is applied to this new displacement before pooling to each feature scale and sampling the features. For the Shift intervention, the fixed offsets are (⌊Ny/3⌋,⌊Nx/5⌋)( N_y/3 , N_x/5 ) in 2D and (⌊Nz/5⌋,⌊Ny/4⌋,⌊Nx/3⌋)( N_z/5 , N_y/4 , N_x/3 ) in 3D. The displacement values continue to evolve with the predicted state at every rollout step. Final intervention errors are reported in Table 2; Figure 2c–d shows the corresponding error ratios and rollout profiles. I.2 Independently trained controls All variants share the core architecture within each dimension. Disabling carrier features sets their input channels to zero; disabling routing sets the latent displacement to zero. The centered and routed branches remain in place, giving the four factorial variants identical parameter counts. The Shifted equation-specified route control applies the same deterministic shift as above throughout both training and testing. Learned route. The learned-route displacement head has a pointwise input projection, one local residual block, and a pointwise output projection. Its width is 24 in 2D and 12 in 3D; it outputs two or three displacement components. The final projection is initialized to zero. For each component, dϕ,j(q)=ajtanh([Hϕ(xμ(q))]j),aj=max3RMStrain(dj),Δxj.d_φ,j(q)=a_j \! ([H_φ(x_μ(q))]_j ), a_j= \3\,RMS_train(d_j), x_j\. (50) Here djd_j is the equation-specified departure displacement on training states and Δxj x_j is the full-grid spacing. The 2D scale estimate uses 64 training states sampled with replacement using seed 23071; the 3D estimate uses 16 evenly spaced entries in the ordered list of training trajectory/frame pairs. The head receives the same normalized input as the forecasting core and is optimized jointly with it through the forecasting loss. State advection + residual. The State advection + residual control first computes qadv=ℐ[q]([q])q_adv=I[q]( [q]) and then applies q^+=q+Πμ(qadv−q+σΔq,μRθ(xμ(q))). q^+=q+ _μ\! (q_adv-q+ _ q,μR_θ(x_μ(q)) ). (51) The centered network receives the current state, carrier, and context, and all latent displacements are set to zero. It is trained from scratch with the same loss, training budget, and checkpoint rule. The task projection acts on the complete state increment, including qadv−q_adv-q. The diagnostic carrier and departure map are recomputed from the current prediction at every step. I.3 Visualization sample and transect For Figure 2e, the selected test window maximizes the reduction in final relative L2L^2 error from Carrier feature only to GRNO among the twelve prescribed windows. Within this window, the vertical line maximizes the reduction in absolute field error averaged over h=1,…,40h=1,…,40 and the vertical grid coordinate. Both independently trained models use the same window and line. The three concentration maps use a common color scale, and the same reference concentration contour is overlaid as a dashed line on each map. Appendix J Additional Visualizations J.1 Additional qualitative forecasting results Figures 3–7 compare reference fields and pointwise absolute errors at ten horizons spanning the 40-step free rollout, using GRNO and the strongest baseline in Table 1. Figure 3: Long-horizon qualitative comparison on Kolmogorov flow. The first row shows the reference vorticity field at ten horizons spanning the 40-step free rollout. The second and third rows show the pointwise absolute errors of GRNO and Flowers, respectively. The error scale is shared by both models and all displayed horizons. Figure 4: Long-horizon qualitative comparison on dissipative SQG. The first row shows the reference scalar field at ten horizons spanning the 40-step free rollout. The second and third rows show the pointwise absolute errors of GRNO and ReViT, respectively. The error scale is shared by both models and all displayed horizons. Figure 5: Long-horizon qualitative comparison on Vlasov–Poisson. The first row shows the reference phase-space distribution at ten horizons spanning the 40-step free rollout. The second and third rows show the pointwise absolute errors of GRNO and Flowers, respectively. The error scale is shared by both models and all displayed horizons. Figure 6: Long-horizon qualitative comparison on Keller–Segel. The first row shows the reference density field on the fixed visualization slice at ten horizons spanning the 40-step free rollout. The second and third rows show the pointwise absolute errors of GRNO and Flowers, respectively. The error scale is shared by both models and all displayed horizons. Figure 7: Long-horizon qualitative comparison on porous convection. The first row shows the reference concentration field on the fixed visualization slice at ten horizons spanning the 40-step free rollout. The second and third rows show the pointwise absolute errors of GRNO and Flowers, respectively. The error scale is shared by both models and all displayed horizons. J.2 Qualitative effects of route interventions Figures 8–12 show the fixed-weight route interventions at h=40h=40. For each system, the fitted GRNO checkpoint is held fixed, and only the departure displacement is changed. Figure 8: Route interventions on Kolmogorov flow at h=40h=40. The top row shows the reference field and predictions using the correct, zero, reversed, and shifted routes. The bottom row shows the corresponding pointwise absolute errors. All predictions use the same fitted GRNO checkpoint; only the departure route is changed. Figure 9: Route interventions on dissipative SQG at h=40h=40. The top row shows the reference field and predictions using the correct, zero, reversed, and shifted routes. The bottom row shows the corresponding pointwise absolute errors. All predictions use the same fitted GRNO checkpoint; only the departure route is changed. Figure 10: Route interventions on Vlasov–Poisson at h=40h=40. The top row shows the reference phase-space distribution and predictions under the Correct, Zero, Reverse, and Shift interventions. The bottom row shows the corresponding pointwise absolute errors. All predictions use the same fitted GRNO checkpoint; only the departure route is changed. Figure 11: Route interventions on Keller–Segel at h=40h=40. The top row shows the reference field and predictions on the fixed visualization slice using the correct, zero, reversed, and shifted routes. The bottom row shows the corresponding pointwise absolute errors. All predictions use the same fitted GRNO checkpoint; only the departure route is changed. Figure 12: Route interventions on porous convection at h=40h=40. The top row shows the reference concentration field and predictions on the fixed visualization slice using the correct, zero, reversed, and shifted routes. The bottom row shows the corresponding pointwise absolute errors. All predictions use the same fitted GRNO checkpoint; only the departure route is changed.