Paper deep dive
Read, Write, Relax: Why Neural PDE Surrogates Need Both Global and Local Processing
Anuj Kumar, Heiko Zimmermann, Josiah Bjorgaard, Jacan Chaplais, Nikolaos Bouklas, Matteo Salvador, Alexander Lavin
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 92%
Last extracted: 8/25/2026, 7:05:15 AM
Summary
The paper introduces Read-Write-Relax (RWR), a unified neural surrogate architecture for mesh-based Partial Differential Equation (PDE) simulations that combines global latent-attention mechanisms with local message-passing operations. The authors demonstrate that global models (like latent-attention) act as spatial low-pass filters, correcting low-frequency errors, while local models (like MeshGraphNets) act as high-pass filters, correcting high-frequency errors. By interleaving these two processes, RWR lowers error across the entire spectrum, achieving superior accuracy, data efficiency, and scalability on industrial-grade benchmarks compared to existing single-family approaches.
Entities (9)
Relation Signals (10)
Read-Write-Relax → combines → Message Passing
confidence 95% · We introduce Read-Write-Relax (RWR), which interleaves latent attention with message-passing relaxation under a unified formulation.
Read-Write-Relax → combines → Latent Attention
confidence 95% · We introduce Read-Write-Relax (RWR), which interleaves latent attention with message-passing relaxation under a unified formulation.
Read-Write-Relax → evaluatedon → Branched Pipe
confidence 95% · All analyses of Sec. 5 run on this benchmark [Branched pipe].
Read-Write-Relax → evaluatedon → Centrifugal Pump
confidence 95% · Centrifugal pump... both appear in the benchmarks of Sec. 5.6.
Latent Attention → actsas → Spatial Low-Pass Filter
confidence 92% · In global approaches, latent-token attention acts as a spatial low-pass filter...
MeshGraphNets → ismemberof → Local Models
confidence 92% · Local models such as MeshGraphNets (MGN) [34] update each node from its mesh neighbors...
Message Passing → corrects → High-Frequency Errors
confidence 90% · local message passing lacks the global reach... one corrects errors at the lower end of the spectrum, the other at the higher end...
Read-Write-Relax → outperforms →
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Recent mesh-based simulation advances have, in no small part, relied on neural surrogates of two distinct families: global models that route information through a small set of latent tokens, and local models that perform message passing across mesh edges. Consistent with both classes is the inability to perform beyond low-dimensional problems and small-scale or oversimplified meshes, the simulation regimes where industrial problems reside. Our work shows this explicitly and presents a unified formulation. In global approaches, latent-token attention acts as a spatial low-pass filter, while local message passing lacks the global reach necessary to propagate information across large mesh spaces. Viewed through the error, the two operators are the halves of a multigrid cycle: one corrects errors at the lower end of the spectrum, the other at the higher end, and neither can do the other's job. We introduce Read-Write-Relax (RWR), which interleaves latent attention with message-passing relaxation under a unified formulation. The interleaved processor lowers error across the entire spectrum, making RWR the most accurate model in nearly every comparison across our industrial and public benchmarks. It is also markedly data-efficient in the scarce-data regimes, accurate on the engineering quantities of interest, and scales full-field predictions to challenging, large-scale problems.
Tags
Links
- Source: https://arxiv.org/abs/2608.21677v1
- Canonical: https://arxiv.org/abs/2608.21677v1
Trouble viewing inline? Open PDF directly →
Full Text
71,709 characters extracted from source content.
Expand or collapse full text
Read, Write, Relax: Why Neural PDE Surrogates Need Both Global and Local Processing Read, Write, Relax: Why Neural PDE Surrogates Need Both Global and Local Processing Anuj Kumar Heiko Zimmermann Josiah Bjorgaard Jacan Chaplais Nikolaos Bouklas Matteo Salvador Alexander Lavin August 2026 Abstract Recent mesh-based simulation advances have, in no small part, relied on neural surrogates of two distinct families: global models that route information through a small set of latent tokens, and local models that perform message passing across mesh edges. Consistent with both classes is the inability to perform beyond low-dimensional problems and small-scale or oversimplified meshes, the simulation regimes where industrial problems reside. Our work shows this explicitly and presents a unified formulation. In global approaches, latent-token attention acts as a spatial low-pass filter, while local message passing lacks the global reach necessary to propagate information across large mesh spaces. Viewed through the error, the two operators are the halves of a multigrid cycle: one corrects errors at the lower end of the spectrum, the other at the higher end, and neither can do the other’s job. We introduce Read-Write-Relax (RWR), which interleaves latent attention with message-passing relaxation under a unified formulation. The interleaved processor lowers error across the entire spectrum, making RWR the most accurate model in nearly every comparison across our industrial and public benchmarks. It is also markedly data-efficient in the scarce-data regimes, accurate on the engineering quantities of interest, and scales full-field predictions to challenging, large-scale problems. Introduction Learned surrogates for mesh-based simulation have converged on two designs. Global models compress the discretized field into a small set of latent tokens, process the tokens, and expand the result back onto the mesh, a pattern shared by Perceiver-style encoders [22, 21], Transolver’s physics attention [54], and FLARE’s low-rank routing [36]; we refer to these global models as latent-attention models. Local models such as MeshGraphNets (MGN) [34] update each node from its mesh neighbors, realizing learnable stencils that resolve structure down to the mesh scale but communicate only within a receptive field that grows by one hop per layer [3]. Steady-state problems stress both designs at once. At convergence, the solution everywhere is shaped by the boundary and operating conditions together with the geometric parameterization, so a surrogate must propagate information across the entire domain in a single forward pass rather than over many rollout steps, while simultaneously resolving the steep local gradients of the solution. Global communication and local resolution are therefore both required, and each family delivers one while falling short on the other. The local half carries particular weight in practice, because the outputs that engineering decisions rest on are derivative-based, among them wall shear stress, vorticity, and integrated losses computed from them; we call them engineering quantities. The way surrogates for Computer-Aided Engineering (CAE) are trained and evaluated magnifies both shortfalls, the local family’s limited communication and the global family’s limited resolution. On the training side, every sample is the output of an expensive solver run, with industrial-grade Reynolds-Averaged Navier–Stokes (RANS) and hybrid RANS–large-eddy simulations costing hundreds to thousands of core-hours per geometry, so CAE datasets number in the hundreds of samples rather than the millions common elsewhere in deep learning [7, 4, 5]. Data efficiency is therefore a requirement rather than a convenience. On the evaluation side, progress is predominantly reported on bulk metrics such as the Mean Absolute Error (MAE), relative L2L_2 error, and R2R^2, and these are documented to be poor proxies for the engineering quantities just defined [44, 7]. We take up these requirements on two fronts, an analysis of why each family fails and an architecture built on the answer. The analysis compares a message-passing model, MGN, with a latent-attention model we call the Read-Write Perceiver (RWP), tracking their error spectra over the course of training. The two designs fail from opposite ends of the spectrum, and read through multigrid [9, 10], they occupy the two complementary roles of a cycle. Read-Write-Relax (RWR), our architecture, interleaves the two mechanisms and takes its name from this reading. The read and the write act as the restriction and prolongation of a coarse correction in latent space, and the message-passing sweeps act as the learned relaxation. The combined processor recovers the entire spectrum, carries the accuracy into the gradients and wall quantities that make up the engineering quantities, and holds where it matters most in practice, at small training budgets, where the latent-attention baselines degrade fastest. Contributions 1. The analysis. We place recent global models in one family, Bottlenecked Self-Attention (BSA), of which our RWP is itself a member, and characterize MGN, RWP, and RWR under a single spectral lens on an industrial steady-flow benchmark. This yields a unified analysis of their complementary failure modes across the error spectrum (Sec. 5). 2. The architecture. RWR interleaves latent attention with message-passing relaxation under a unified configuration (m1,r,m2)(m_1,r,m_2) that prescribes how much local and global processing the model performs, with MGN and RWP recovered as special cases. Geometry-anchored queries and boundary-only encoding, which lets the attention read only the information-dense boundary points, keep the design scalable to large industrial problems, and the resulting model is markedly data-efficient, delivering accurate engineering quantities from far fewer training samples than the latent-attention baselines. In addition, across the benchmarks, RWR is the most accurate model in nearly every comparison (Secs. 3, 5). Related Work Attention and message passing on graphs The conflict between local message passing and global information exchange is well documented in graph learning. Message passing cannot reach beyond its receptive field, and even within it, information from exponentially growing neighborhoods is compressed into fixed-size vectors [3], a failure mode known as over-squashing, later linked to negatively curved edges in the graph [46]. On the other hand, graph transformers discard message passing altogether and attend over all node pairs, relying on spectral or structural encodings to reinstate the topology that attention alone would ignore [26, 58]. At the opposite and cheapest end, a single virtual node grafted onto a message-passing network acts as a global scratch space [16], a mechanism that provably approximates linear attention [11]. Hybrid approaches define the middle ground, where some run a Graph Neural Network (GNN) first and a transformer afterward, so that local structure is summarized before global reasoning [55], while others run message passing and global attention in parallel within every layer [40, 43]. Our architecture follows this lineage but alternates the two mechanisms across the depth of the network, with the global pathway realized as cross-attention to a small latent set rather than attention over all nodes. Local and global processing in PDE surrogates The same problem arises in mesh-based PDE surrogates. MGN-style processors [34] need message-passing depth proportional to the mesh diameter, which is particularly limiting for steady and elliptic problems where boundary information influences the interior globally [50] and accurate prediction requires long-range exchange that deep stacks capture poorly and train slowly [41, 17]. The same issue appears in time-dependent regimes: parabolic problems require stable long-horizon propagation, and hyperbolic problems require sharp, non-dissipative transport of features like shock waves, both hard for a processor whose receptive field scales only with depth. Early responses built multiscale hierarchies [15, 12], and EAGLE introduced global attention over a pooled, clustered mesh [23]. A second family adopts a stacked composition that aggregates locally once and then processes globally with transformers, as in GINO [30] and GAOT [53]. A third line of work keeps both mechanisms inside the processor. GITO and FRGT combine the two mechanisms in parallel and stacked form, respectively [39, 14]. Most relevant is MeshTransolver, which wraps a Transolver-style token-attention core between message-passing pre- and post-processing stages [13], a sandwich composition that concurrent work adopts as well [19]. These designs fix the composition in advance as parallel, stacked, or a single wrapped core, and are validated by end-task accuracy, largely leaving open why the combination outperforms either mechanism alone. RWR alternates relaxation and latent-attention iterations throughout the network and exposes the mix as a configuration, and the analysis of Sec. 5 supplies the missing account in spectral terms. Latent-token surrogates A parallel line of latent-token surrogates, including Perceiver IO [21], Transolver and its successors [54, 32, 1], LNO [52], UPT [2], AROMA [42], and FLARE [36], compresses the field into a small set of tokens and processes it there. We return to these models in the analysis (Sec. 5), where we show that they instantiate a single BSA template together with our own RWP. The Model: Read-Write-Relax Figure 1: One RWR block. The encoder and decoder are pointwise MLPs; the processor runs m1m_1 message-passing sweeps (relax), r read-write iterations against L persistent latent tokens, and m2m_2 further sweeps; nBn_B such blocks are stacked. Setting r=0r=0 recovers MGN, and setting m1=m2=0m_1=m_2=0 recovers RWP. MP denotes message passing, CA cross-attention, and SA self-attention. Our surrogate follows the encode-process-decode organization of MGN [34]: node and edge features are lifted by pointwise encoders, transformed by a processor, and projected to output fields by a decoder. RWR is the processor. It consists of nBn_B blocks, each parameterized by a configuration (m1,r,m2)(m_1,r,m_2) that runs m1m_1 message-passing sweeps, r read-write iterations against a set of L latent tokens, and m2m_2 further sweeps (Fig. 1). The configuration prescribes how much local and global processing the model performs. In particular, r=0r=0 leaves only message passing and recovers MGN, while setting m1=m2=0m_1=m_2=0 removes the relaxation sweeps entirely and yields the RWP, the latent-attention model we analyze alongside RWR throughout the paper. All three models in our experiments are therefore the same implementation at different configurations. Let U∈ℝN×CU ^N× C denote the node features, P their positional encodings, and Z∈ℝL×CZ ^L× C the latents. One read-write iteration applies read: : Z←Z+CA(Q=Z,KV=U|ℛ+P|ℛ) Z← Z+CA\! (Q=Z,\;KV=U|_R+P|_R ) (1) write: : U←U+CA(Q=U+P,KV=Z), U← U+CA\! (Q=U+P,\;KV=Z ), where CACA is pre-norm cross-attention with a feed-forward layer and ℛR is the read set. An optional latent self-attention tower can be inserted between the read and the write. The latents persist across the r iterations of a block, so each read refines the state accumulated by earlier ones rather than rebuilding it, a design shared with the read-process-write generator of 20. The relaxation sweeps are standard message-passing updates with residual connections, in which edge features are updated from their endpoint nodes, and node features are updated from their aggregated incident edges [34]. The names record the multigrid reading introduced above. The read acts as a restriction onto an L-dimensional coarse space, the write as a prolongation back to the mesh, and together with the latent update they form a coarse correction; the message-passing sweeps act as learned relaxation. Geometry-anchored queries The latent queries are anchored in the geometry rather than learned as free vectors. Each latent possesses a learnable anchor coordinate in physical space, and its query is produced by passing this coordinate through the same random-Fourier-feature embedding [45] that encodes node positions for the read and the write, so queries and keys meet in a single coordinate space and each latent costs three parameters instead of a full feature vector. One embedding is shared across all blocks, which ties the frequency content of the positional encoding across the processor. Boundary-only encoding Because the read and the write are separate attention maps, the read set ℛR need not cover the full mesh. For steady problems we can optionally restrict it to the boundary points, where the geometry, the boundary conditions, and hence the identity of the problem are concentrated. The choice has a classical justification, since for elliptic problems the interior solution is determined by boundary data through the Green representation [25], an observation that boundary-embedded neural operators exploit directly [50], and it reduces the read cost from O(NL)O(NL) to O(NbL)O(N_bL) with Nb≪N_b N boundary points. The write continues to address every node, since the output lives on the full mesh. Datasets We evaluate the surrogate models on two industrial steady-flow problems (Fig. 2), the second of which we use at two extents. Figure 2: The two simulation environments, branched pipe (left) and centrifugal pump (right), showing a representative velocity field. Branched pipe Steady flow through a Y-shaped Heating-Ventilation-Air-Conditioning (HVAC) duct junction whose per-channel flow rates are controlled by baffle plates; the sharp turning angles at the baffles make the flow at the junction highly turbulent. Samples vary the duct geometry, the baffle angles, and the inlet velocity, and the modeled fields are the pressure and the two in-plane velocity components (p, uxu_x, uzu_z) on meshes of roughly 60k cells. All analyses of Sec. 5 run on this benchmark. Centrifugal pump Steady flow through a centrifugal pump, in which a rotating impeller draws fluid in through the eye of a stationary casing and discharges it radially through the volute; the rotating blades raise both the pressure and the velocity of the fluid, and the rotation is modeled with a moving reference frame. Samples vary the inlet and impeller radii, the blade pitch angle, and the inlet velocity, and the modeled fields are the three velocity components and the pressure. We use it at two extents, the casing region alone, the stationary zone of the simulation with roughly 270k cells, and the full production mesh including the rotating impeller zone with roughly 1.4M cells; both appear in the benchmarks of Sec. 5.6. Both datasets are generated with a finite-volume RANS setup. Simulation configurations, design-parameter ranges, and dataset splits are given in App. A. Results and Discussion All analyses in this section, except for the ones in Sec. 5.6, use the branched pipe benchmark. MGN, RWP, and RWR are trained as configurations of the same RWR processor at matched parameter counts and identical training setups; where the comparison concerns the latent-attention family more broadly, we also train Transolver and GeoTransolver in the same pipeline under the same budget. Errors are reported in physical space and spectrally. We use the error spectrum Eerr(k)E_err(k), computed on a mid-plane slice of the domain, and gradient and wall metrics for the engineering quantities; band-resolved errors in the convention of PDEBench [44] are reported in the appendix (Table 8). Precise metric definitions (App. B), together with model configurations and training details (App. D.4), are given in the appendix as well. 5.1 One Global Family: BSA The ideal global mixer is full self-attention over all N mesh points, which is message passing on the complete graph and couples every point to every other in a single layer. Its O(N2)O(N^2) cost is prohibitive at industrial mesh sizes, so the recent global models buy the coupling at a discount by routing it through L≪NL N latent tokens, and the discount has one common form. In the notation of Sec. 3, one layer of the family computes Z=AV(U),Z′=g(Z),U←U+BZ′Z=A\,V(U), Z =g(Z), U← U+B\,Z (2) with a row-stochastic read A∈ℝL×NA ^L× N, a value projection V, a latent processor g, and a row-stochastic write B∈ℝN×LB ^N× L. We call this template Bottlenecked Self-Attention (BSA). The induced point-to-point mixing matrix W=BAW=BA has rank at most L, a bound FLARE states for its own mixing operator [36] and which holds across the family: each layer can represent the mesh with at most L spatial basis fields. The template thus gives the family a common vocabulary, a read that compresses the field into tokens, a latent update, and a write that broadcasts the update back, and the recent global models differ only in how each of these three roles is parameterized. Table 1 places them in the template. Perceiver IO reads once into latents and decodes once at the end [21]. Transolver’s slices are latents, since its slice weights are a softmax read, its slice tokens the weighted accumulation of the points, self-attention runs over the L slice tokens, and the same tied weights broadcast the result back [54]. LNO’s physics-cross-attention projects the field onto L learnable position embeddings and inverts the projection at the output [52]. UPT aggregates the mesh onto supernodes and pools them into latent tokens that a transformer processes end to end [2]. AROMA encodes into L tokens refined by a latent diffusion transformer, decoding through query-local neural fields [42]. FLARE routes through per-head latent tokens with no latent processing at all [36]. Our RWP belongs to the same family, built for repeated use inside a processor with a decoupled read and write and latents that persist across iterations. The kind of update this rank-bounded template permits is illustrated empirically by the training dynamics below. read process write write-back freq. Perceiver IO queries SA out queries once Transolver family slice weights SA tied weights layer LNO learned PE SA inverse proj. once UPT MP + pool SA perceiver once AROMA cross-attn DiT query-local once FLARE queries/head — roles swapped layer RWP (ours) geo-anchored SA (opt.) decoupled iter. Table 1: The BSA family, Eq. (2). All members route the field update through L≪NL N latent tokens and inherit the rank bound; they differ in how the read and write are parameterized and how often the write recurs (once at decode, every layer, or every iteration). RWP’s latents persist across iterations within a block, and its read may further be restricted to the boundary points. The Transolver family row covers Transolver, Transolver++, and GeoTransolver, which share the slice read and write. SA denotes self-attention over the latent tokens, DiT a diffusion transformer, PE a positional embedding, and MP message passing; in FLARE the read and write swap the query and key-value roles between the same projections. Figure 3: Normalized error spectra over training on held-out branched pipe cases. Each panel shows Eerr(k)/Eerr(0)(k) E_err(k)/E_err^(0)(k), the error amplitude at wavenumber k normalized by the initial error of RWP, for the two velocity components (rows) under MGN, RWP, and RWR (columns) at five training epochs. Dashed lines mark the low/mid/high band cutoffs. 5.2 Latent Attention Behaves as a Spatial Low-Pass Filter Fig. 3 tracks the error spectrum of each model over the course of training. The RWP column shows a strong and persistent tilt. Error in the low band contracts by more than an order of magnitude and continues to contract as training proceeds, while the mid- and high-band errors stall early and settle several times higher. Whatever the model learns, it learns overwhelmingly at large spatial scales. The latent-attention processor behaves as a spatial low-pass filter. The band-resolved training curves (Fig. 6 in the appendix) show the same picture, with RWP’s high band separating from its low band early in training and never rejoining it. We attribute this behavior to several compounding factors rather than to a single cause. The first is the bottleneck itself, since a BSA update spans at most L spatial basis fields per layer, and covering the domain with L read profiles forces every latent to average over a finite region, attenuating the scales finer than that region. The second is attention, independent of any bottleneck, since a row-stochastic attention matrix is a kernel smoother [47] and provably a low-pass filter under repeated application [51, 33]. The third is the loss, since the spectra of PDE solutions decay with wavenumber, so the low wavenumbers dominate a Euclidean objective while the high wavenumbers contribute almost nothing, a documented driver of over-smoothed neural-operator predictions [31, 37, 24]. The fourth, not independent of the third, is spectral bias, since networks fit low frequencies before high ones [38, 56]. The last three factors apply to any attention-based surrogate trained with a mean-squared objective, so our measurements establish consistency with this account rather than isolating the bottleneck’s contribution. If the bottleneck were the sole cause, widening it should recover the missing scales. Sweeping the number of latents from 8 to 512 with all other settings fixed leaves every band almost flat (Fig. 5 in the appendix), consistent with the saturation reported for Perceiver IO, the Set Transformer, and low-rank spatial attention [27, 21, 57]. Training does not exploit the capacity a wider bottleneck offers, as the loss- and optimization-side factors predict; the practical route to the high band is a second operator class rather than a wider bottleneck. 5.3 Message Passing Fails from the Other End The MGN column of Fig. 3 shows a different failure. The error contracts quickly over the first 150 epochs, with the larger early gains toward the middle and higher wavenumbers, and then plateaus; from epoch 300 onward the curves are nearly indistinguishable, and no band, including the lowest, improves further. Where RWP’s low band keeps contracting throughout training, MGN stalls everywhere, at an error several times higher across the spectrum. This is the failure the message-passing literature predicts for a steady problem. A message-passing layer is a local stencil that resolves structure near the mesh scale cheaply, which is consistent with the early gains sitting in the higher bands, but information travels only one edge per layer, and on steady and elliptic problems, where boundary data shapes the interior globally, the number of sweeps must grow with the ratio of domain size to mesh spacing before solutions stop degrading [8]. On our meshes the graph diameter far exceeds the depth of a typical MGN processor, so the domain-scale content is out of reach; nor is deepening the stack a clean fix, since squeezing an exponentially growing neighborhood through fixed-width messages degrades the information that does arrive. Since the nonlinear couplings of the flow tie the wavenumber bands together, the stalled low band plausibly holds back the rest of the spectrum, though our runs do not isolate this mechanism. The two pure processors therefore fail from opposite ends, one shedding the scales it cannot afford, the other starving the scales it cannot reach, and the remedy each needs is what the other provides. 5.4 Combining Local and Global Processing The RWR column of Fig. 3 completes the picture. The interleaved processor contracts the error across the entire spectrum throughout training, with neither the plateau of MGN nor the high-band stagnation of RWP. Table 2 compares the fully trained RWR with the three latent-attention baselines at matched training budget, and the errors split along the expected lines. The upper rows report each velocity component and its spatial gradients, the lower rows the derived vorticity and wall shear stress, and the bracketed factors give each baseline’s error relative to RWR. Transolver GeoTransolver RWP RWR uxu_x 0.078 [1.01] 0.098 [1.28] 0.090 [1.17] 0.077 uxu_x grad. 0.174 [1.18] 0.226 [1.54] 0.168 [1.14] 0.147 uzu_z 0.080 [1.11] 0.104 [1.45] 0.091 [1.27] 0.072 uzu_z grad. 0.207 [1.27] 0.286 [1.75] 0.218 [1.33] 0.164 ωy _y 0.143 [1.23] 0.212 [1.82] 0.138 [1.19] 0.117 τw,x _w,x 0.170 [1.30] 0.246 [1.88] 0.157 [1.20] 0.131 τw,z _w,z 0.231 [1.42] 0.363 [2.23] 0.247 [1.51] 0.163 Table 2: Held-out test errors on the branched pipe at matched training budget: relative L2L_2 error of each velocity component and of its spatial gradients, and relative L2L_2 error of the derived quantities, the in-plane vorticity ωy _y and the wall shear stress components τw,i=μ∂ui/∂n _w,i=μ\,∂ u_i/∂ n evaluated at wall boundary points. Bracketed values give the ratio of each model’s error to RWR’s. Metric definitions are given in App. B. Figure 4: Held-out test error on the branched pipe as a function of the training-set size for the bulk relative L2L_2 (top row), the gradient relative L2L_2 (middle row), and the wall shear stress (bottom row) of the two velocity components. Moving from the bulk errors to the gradients and on to the derived wall quantities, the latent-attention models fall progressively further behind RWR, and Fig. 3 shows the same progression spectrally, with the deficit growing toward the high wavenumbers. In the opposite direction the ordering softens, since the strongest latent-attention baselines still match RWR at the lowest wavenumbers (Table 8 in the appendix). Both halves are consistent with the preceding subsections. The latent-attention path remains the right tool for the smooth content, relaxation is what buys the fine scales, and the interleaved processor is the only one that holds both ends at once. The gradient and wall quantities deserve the emphasis, because they are the engineering quantities defined in the introduction, and bulk errors are documented to be poor proxies for them. It is precisely here that the combination delivers its largest margins, and the pointwise diagnostics (Figs. 7 and 8 in the appendix) read the same, with RWR attaining the highest wall-normal derivative parity of the four models and the lowest near-wall error on the velocity components. The results so far establish that both operator classes are needed, but not that they must alternate, since a stacked composition that runs all relaxation before or after the latent-attention iterations contains the same two mechanisms. Comparing four allocations of an identical budget of eight relaxation sweeps and four latent-attention iterations, the two interleaved configurations improve on the two stacked orderings on every field, though all four land within a narrow range (Table 9 in the appendix). Once both operator classes are present at a sensible budget, the composition refines the picture rather than redrawing it. 5.5 RWR Is Data-Efficient The comparisons so far hold the training set fixed at one hundred samples. Fig. 4 varies the budget from 50 to 300 and tracks the bulk, gradient, and wall shear errors of the four models. RWR attains the lowest error on every quantity at every budget, and its curves are the flattest, especially for the engineering quantities, so the margins are widest where data is scarcest. Trained on 50 samples, RWR already matches or improves on every latent-attention baseline trained on twice as much data on the wall shear stress components. Since a training sample here is a full solver run, the scarce side of the sweep is where these surrogates actually operate, and it is there that RWR pays the most. 5.6 Benchmarks We train Transolver, GeoTransolver, RWP, and RWR on the branched pipe and the two centrifugal pump datasets, each at three training-set sizes. Table 4 reports the mean test Normalized MAE (NMAE) over the output fields; per-field values, parameter counts, memory, and training cost are given in Table 12 in the appendix. The ordering of Fig. 4 on the branched pipe carries over to the industry-scale meshes, RWR attains the lowest mean error in every cell, and the margins widen toward the smallest budgets, extending the data-efficiency picture of Sec. 5.5 across mesh scales. Qualitative predictions with pointwise error maps on held-out test cases of all three datasets are shown in Figs. 9–14 in the appendix. Branched pipe Pump casing Pump full Model (60k) (270k) (1.4M) Transolver 11.3 || 7.5 || 4.8 15.6 || 10.9 || 7.9 37.1 || 23.2 || 15.9 GeoTransolver 10.5 || 8.3 || 5.5 12.1 || 09.3 || 6.9 32.2 || 21.1 || 14.1 RWP 10.3 || 7.3 || 5.0 11.4 || 08.1 || 6.3 21.8 || 17.3 || 14.4 RWR 08.4 || 6.9 || 4.7 10.8 || 07.7 || 5.4 19.3 || 15.2 || 12.7 Table 3: Benchmarks on the branched pipe and pump datasets. Each cell reports the mean test NMAE (×10−2× 10^-2) over the output fields at 50 || 100 || 300 training samples, for models trained in our pipeline; model and training details are given in App. E. Bold marks the best value per dataset and budget. Geo-FNO Model AirfRANS AhmedML airfoil pipe Darcy MGN 141 122 60.4 205 98.7 MeshTransolver 19.6 9.8 2.32 0.403 0.678 Transolver 21.4 9.3 1.99 0.643 0.356 GeoTransolver 20.9 8.6 2.15 0.845 0.412 RWP 29.9 10.8 1.78 0.554 0.218 RWR 15.7 7.2 1.58 0.409 0.190 Table 4: Best validation MSE for additional benchmarks (×10−3× 10^-3) on normalized fields. Bold marks the best value per dataset. We trained models additionally on AirfRANS [7], the surface mesh component of the AhmedML [4] dataset, and the airfoil, pipe, and Darcy problems of the Geo-FNO suite [28, 29]. These datasets provide a combination of irregular/regular 2D meshes and a 3D surface mesh to complement the 3D meshes of Table 4. Table 4 reports the best validation Mean Squared Error (MSE) on normalized fields for MGN, MeshTransolver, Transolver, GeoTransolver, RWP, and RWR; datasets, model configurations, and the full protocols are given in App. F. The two tables read consistently. Since these benchmarks use a single random seed, we make only general claims about the comparisons. RWR is the best-performing model on all datasets of Table 4 and on all additional benchmarks except the Geo-FNO pipe, where MeshTransolver is marginally better. As discussed in Sec. 5.4 and demonstrated in Fig. 3, we expect both RWP and MGN to perform worse than RWR. Indeed, MGN is consistently and significantly the worst-performing model on all benchmarks. RWP performs below RWR in all cases, yet remains in the performance range of Transolver and GeoTransolver, exceeding both on the three Geo-FNO problems. Finally, we see that the MGN-Transolver hybrid structure of MeshTransolver offers competitive performance on these additional benchmarks. Importantly, these comparisons are on bulk metrics alone. Table 2 showed that where RWR improves the bulk error over a pure latent-attention model, its gains on the gradient and wall metrics are larger still, so we expect the margins reported here to widen further on those engineering quantities (Sec. 5.4). Future work will demonstrate detailed comparisons of additional engineering quantities on industrial-scale datasets for further applications using the RWR architecture. Conclusion Global latent-token surrogates and local message passing fail on complementary ends of the error spectrum. We placed the recent global models in a single bottlenecked self-attention template and showed that on an industrial steady-flow benchmark they learn overwhelmingly at large spatial scales, a behavior consistent with several compounding causes, and one that widening the bottleneck does not repair. Message passing fails in the opposite way, gaining early on the fine scales but stalling across the spectrum when its receptive field cannot cover the domain. RWR interleaves the two mechanisms as multigrid alternates coarse correction and relaxation, exposes the mix as a single configuration (m1,r,m2)(m_1,r,m_2) with MGN and RWP as special cases, and contracts the error across the entire spectrum, with its largest margins on the engineering quantities and at the smallest training budgets, where CAE surrogates actually operate. On public benchmarks it is the most accurate model in nearly every comparison, and on industry-scale meshes it attains the lowest error at every training budget. Limitations Our analysis establishes consistency rather than causal isolation, since the compounding factors behind the low-pass behavior are not disentangled by our experiments, and narrowing them down requires controlled interventions on the loss, the attention, and the bottleneck individually. The behavior also attenuates as the training set grows, so part of it belongs to the data regime rather than to the architecture alone, and separating the two requires sweeps beyond the budgets we report. The spectral analysis rests on a single steady RANS benchmark, and our datasets are steady throughout, so transient rollouts stress failure modes we do not measure. The multigrid reading is an analogy of roles rather than a convergence theorem. Finally, our comparisons with hybrid models that combine message passing with attention are limited to our own implementation of the MeshTransolver-style sandwich, whose reference code is not public, and a broader comparison with this class remains open. Acknowledgements The authors acknowledge the contributions of Diego Andrade, Jau-Uei Chen, and Oliver Littlewood for industrially-motivated simulation dataset preparation. References References [1] C. Adams, R. Ranade, S. Nidhan, M. A. Nabian, D. Akhare, R. Cherukuri, and S. Choudhry (2025) GeoTransolver: learning physics on irregular domains using multi-scale geometry aware physics attention transformer. External Links: 2512.20399 Cited by: Appendix F, §2. [2] B. Alkin, A. Fürst, S. Schmid, L. Gruber, M. Holzleitner, and J. Brandstetter (2024) Universal physics transformers: a framework for efficiently scaling neural operators. In Advances in Neural Information Processing Systems, Cited by: §2, §5.1. [3] U. Alon and E. Yahav (2021) On the bottleneck of graph neural networks and its practical implications. In International Conference on Learning Representations, Cited by: §1, §2. [4] N. Ashton, D. C. Maddix, S. Gundry, and P. M. Shabestari (2024) AhmedML: high-fidelity computational fluid dynamics dataset for incompressible, low-speed bluff body aerodynamics. External Links: 2407.20801 Cited by: Appendix F, §1, §5.6. [5] N. Ashton, C. Mockett, M. Fuchs, L. Fliessbach, H. Hetmann, T. Knacke, N. Schonwald, V. Skaperdas, G. Fotiadis, A. Walle, B. Hupertz, and D. Maddix (2024) DrivAerML: high-fidelity computational fluid dynamics dataset for road-car external aerodynamics. External Links: 2408.11969 Cited by: §1. [6] G. K. Batchelor (1951) Pressure fluctuations in isotropic turbulence. Mathematical Proceedings of the Cambridge Philosophical Society 47 (2), p. 359–374. Cited by: Appendix C. [7] F. Bonnet, J. A. Mazari, P. Cinnella, and P. Gallinari (2022) AirfRANS: high fidelity computational fluid dynamics dataset for approximating Reynolds-Averaged Navier–Stokes solutions. In Advances in Neural Information Processing Systems, Cited by: Appendix C, Appendix F, §1, §5.6. [8] J. Brandstetter, D. E. Worrall, and M. Welling (2022) Message passing neural PDE solvers. In International Conference on Learning Representations, Cited by: §5.3. [9] A. Brandt (1977) Multi-level adaptive solutions to boundary-value problems. Mathematics of Computation 31 (138), p. 333–390. Cited by: §1. [10] W. L. Briggs, V. E. Henson, and S. F. McCormick (2000) A multigrid tutorial. 2nd edition, SIAM. Cited by: §1. [11] C. Cai, T. S. Hy, R. Yu, and Y. Wang (2023) On the connection between MPNN and graph transformer. In Proceedings of the 40th International Conference on Machine Learning, PMLR. Cited by: §2. [12] Y. Cao, M. Chai, M. Li, and C. Jiang (2023) Efficient learning of mesh-based physical simulation with bi-stride multi-scale graph neural network. In Proceedings of the 40th International Conference on Machine Learning, PMLR. Cited by: §2. [13] G. Curtosi, C. M. Ruiz Ruiz, F. Cavaliere, and X. Larráyoz Izcara (2026) Crash assessment via mesh-based graph neural networks and physics-aware attention. External Links: 2605.11784 Cited by: §2. [14] G. Duthé et al. (2025) Graph transformers for inverse physics: reconstructing flows around arbitrary 2D airfoils. External Links: 2501.17081 Cited by: §2. [15] M. Fortunato, T. Pfaff, P. Wirnsberger, A. Pritzel, and P. Battaglia (2022) MultiScale meshgraphnets. In 2nd AI for Science Workshop at the 39th International Conference on Machine Learning, Cited by: §2. [16] J. Gilmer, S. S. Schoenholz, P. F. Riley, O. Vinyals, and G. E. Dahl (2017) Neural message passing for quantum chemistry. In Proceedings of the 34th International Conference on Machine Learning, PMLR. Cited by: §2. [17] R. J. Gladstone, H. Rahmani, V. Suryakumar, H. Meidani, M. D’Elia, and A. Zareei (2024) Mesh-based GNN surrogates for time-independent PDEs. Scientific Reports 14, p. 3394. Cited by: §2. [18] T. Gotoh and D. Fukayama (2001) Pressure spectrum in homogeneous turbulence. Physical Review Letters 86 (17), p. 3775–3778. Cited by: Appendix C. [19] M. M. Iparraguirre, I. Alfaro, D. Gonzalez, and E. Cueto (2026) MeshGraphNet-transformer: scalable mesh-based learned simulation for solid mechanics. External Links: 2601.23177 Cited by: §2. [20] A. Jabri, D. J. Fleet, and T. Chen (2023) Scalable adaptive computation for iterative generation. In Proceedings of the 40th International Conference on Machine Learning, PMLR. Cited by: §3. [21] A. Jaegle, S. Borgeaud, J. Alayrac, C. Doersch, C. Ionescu, D. Ding, S. Koppula, D. Zoran, A. Brock, E. Shelhamer, et al. (2022) Perceiver IO: a general architecture for structured inputs and outputs. In International Conference on Learning Representations, Cited by: §1, §2, §5.1, §5.2. [22] A. Jaegle, F. Gimeno, A. Brock, A. Zisserman, O. Vinyals, and J. Carreira (2021) Perceiver: general perception with iterative attention. In Proceedings of the 38th International Conference on Machine Learning, PMLR. Cited by: §1. [23] S. Janny, A. Beneteau, M. Nadri, J. Digne, N. Thome, and C. Wolf (2023) EAGLE: large-scale learning of turbulent fluid dynamics with mesh transformers. In International Conference on Learning Representations, Cited by: §2. [24] S. Khodakarami, V. Oommen, A. Bora, and G. E. Karniadakis (2025) Mitigating spectral bias in neural operators via high-frequency scaling for physical systems. External Links: 2503.13695 Cited by: §5.2. [25] R. Kress (2014) Linear integral equations. 3rd edition, Springer. Cited by: §3. [26] D. Kreuzer, D. Beaini, W. L. Hamilton, V. Létourneau, and P. Tossou (2021) Rethinking graph transformers with spectral attention. In Advances in Neural Information Processing Systems, Cited by: §2. [27] J. Lee, Y. Lee, J. Kim, A. R. Kosiorek, S. Choi, and Y. W. Teh (2019) Set transformer: a framework for attention-based permutation-invariant neural networks. In Proceedings of the 36th International Conference on Machine Learning, PMLR. Cited by: §5.2. [28] 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. Cited by: Appendix F, §5.6. [29] 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, Cited by: Appendix F, §5.6. [30] Z. Li, N. B. Kovachki, C. Choy, B. Li, J. Kossaifi, S. P. 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, Cited by: §2. [31] P. Lippe, B. S. Veeling, P. Perdikaris, R. E. Turner, and J. Brandstetter (2023) PDE-refiner: achieving accurate long rollouts with neural PDE solvers. In Advances in Neural Information Processing Systems, Cited by: §5.2. [32] H. Luo, H. Wu, H. Zhou, L. Xing, Y. Di, J. Wang, and M. Long (2025) Transolver++: an accurate neural solver for PDEs on million-scale geometries. In Proceedings of the 42nd International Conference on Machine Learning, PMLR. Cited by: §2. [33] N. Park and S. Kim (2022) How do vision transformers work?. In International Conference on Learning Representations, Cited by: §5.2. [34] T. Pfaff, M. Fortunato, A. Sanchez-Gonzalez, and P. Battaglia (2021) Learning mesh-based simulation with graph networks. In International Conference on Learning Representations, Cited by: §1, §2, §3, §3. [35] S. B. Pope (2000) Turbulent flows. Cambridge University Press. Cited by: Appendix C. [36] V. Puri, A. Joglekar, S. D. G. Bandreddi, K. Ferguson, Y. Chen, Y. J. Zhang, and L. B. Kara (2025) FLARE: fast low-rank attention routing engine. External Links: 2508.12594 Cited by: §1, §2, §5.1, §5.1. [37] S. Qin, F. Lyu, W. Peng, D. Geng, J. Wang, X. Tang, S. Leroyer, N. Gao, X. Liu, and L. Wang (2024) Toward a better understanding of Fourier neural operators from a spectral perspective. External Links: 2404.07200 Cited by: §5.2. [38] N. Rahaman, A. Baratin, D. Arpit, F. Draxler, M. Lin, F. A. Hamprecht, Y. Bengio, and A. Courville (2019) On the spectral bias of neural networks. In Proceedings of the 36th International Conference on Machine Learning, PMLR. Cited by: §5.2. [39] M. Ramezankhani, J. M. Patel, A. Deodhar, and D. Birru (2025) GITO: graph-informed transformer operator for learning complex partial differential equations. External Links: 2506.13906 Cited by: §2. [40] L. Rampášek, M. Galkin, V. P. Dwivedi, A. T. Luu, G. Wolf, and D. Beaini (2022) Recipe for a general, powerful, scalable graph transformer. In Advances in Neural Information Processing Systems, Cited by: §2. [41] W. Ripken, L. Coiffard, F. Pieper, and S. Dziadzio (2023) Multiscale neural operators for solving time-independent PDEs. In The Symbiosis of Deep Learning and Differential Equations I Workshop at NeurIPS, Cited by: §2. [42] L. Serrano, T. X. Wang, E. Le Naour, J. Vittaut, and P. Gallinari (2024) AROMA: preserving spatial structure for latent PDE modeling with local neural fields. In Advances in Neural Information Processing Systems, Cited by: §2, §5.1. [43] H. Shirzad, A. Velingker, B. Venkatachalam, D. J. Sutherland, and A. K. Sinop (2023) Exphormer: sparse transformers for graphs. In Proceedings of the 40th International Conference on Machine Learning, PMLR. Cited by: §2. [44] M. Takamoto, T. Praditia, R. Leiteritz, D. MacKinlay, F. Alesiani, D. Pflüger, and M. Niepert (2022) PDEBench: an extensive benchmark for scientific machine learning. In Advances in Neural Information Processing Systems, Cited by: Appendix B, §1, §5. [45] M. Tancik, P. P. Srinivasan, B. Mildenhall, S. Fridovich-Keil, N. Raghavan, U. Singhal, R. Ramamoorthi, J. T. Barron, and R. Ng (2020) Fourier features let networks learn high frequency functions in low dimensional domains. In Advances in Neural Information Processing Systems, Cited by: §3. [46] J. Topping, F. Di Giovanni, B. P. Chamberlain, X. Dong, and M. M. Bronstein (2022) Understanding over-squashing and bottlenecks on graphs via curvature. In International Conference on Learning Representations, Cited by: §2. [47] Y. H. Tsai, S. Bai, M. Yamada, L. Morency, and R. Salakhutdinov (2019) Transformer dissection: a unified understanding of transformer’s attention via the lens of kernel. In Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing, Cited by: §5.2. [48] Y. Tsuji and T. Ishihara (2003) Similarity scaling of pressure fluctuation in turbulence. Physical Review E 68, p. 026309. Cited by: Appendix C. [49] D. G. Vlaykov and M. Wilczek (2019) On the small-scale structure of turbulence and its impact on the pressure field. Journal of Fluid Mechanics 861, p. 422–446. Cited by: Appendix C. [50] H. Wang, J. Li, A. Dwivedi, K. Hara, and T. Wu (2024) BENO: boundary-embedded neural operators for elliptic PDEs. In International Conference on Learning Representations, Cited by: §2, §3. [51] P. Wang, W. Zheng, T. Chen, and Z. Wang (2022) Anti-oversmoothing in deep vision transformers via the Fourier domain analysis: from theory to practice. In International Conference on Learning Representations, Cited by: §5.2. [52] T. Wang and C. Wang (2024) Latent neural operator for solving forward and inverse PDE problems. In Advances in Neural Information Processing Systems, Cited by: §2, §5.1. [53] S. Wen S. Mousavi et al. (2025) Geometry aware operator transformer as an efficient and accurate neural surrogate for PDEs on arbitrary domains. External Links: 2505.18781 Cited by: §2. [54] 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, PMLR. Cited by: Appendix F, §1, §2, §5.1. [55] Z. Wu, P. Jain, M. A. Wright, A. Mirhoseini, J. E. Gonzalez, and I. Stoica (2021) Representing long-range context for graph neural networks with global attention. In Advances in Neural Information Processing Systems, Cited by: §2. [56] Z. J. Xu, Y. Zhang, T. Luo, Y. Xiao, and Z. Ma (2020) Frequency principle: Fourier analysis sheds light on deep neural networks. Communications in Computational Physics 28 (5), p. 1746–1767. Cited by: §5.2. [57] Z. Yang, H. Xin, T. Du, and L. Liu (2026) Simple yet effective: low-rank spatial attention for neural operators. External Links: 2604.03582 Cited by: §5.2. [58] C. Ying, T. Cai, S. Luo, S. Zheng, G. Ke, D. He, Y. Shen, and T. Liu (2021) Do transformers really perform badly for graph representation?. In Advances in Neural Information Processing Systems, Cited by: §2. Appendix A Datasets and Simulations Both environments are internal-flow cases solved with the Finite Volume Method (FVM), using a RANS turbulence model, a velocity-inlet and pressure-outlet boundary-condition pair, and no-slip walls. The branched pipe applies a one-seventh power-law velocity profile at the inlet; the pump applies a uniform inlet velocity and models the impeller rotation with the Moving Reference Frame (MRF) method. Both are steady-state simulations solved with an adaptive time-stepping strategy. Table 5 summarizes the simulation configurations, Table 6 the design parameters sampled to generate each dataset, and Table 7 the modeling task and dataset specification shared by all models. Branched pipe Centrifugal pump (casing) (full) cell number ∼60 60k ∼270 270k ∼1.4 1.4M fluid air water water turbulence model RANS RANS RANS boundary conditions velocity inlet, pressure outlet, no-slip wall solver FVM (Fluent) special treatment — MRF MRF Table 5: Summary of simulation configurations. Parameter Range Branched pipe angle of middle baffle plate 0–60 (deg) aperture length of inlet baffle plate 20–37 (m) angle between ducts 75–120 (deg) curvature angle of duct branch 90–180 (deg) inlet mean velocity 3–10 (m/s) Centrifugal pump inlet radius 160–200 (m) impeller radius 80–100 (m) blade pitch angle 5–40 (deg) inlet flow velocity 1.5–3.5 (m/s) Table 6: Design parameters used to generate the simulation datasets, with the range each parameter is sampled over. Branched pipe Centrifugal pump modeling task steady-state steady-state total samples 300 300 training-set sizes 50 / 100 / 300 samples train/val/test split 0.7/0.15/0.15 random seed 42 field normalization per-component standardization input features coordinates, inlet/outlet boundary conditions input / processing domain full mesh modeled fields p, uxu_x, uzu_z uxu_x, uyu_y, uzu_z, p Table 7: Modeling task and dataset specification, common to all models. The two pump extents share the same design of experiments and splits. Appendix B Metric Definitions All metrics compare a predicted field u against the reference solution u∗u^* through the error e=u−u∗e=u-u^*, in physical units after inverting the training normalization. Cell-centered data are converted to point data by volume-weighted averaging before all point-based metrics, with R2R^2 the one exception, and all means are unweighted. The analysis metrics are computed per case, field, and model on the held-out test split, capped to eight cases, and reported as means over these cases; the benchmark tables instead use the full 45-case test split. Truth and prediction always pass through the same operators, so operator error largely cancels in the relative metrics. Relative L2L_2 mean(e2)/mean(u∗2) mean(e^2)/mean(u^*2) over all volume mesh points. R2R^2 1−mean(e2)/var(u∗)1-mean(e^2)/var(u^*), computed on the raw cell-centered arrays. Gradients Spatial gradients use a k-nearest-neighbor weighted least-squares operator rather than a mesh-based filter. For each point, the k=20k=20 nearest neighbors define displacements djd_j and value differences δvjδ v_j, and the gradient g solves the ridge-regularized weighted normal equations (∑jwjdjdj⊤+ϵI)g=∑jwjdjδvj ( _jw_jd_jd_j +ε I )\,g= _jw_jd_j\,δ v_j with wj=1/|dj|w_j=1/|d_j|. The operator reproduces linear fields to machine precision, and the same operator is applied to truth and prediction. The gradient relative L2L_2 is mean|∇e|2/mean|∇u∗|2 mean|∇ e|^2/mean|∇ u^*|^2 with |⋅||·| the Euclidean norm. Error spectrum Spectral quantities live on a two-dimensional mid-plane slice. A y-normal plane at the domain’s mid-y is sampled on a uniform raster 512 pixels wide with square pixels, points outside the fluid are masked and zero-filled, and the two-dimensional FFT power is binned into integer radial wavenumber rings k=round(kx2+kz2)k=round ( k_x^2+k_z^2 ), in cycles per domain extent. Eerr(k)E_err(k) is the ring-summed power spectrum of the error field and Etrue(k)E_true(k) that of the reference; wavenumbers above k=92k=92, where raster-interpolation noise floors the spectrum, are excluded, as is the DC ring. Band-resolved fRMSE fRMSE, the RMSE in Fourier space, is the band-relative error ∑k∈bandEerr(k)/∑k∈bandEtrue(k) _k E_err(k)/ _k E_true(k), 0 for a perfect prediction and 1 when the error is as energetic as the truth within the band. The cutoffs follow the PDEBench bins [44], 0≤k≤40≤ k≤ 4 (low), 5≤k≤125≤ k≤ 12 (mid), and 13≤k≤9213≤ k≤ 92 (high). Vorticity ωy=∂ux/∂z−∂uz/∂x _y=∂ u_x/∂ z-∂ u_z/∂ x, computed with the same gradient operator and scored with the relative L2L_2 formula as a scalar field. Wall shear stress The reported components are the wall-normal derivative proxies τw,i=μ(∇ui⋅n^) _w,i=μ\,(∇ u_i· n) for i∈x,zi∈\x,z\, with μ the constant laminar viscosity and n n the unit wall normals from the solver’s face-area vectors; volume gradients are interpolated to each wall point as the mean over its four nearest volume points. This proxy is the dominant term of the wall shear stress for these components, and its truth values track the solver’s stored wall shear. Scores are the relative L2L_2 and R2R^2 over all wall boundary points, and the parity plots of Fig. 7 show predicted against true μ∂ui/∂nμ\,∂ u_i/∂ n pooled over cases. Error versus wall distance Wall distance dwd_w is the solver’s stored wall-distance field. Points fall into twelve quantile bins of dwd_w, once over the full range and once over the nearest-to-wall fifth of points, and each bin reports the mean absolute error, averaged over the test cases. NMAE Mean absolute error over all mesh cells, normalized by the standard deviation of the ground-truth field. Training-dynamics normalization The training-dynamics figures evaluate five evenly spaced checkpoints on four cases. Every curve is divided by RWP’s first-checkpoint value, per field and wavenumber or band, so a value of 1 means as wrong as RWP at initialization. Appendix C Pressure and the Low-Pass Bias Across our experiments, pressure consistently behaves as the easiest field for the latent-attention models. In Table 8 the strongest latent-attention baseline matches or exceeds RWR on every pressure metric, and its band-resolved errors rise only mildly from the low to the high band on pressure while they steepen sharply on the velocity components. The same pattern has been reported for steady RANS airfoil flows, where surrogate models across architectures inferred surface pressure and lift coefficient far more accurately than wall shear stress and drag coefficient [7]. The behavior has a physical reading. In incompressible flow, taking the divergence of the momentum equation yields a Poisson equation for the pressure, ∇2p=−ρ(∂ui/∂xj)(∂uj/∂xi)∇^2p=-ρ\,(∂ u_i/∂ x_j)(∂ u_j/∂ x_i), whose Green’s-function solution expresses the pressure at a point as a domain-wide, 1/|x−y|1/|x-y|-weighted integral of the velocity-gradient source [35]. The inverse Laplacian attenuates the high-wavenumber content of this source as 1/k21/k^2, and pressure is accordingly characterized as a long-range field [49]. Consistently, in high-Reynolds-number turbulence the inertial-range pressure spectrum falls as k−7/3k^-7/3, steeper than the k−5/3k^-5/3 velocity spectrum [6, 18], although the exponent is only cleanly observed at high Reynolds number [48] and inertial-range scalings do not strictly apply to the steady RANS mean fields of our benchmark. The elliptic-smoothing mechanism, however, is Reynolds-independent, since the divergence of the steady RANS momentum equation yields the same Poisson structure for the mean pressure, with a source built from mean-velocity gradients and Reynolds-stress divergences. A model biased toward low wavenumbers therefore degrades pressure least and the velocity components most, whose spectra keep more of their content in the bands the bottleneck sheds. Appendix D Additional Diagnostics Figure 5: Band-resolved test error of RWP as a function of the number of latents L on the branched pipe. Each panel shows the fRMSE over the low, mid, and high wavenumber bands for pressure and the two velocity components, with L swept from 8 to 512 and all other settings fixed. Figure 6: Band-resolved training dynamics on the branched pipe. Each panel tracks the fRMSE of the low, mid, and high wavenumber bands over training for the two velocity components (rows) under MGN, RWP, and RWR (columns), normalized by the initial band value of RWP. D.1 Full Branched Pipe Metrics Table 8 reports the complete per-field metrics behind Table 2, the bulk and gradient relative L2L_2 errors together with the band-resolved fRMSE over the low, mid, and high wavenumber bands, for pressure and both velocity components. rel. L2L_2 fRMSE Field Model field grad. low mid high p Transolver 0.090 [0.85] 0.165 [0.95] 0.118 [0.99] 0.129 [0.94] 0.133 [0.87] GeoTransolver 0.135 [1.28] 0.212 [1.22] 0.159 [1.33] 0.168 [1.22] 0.174 [1.14] RWP 0.134 [1.27] 0.199 [1.14] 0.155 [1.30] 0.177 [1.29] 0.191 [1.25] RWR 0.106 0.174 0.119 0.137 0.153 uxu_x Transolver 0.078 [1.01] 0.174 [1.18] 0.035 [0.95] 0.070 [1.10] 0.129 [1.37] GeoTransolver 0.098 [1.28] 0.226 [1.54] 0.057 [1.54] 0.073 [1.15] 0.107 [1.13] RWP 0.090 [1.17] 0.168 [1.14] 0.052 [1.39] 0.090 [1.42] 0.149 [1.59] RWR 0.077 0.147 0.037 0.063 0.094 uzu_z Transolver 0.080 [1.11] 0.207 [1.27] 0.060 [0.72] 0.086 [0.94] 0.150 [1.20] GeoTransolver 0.104 [1.45] 0.286 [1.75] 0.072 [0.87] 0.092 [1.01] 0.136 [1.08] RWP 0.091 [1.27] 0.218 [1.33] 0.092 [1.11] 0.121 [1.32] 0.199 [1.59] RWR 0.072 0.164 0.083 0.092 0.125 Table 8: Held-out test errors on the branched pipe at matched training budget for all output fields: relative L2L_2 error of each field and of its first spatial gradients, and band-resolved fRMSE over the low, mid, and high wavenumber bands. Bracketed values give the ratio of each model’s error to RWR’s; bold marks the best value per metric and field. The pressure results are discussed in App. C. D.2 Composition Ablation A stacked composition that runs all relaxation sweeps before or after the latent-attention iterations contains the same two mechanisms as an interleaved one. Table 9 compares four allocations of an identical budget of eight relaxation sweeps and four latent-attention iterations, from the two stacked orderings to progressively finer interleavings. Steps p uxu_x uzu_z (8,4,0)×1(8,4,0)×1 5.8 5.8 4.8 (0,4,8)×1(0,4,8)×1 6.0 5.3 4.7 (2,2,2)×2(2,2,2)×2 5.2 5.2 4.5 (1,1,1)×4(1,1,1)×4 5.1 5.1 4.5 Table 9: Composition ablation on the branched pipe at a fixed budget of eight relaxation sweeps and four latent iterations: held-out test normalized mean absolute error (NMAE, ×10−2× 10^-2) of pressure and the two velocity components. A configuration (m1,r,m2)×nB(m_1,r,m_2)×n_B runs m1m_1 relaxation sweeps, r latent iterations, and m2m_2 further sweeps per block, repeated over nBn_B blocks: (8,4,0)(8,4,0) places all relaxation before the latent path, (0,4,8)(0,4,8) after it, and the remaining rows interleave the two at increasing granularity. D.3 Wall Diagnostics Fig. 7 shows parity plots of the wall-normal velocity derivatives over wall cells, and Fig. 8 the mean absolute error of each field as a function of wall distance, with a near-wall zoom. Figure 7: Parity plots of the wall-normal derivatives ∂ux/∂n∂ u_x/∂ n and ∂uz/∂n∂ u_z/∂ n over wall cells on held-out branched pipe cases, with the coefficient of determination R2R^2 per model. Figure 8: Mean absolute error of pressure and the velocity components as a function of wall distance dwd_w on held-out branched pipe cases (top row), with a near-wall zoom (bottom row). D.4 Analysis Training Details All analysis models are trained on the branched pipe in the same pipeline, under the 100-sample branched pipe training protocol of the main benchmarks (App. E). MGN, which is absent from the benchmarks, trains with an initial learning rate of 3×10−43× 10^-4 for the same 550 epochs as RWP and RWR. Table 10 lists the architectural configuration of RWP and RWR. The latent queries are geometry-anchored with learned anchor coordinates, passed through the shared random-Fourier-feature embedding (σ=0.5σ=0.5, learnable coefficients, dimension equal to the hidden width), and the boundary-only read is enabled in all runs. MGN runs 15 message-passing sweeps at hidden width 128. Transolver is our reimplementation with 8 layers, hidden width 172, 4 attention heads, 32 slices, and MLP ratio 2; GeoTransolver is our reimplementation with 4 layers, hidden width 144, 8 heads, 32 slices, and ball-query radii 0.25,1.0\0.25,1.0\ on normalized coordinates with 8,32\8,32\ neighbors. RWP RWR steps (m1,r,m2)×nB(m_1,r,m_2)× n_B (0,2,0)×2(0,2,0)× 2 (2,2,2)×2(2,2,2)× 2 latents L 128 128 hidden width 192 128 attention heads 4 4 latent self-attention none none boundary-only read yes yes weight sharing across r no no Table 10: Architectural configuration of RWP and RWR for the branched pipe analysis. Hidden widths are chosen to approximately match parameter counts across MGN, RWP, and RWR. Appendix E Main Benchmark Results Table 12 reports the per-field NMAE, parameter counts, peak memory, and training cost behind Table 4. Training protocol All models are trained using the same pipeline under the same training protocol. Each dataset contains 300300 simulations, split 70/15/1570/15/15 into train/validation/test (210/45/45210/45/45 samples). The data-limited variants use only the first 5050 or 100100 samples, split the same way (35/7/835/7/8 and 70/15/1570/15/15). Inputs and outputs are cell-centered mesh fields for the branched pipe data set and node-centered for the pump data sets. Input fields are standardized per component, with node coordinates appended. We minimize the mean-squared error with AdamW (batch size 11, gradient-norm clipping at 2.02.0) under a linear schedule that decays the learning rate to one-tenth of its base value over the first (nepoch−100n_epoch-100) epochs and holds it constant thereafter. The base learning rate is 10−310^-3 for the branched pipe and 3×10−43× 10^-4 for the pump. The training budget on the branched pipe is 550550 epochs for RWP and RWR and 650650 epochs for the attention baselines. For the pump-casing- and full-pump experiments all models are trained for 600600 and 500500 epochs respectively. For all completed runs we report metrics computed on the test set. Each run uses a single A100 80 GB GPU with gradient checkpointing enabled. Model configurations RWR and RWP share a MeshGraphNet backbone (two blocks with hidden width 128128 on the branched pipe and casing, and width 9696 on the full pump) wrapping a latent-attention bottleneck: RWR interleaves message-passing sweeps with the latent read/write (steps=[2,2,2] steps=[2,2,2]), while RWP performs the read/write alone (steps=[0,2,0] steps=[0,2,0]). Transolver consists of eight layers with hidden width 256256 and uses eight heads. GeoTransolver consists of four layers with width 256256, eight heads, and local geometric features gathered within radii 0.25,1.0\0.25,1.0\. Exact parameter counts per configuration are listed in Table 12. Sweeps and selection For every model family we sweep its principal capacity- and locality/geometry-governing hyperparameters and report the best model variant in Table 11 and their corresponding performance metrics in Table 12. For RWR and RWP we vary the processor depth (two or three blocks), the latent-token count L, and whether the read attends to the full volume or only the boundary surface. Transolver and GeoTransolver are swept over their slice count and, in the data-limited setting, their hidden width. Table 11: Architectural details of best performing models per family for each mesh and data setting. Transolver GeoTransolver RWR RWP #samples #params #slices width #params #slices width #params #latents blocks read width #params #latents blocks read width Branched pipe (∼ 60k cells) 50 3.58M 32 172 4.02M 32 144 3.58M 064 2 bdy 128 3.87M 512 2 bdy 128 100 3.58M 32 172 4.02M 16 144 2.79M 128 2 bdy 128 5.65M 128 3 bdy 128 300 7.72M 32 256 9.53M 32 256 5.30M 128 3 bdy 128 2.52M 128 3 bdy 128 Pump casing (∼ 270k cells) 50 7.89M 32 256 9.53M 32 256 5.30M 128 3 bdy 128 5.65M 128 3 bdy 128 100 7.89M 32 256 9.53M 32 256 5.30M 128 3 bdy 128 5.65M 128 3 bdy 128 300 7.72M 32 256 9.52M 16 256 5.30M 128 3 bdy 128 2.52M 128 3 bdy 128 Pump full (∼1.4 1.4M cells) 50 1.10M 32 096 2.46M 32 096 2.03M 512 2 bdy 128 1.73M 128 2 bdy 128 100 1.10M 32 096 2.46M 32 096 2.02M 128 2 bdy 128 1.73M 128 2 bdy 128 300 1.10M 32 096 1.99M 32 080 2.03M 512 2 bdy 096 1.73M 128 2 all 128 Table 12: Performance metrics of best models per family for each mesh and data setting. We report the NMAE per output field, across-field mean NMAE, parameter count, peak memory, average wall-clock time per epoch, and the number of epochs to reach a mean R2≥0.95R^2≥ 0.95 (n/r: not reached). The branched pipe is quasi-2D, so uyu_y is absent. NMAE (×10−2× 10^-2) Model p uxu_x uyu_y uzu_z mean #params peak memory average time/epoch #epochs until avg. R2≥0.95R^2≥ 0.95 Branched pipe (∼ 60k cells) 50 samples Transolver 10.3 12.7 10.8 11.3 3.58M 08.1 GB 8s 470 GeoTransolver 9.5 11.8 10.2 10.5 4.02M 08.0 GB 7s 370 RWP 10.4 10.7 9.7 10.3 3.87M 04.2 GB 6s 290 RWR 8.7 8.9 7.7 8.4 3.58M 13.5 GB 11s 160 100 samples Transolver 8.5 7.1 7.0 7.5 3.58M 06.9 GB 16s 220 GeoTransolver 9.6 7.8 7.5 8.3 4.02M 07.4 GB 14s 190 RWP 8.7 6.7 6.4 7.3 5.65M 04.6 GB 15s 140 RWR 8.7 6.3 5.8 6.9 2.79M 10.1 GB 23s 110 300 samples Transolver 4.8 5.3 4.5 4.8 7.72M 13.5 GB 53s 110 GeoTransolver 5.9 5.6 5.0 5.5 9.53M 15.6 GB 51s 140 RWP 5.0 5.4 4.5 5.0 2.52M 03.6 GB 40s 050 RWR 5.0 5.0 4.3 4.7 5.30M 26.7 GB 95s 040 Pump casing (∼ 270k cells) 50 samples Transolver 06.4 12.0 30.9 13.0 15.6 7.89M 28.5 GB 44s 230 GeoTransolver 04.6 9.3 24.5 10.0 12.1 9.53M 31.8 GB 43s 160 RWP 03.9 8.8 23.5 9.5 11.4 5.65M 14.2 GB 30s 080 RWR 04.3 8.8 21.0 9.1 10.8 5.30M 21.8 GB 57s 070 100 samples Transolver 03.9 8.2 23.1 8.3 10.9 7.89M 28.6 GB 87s 120 GeoTransolver 03.1 7.0 20.0 7.1 9.3 9.53M 31.8 GB 85s 080 RWP 02.6 6.3 17.2 6.3 8.1 5.65M 14.1 GB 59s 030 RWR 02.6 6.1 15.5 6.3 7.7 5.30M 21.9 GB 111s 040 300 samples Transolver 02.5 6.0 16.8 6.1 7.9 7.72M 31.4 GB 265s 040 GeoTransolver 02.2 5.3 14.7 5.3 6.9 9.52M 29.6 GB 241s 040 RWP 02.1 5.1 13.0 5.2 6.3 2.52M 09.9 GB 156s 020 RWR 01.7 4.5 10.9 4.5 5.4 5.30M 22.0 GB 331s 020 Pump full (∼1.4 1.4M cells) 50 samples Transolver 33.0 27.4 62.8 25.1 37.1 1.10M 75.2 GB 170s n/r GeoTransolver 26.8 24.7 54.9 22.4 32.2 2.46M 79.2 GB 148s n/r RWP 09.1 19.8 40.2 18.0 21.8 1.73M 33.7 GB 91s n/r RWR 07.8 17.3 36.0 15.9 19.3 2.03M 59.4 GB 182s n/r 100 samples Transolver 09.1 21.3 43.2 19.1 23.2 1.10M 75.5 GB 331s n/r GeoTransolver — — — — — 2.46M 79.2 GB 285s n/r RWP 06.8 15.9 31.7 15.0 17.3 1.73M 34.3 GB 177s 200 RWR 06.4 13.4 28.5 12.5 15.2 2.02M 59.1 GB 340s 140 300 samples Transolver 05.9 14.3 30.1 13.4 15.9 1.10M 75.7 GB 986s 160 GeoTransolver 05.3 12.3 27.6 11.4 14.1 1.99M 73.8 GB 830s 170 RWP 05.4 12.7 27.5 11.9 14.4 1.73M 44.9 GB 968s 100 RWR 05.7 10.8 24.2 9.9 12.7 2.03M 59.6 GB 985s 060 Appendix F Additional Benchmarks Datasets AirfRANS: two-dimensional RANS airfoil flows with roughly 180k nodes per case, velocity, pressure, and turbulent viscosity as targets, and per-case inlet velocity and angle of attack as global parameters [7]. AhmedML: three-dimensional surface meshes of roughly 1.1M points over 500 Ahmed-body variants, surface pressure and wall shear stress as targets, eight global geometry parameters, and a 400/50/50 split [4]. Geo-FNO airfoil and pipe [28] and Darcy flow [29], with the standard splits of roughly 1000 training and 200 test cases used by published Transolver-family reproductions. Input parity Every model receives the same per-case information. The MGN-family models encode geometry through relative edge displacements, RWR additionally through the positional embeddings of its global stage, and Transolver and GeoTransolver receive coordinates concatenated per token. Per-case global parameters are injected into every message-passing step for the MGN-family models, broadcast-concatenated per token for RWP and Transolver, and passed through the dedicated global-context tokenizer for GeoTransolver. The read attends to all nodes on every dataset; on AhmedML the boundary-only read is enabled but coincides with the full read, since every node lies on the surface. Model Hyperparameters RWR uses 3 blocks on AirfRANS and AhmedML and 2 blocks on the Geo-FNO datasets; RWP uses 4 and 3 blocks, respectively. For the read-write stage, we use 128 latents on AirfRANS and AhmedML and 64 latents on the Geo-FNO datasets, at a hidden width of 128 and with 4 attention heads throughout. Transolver and GeoTransolver are run with 8 layers, hidden width 256, MLP ratio 2, and 32 slices, and 4 layers with hidden width 256 and 32 slices, respectively [54, 1]. GeoTransolver’s multi-scale ball-query mechanism supports only three-dimensional coordinates and is active on the AhmedML surface dataset; on the two-dimensional datasets it runs without it. Our implementation of MeshTransolver uses a 2–2–2 structure of MGN–Transolver–MGN layers, with the two Transolver layers using 4 attention heads, 32 slices, and a hidden width of 128 on all datasets. All MGN-based models use the same 2-layer MLP layer widths of [128,128] on all datasets. Training protocol AdamW with a base learning rate of 10−310^-3 and weight decay 0.01, linear learning-rate decay to 0.1×0.1× over the full run, and gradient-norm clipping at 2.0. Models train in bfloat16 on AirfRANS and AhmedML and in fp32 on the Geo-FNO datasets, with per-block activation checkpointing where required to fit memory. The batch size is 1 on AirfRANS and AhmedML; on the Geo-FNO datasets it is 8 for the MGN-family models and 1 for Transolver and GeoTransolver. Training runs 100 epochs on AirfRANS and AhmedML and 800 epochs on the Geo-FNO datasets; on the Geo-FNO airfoil each model reports the better of two settings, (10−3CLOSE(10^-3, 800 epochs)) and (5×10−4CLOSE(5× 10^-4, 1000 epochs)). The loss is the MSE on standardized fields, each run reports its best validation loss, and all results are single seed. Appendix G Model Predictions Fig. 9–14 show model predictions and signed pointwise errors for two held-out test cases per dataset, for the pressure and the velocity magnitude. Branched pipe fields are shown on the mid-y slice, pump-casing fields on the y=0.04y=0.04 m plane, and full-pump fields on the impeller surface. Each model panel reports its R2R^2 and NMAE for the shown field and case. Figure 9: Predicted velocity magnitude on the branched pipe for two held-out test cases (top and bottom groups), shown on the mid-y slice. The left column shows the ground truth; each model column shows the prediction, with its R2R^2 and NMAE, above the signed pointwise error. Figure 10: Predicted pressure on the branched pipe for two held-out test cases (top and bottom groups), shown on the mid-y slice. The left column shows the ground truth; each model column shows the prediction, with its R2R^2 and NMAE, above the signed pointwise error. Figure 11: Predicted velocity magnitude on the pump casing for two held-out test cases (top and bottom groups), shown on the y=0.04y=0.04 m plane. The left column shows the ground truth; each model column shows the prediction, with its R2R^2 and NMAE, above the signed pointwise error. Figure 12: Predicted pressure on the pump casing for two held-out test cases (top and bottom groups), shown on the y=0.04y=0.04 m plane. The left column shows the ground truth; each model column shows the prediction, with its R2R^2 and NMAE, above the signed pointwise error. Figure 13: Predicted velocity magnitude on the full pump mesh for two held-out test cases (top and bottom groups), shown on the impeller surface. The left column shows the ground truth; each model column shows the prediction, with its R2R^2 and NMAE, above the signed pointwise error. Figure 14: Predicted pressure on the full pump mesh for two held-out test cases (top and bottom groups), shown on the impeller surface. The left column shows the ground truth; each model column shows the prediction, with its R2R^2 and NMAE, above the signed pointwise error.