Paper deep dive
Hybrid Lagrangian-Eulerian Model for Lagrangian Fluid Simulation
Ruoyan Li, Wei Wang, Yizhou Sun
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:Pure Lagrangian neural simulators offer geometric flexibility and exact advection, making them well-suited for modeling moving domains and free surfaces. However, the absence of a fixed global reference frame introduces two severe limitations: a spatial bottleneck, in which model capacity is wasted on uniform regions because the dense particle neighborhoods required for stable gradients are applied indiscriminately, and rapid temporal drift, caused by purely local message passing that lacks a global anchor. Inspired by classical hybrid numerical solvers, we propose a Hybrid Lagrangian-Eulerian neural simulator that augments Lagrangian dynamics with an Eulerian representation. To address the spatial bottleneck, we introduce adaptive downsampling that eliminates kinematic redundancy, preserving micro-scale details on particles while aggregating compressed features onto Eulerian nodes to resolve large-scale dynamics. To counter temporal drift, we employ a cross-attention mechanism that queries these Eulerian features, using the fixed grid as a stable spatial anchor to correct trajectory deviations at every timestep. Comprehensive experiments show that this hierarchical, cross-attended design substantially suppresses error accumulation, establishing a new state-of-the-art for accuracy and rollout stability in Lagrangian fluid simulation.
Tags
Links
- Source: https://arxiv.org/abs/2608.01164v1
- Canonical: https://arxiv.org/abs/2608.01164v1
Trouble viewing inline? Open PDF directly →
Full Text
91,044 characters extracted from source content.
Expand or collapse full text
Hybrid Lagrangian-Eulerian Model for Lagrangian Fluid Simulation Ruoyan Li University of California, Los Angeles Wei Wang University of California, Los Angeles &Yizhou Sun University of California, Los Angeles Abstract Pure Lagrangian neural simulators offer geometric flexibility and exact advection, making them well-suited for modeling moving domains and free surfaces. However, the absence of a fixed global reference frame introduces two severe limitations: a spatial bottleneck, in which model capacity is wasted on uniform regions because the dense particle neighborhoods required for stable gradients are applied indiscriminately, and rapid temporal drift, caused by purely local message passing that lacks a global anchor. Inspired by classical hybrid numerical solvers, we propose a Hybrid Lagrangian–Eulerian neural simulator that augments Lagrangian dynamics with an Eulerian representation. To address the spatial bottleneck, we introduce adaptive downsampling that eliminates kinematic redundancy, preserving micro-scale details on particles while aggregating compressed features onto Eulerian nodes to resolve large-scale dynamics. To counter temporal drift, we employ a cross-attention mechanism that queries these Eulerian features, using the fixed grid as a stable spatial anchor to correct trajectory deviations at every timestep. Comprehensive experiments show that this hierarchical, cross-attended design substantially suppresses error accumulation, establishing a new state-of-the-art for accuracy and rollout stability in Lagrangian fluid simulation. 1 Introduction Neural network–based physical simulators have recently achieved remarkable progress in fluid dynamics (Wang et al., 2024, 2025). These models are typically built on either an Eulerian formulation, where dynamics are resolved on fixed spatial domains, or a Lagrangian formulation, which models the fluid as a set of discrete particles that carry and update physical quantities along their moving trajectories. Lagrangian representations are highly favored for their physical intuition and geometric flexibility. By tracking discrete particles, these methods naturally handle complex phenomena like moving domains and free surfaces, which are notoriously difficult to capture on fixed grids (Monaghan, 1992; Gingold & Monaghan, 1977). Furthermore, since particles carry their properties as they move, advection is treated exactly and computational effort is concentrated where mass exists, avoiding the numerical diffusion and wasted resolution that plague grid-based schemes. For example, the dam break dataset in Section 5.1 and the stretching droplet dataset in Section 5.4 are particularly well-suited for a Lagrangian representation. An Eulerian approach would be highly inefficient here, as tracking the water splashes would require a prohibitively large and dense grid, and the stretching droplet would demand costly re-meshing at every time step. However, the pure Lagrangian neural simulator introduces fundamental modeling challenges, primarily stemming from the lack of a fixed global reference frame. This absence cascades into two intrinsic challenges. First, without a background grid, Lagrangian methods suffer from a severe spatial and modeling bottleneck. While an Eulerian grid can calculate stable finite differences with as few as two neighbors, particle-based methods rely on kernel averaging that often requires a dense neighborhood of 30–50 particles to yield stable gradients (Gingold & Monaghan, 1977). Consequently, maintaining grid-equivalent fidelity requires a drastically higher particle count (Agertz et al., 2007; Price, 2012; Dehnen & Aly, 2012). Because of this dense packing, local particle motions often exhibit indistinguishable kinematics, forcing the neural solver to expend significant representational capacity on parsing highly correlated, nearly uniform local dynamics rather than capturing the more complex interactions. Second, the lack of a global anchor causes rapid temporal error accumulation. Because particles update their states based solely on local message passing, any slight deviation in a trajectory causes incorrect neighborhood connections in the next timestep. Without a fixed Eulerian grid to anchor and correct these local mistakes, topological errors compound geometrically, causing the simulation to drift significantly. In classical computational physics, this precise dilemma led to the development of hybrid methods, most notably the Particle-in-Cell (PIC) (Francis Harvey Harlow, 1955; Dawson, 1983; Tskhakaya, 2008) method. PIC fundamentally improves Lagrangian simulation by leveraging an Eulerian grid as a global computational anchor. Instead of resolving complex PDE dynamics strictly between unstructured particle neighbors, PIC transfers particle properties onto a fixed background grid. The governing equations are then solved stably and efficiently on this Eulerian grid, and the resulting physical fields are interpolated back to update the particles. This approach anchors the drifting particle trajectories to a robust global structure while preserving the Lagrangian advantages. Inspired by this classical synergy, we investigate how Eulerian representations can be similarly leveraged to enhance the accuracy and robustness of Lagrangian neural simulators. We propose a Hybrid Lagrangian–Eulerian neural model that downsamples indistinguishable particles and aggregates their information onto a set of Eulerian nodes, where the governing equations are solved efficiently. We further employ a cross-attention mechanism that facilitates information exchange between the Lagrangian and Eulerian features. The downsampling operation mitigates the computational burden of high-density representations by eliminating redundancy among particles with indistinguishable kinematics. This architecture also eases modeling difficulty through a hierarchical decomposition, as the downsampler handles micro-scale details while Eulerian aggregation resolves the underlying PDE at a coarse level. Additionally, because Eulerian nodes never move and their features are tied to fixed spatial locations, the Eulerian branch acts as a spatial anchor. Cross-attention re-grounds the Lagrangian update at every step against this stable reference, suppressing the error accumulation that causes pure-Lagrangian surrogates to drift into unphysical particle distributions over long rollouts. Ultimately, our design exploits the Eulerian representation to enhance the accuracy and stability of Lagrangian neural simulators. Our contributions are as follows: (i) Problem Identification: We identify the fundamental limitations of pure Lagrangian neural simulators, specifically, spatial bottlenecks and temporal drift, and explain the benefits of a global Eulerian reference frame; (i) Practical Solution: Inspired by classical hybrid methods like PIC, we present a Hybrid Lagrangian–Eulerian model that unifies the strengths of both representations within a single neural simulator; (i) Experimental Validation: We conduct comprehensive experiments to validate that our model significantly enhances accuracy and stability over baseline methods. 2 Related Work Several prior works have explored Lagrangian fluid simulations. A prominent line of research focuses on purely graph-based methods to model particle interactions. GNS (Sanchez-Gonzalez et al., 2020) applies message-passing graph neural networks to learn the dynamics of fluids, granular media, and deformable bodies. Building on this foundation, several works have introduced physical constraints and solver mechanics to improve these networks. Toshev et al. (2024b) improves upon GNS by incorporating various components from standard SPH solvers. Furthermore, Toshev et al. (2023) introduces an architecture that enforces Euclidean symmetries, including translation, rotation, and reflection invariance, resulting in improved generalization, while Prantl et al. (2022) presents a novel method for guaranteeing linear momentum in the neural simulator. Subsequent works have sought to enhance efficiency and scalability by incorporating particle-grid transfers. Rochman-Sharabi et al. (2025) integrates the material point method into a neural framework, enabling learned particle–grid transfers and updates for efficient emulation of high-fidelity simulations. Xu et al. (2025) integrates a neural simulator with the classical material-point method (Johnson, 1996; Brackbill & Ruppel, 1986; Sulsky et al., 1994; Nairn, 2003) specifically to achieve real-time simulation. To standardize the evaluation of this growing body of literature, LagrangeBench (Toshev et al., 2024a) provides a comprehensive benchmark for particle-based fluid simulation learning, systematically comparing neural models across multiple physical configurations. Finally, Alkin et al. (2024) proposes a framework that generalizes across various discretizations. However, this approach falls under the domain of field learning as it requires ground truth particle coordinates to query velocity and cannot directly predict future particle positions. Due to this fundamental difference in problem formulation, it is not included in our comparative analysis. In the domain of Eulerian simulation, significant progress has been made using deep learning to accelerate solvers on fixed spatial discretizations. This includes (1) grid-based methods (Li et al., 2021, 2023b; Brandstetter et al., 2022; Williams et al., 2023; Taccari et al., 2022; Tran et al., 2023; Huang et al., 2024; Cao, 2021; Lu et al., 2021), and (2) mesh-based techniques designed to generalize across irregular geometries (Pfaff et al., 2021; Wu et al., 2024; Luo et al., 2025; Cao et al., 2023; Li et al., 2023a; Wu et al., 2023; Wang & Wang, 2024; Hao et al., 2023). While these approaches offer efficiency for problems with static topology, they do not address the specific challenges of Lagrangian tracking and dynamic neighborhood configurations central to the particle-based methods discussed in this work. Beyond forward simulation, neural models are widely utilized for foundational models (Hao et al., 2024; McCabe et al., 2024; Ye et al., 2024; Herde et al., 2024; Shen et al., 2024; Zhou et al., ), fluid field reconstruction (Li et al., 2025; Zhong et al., 2023; Mo & Magri, 2024; Yadav et al., 2025; Jing et al., 2024; He et al., 2022), ML-assisted classical solver (List et al., 2022; Sun et al., 2023; Greenfeld et al., 2019; Sappl et al., 2019), aerodynamic shape optimization (Elrefaie et al., 2025, 2024), and inverse design (Behrmann et al., 2019; Teng & Choromanska, 2019; Kruse et al., 2021). While these works fall within the broader area of data-driven neural PDE modeling, they address computational tasks distinct from the sequential forward time-stepping of Lagrangian particles and therefore lie outside the primary scope of our work. 3 Problem Statement We follow the setup in LagrangeBench (Toshev et al., 2024a) and FD-Bench (Wang et al., 2025) to focus on the weakly-compressible Navier–Stokes equation. dρdt=−ρ(∇⋅˙),d˙dt=−1ρ∇p+1Re∇2˙+1ρ, dρdt=-ρ (∇· x ), d xdt=- 1ρ∇ p+ 1Re∇^2 x+ 1ρ F, (1) where ρ is density, ˙ x is velocity, p is pressure, ReRe is the Reynolds number, and F is an external force field. This choice is made solely due to the availability of extensive benchmark datasets. Our design is not tied to this particular equation, and we expect the model to generalize to other fluid scenarios. We adopt a particle-based representation, where the solution is discretized into a set of particles carrying physical properties, denoted by ∈ℝNlag×d x ^N_lag× d, where NlagN_lag represents the number of Lagrangian particles. Here, d denotes the spatial dimension, e.g., 2D or 3D. Our objective is to learn a neural surrogate model Ψ that, given τ consecutive frames of particle positions, predicts the positions in the next frame: t+1=Ψ((t−τ+1):t) x_t+1= ( x_(t-τ+1):t ). Figure 1: Overall framework of the proposed model. We employ an Encoder – Downsampler – Processor – Upsampler – Decoder architecture. Lagrangian inputs are converted to Eulerian graphs, encoded via MLPs, and downsampled into sparse representations. The processor solves coarse PDE dynamics by particle aggregation and cross-attention, after which the upsampler and decoder predict final kinematics. 4 Methodology In this section, we present our Hybrid Lagrangian–Eulerian model, which adopts an Encoder – Downsampler – Processor – Upsampler – Decoder architecture. We describe the role of each component in the following subsections and Figure 1. Preprocess Given the initial particle positions (t−τ+1):t x_(t-τ+1):t, we compute the velocities (t−τ+2):t v_(t-τ+2):t using the backward finite difference method. We then construct the Lagrangian graph lag=(lag,Eclag)G_lag=(V_lag,Eclag), whose feature vector is defined as lagi=[ti‖(t−τ+2):ti‖typei∥Θ]V_lag^i= [ x_t^i\,\|\, v_(t-τ+2):t^i\,\|\,type^i\,\|\, ], where typeitype^i denotes the type of particle i, and Θ represents global variables (both are optional). The edge set ElagE_lag contains directed edge indices and edge attributes, where the attributes include the relative position vector and its norm. The edge indices are determined via a radius graph. We also define a corresponding Eulerian graph. Let ~∈ℝNeu×d x ^N_eu× d denote the positions of grid or mesh points in the spatial domain Ω , where NeuN_eu represents the number of Eulerian nodes. Unlike the Lagrangian particle positions, which evolve over time, the Eulerian nodes remain fixed. We choose NeuN_eu such that Neu≪NlagN_eu N_lag. To obtain Eulerian velocities, we follow Ramachandran et al. (2021) and aggregate the particle velocities using a quintic kernel W:ℝ≥0×ℝ>0→ℝ≥0W:R_≥ 0×R_>0 _≥ 0, defined by W(r,h)=σ2h2[(3−q)5−6(2−q)5+15(1−q)5]+, W(r,h)= _2h^2 [(3-q)^5-6(2-q)^5+15(1-q)^5 ]_+, (2) where q=rhq= rh, σ2=7478π _2= 7478π, and [⋅]+:=max(⋅,0)[\,·\,]_+:= (·,0). The Eulerian velocity aggregation is given by ~j=∑i=1NiW(‖~j−i‖,h)∑i=1NW(‖~j−i‖,h)+ε, v_j= _i=1^N v_i\,W(\| x_j- x_i\|,h) _i=1^NW(\| x_j- x_i\|,h)+ , (3) where ~j v_j denotes the velocity at Eulerian node j. The Eulerian interaction graph is defined as eu=(eu,Eeu)G_eu=(V_eu,E_eu), where eui=[~ti∥~(t−τ+2):ti]V_eu^i= [ x_t^i\,\|\, v_(t-τ+2):t^i ], and EeuE_eu contains directed edge indices that connect each node to its immediate horizontal and vertical neighbors, with edge features encoding the relative position vector and its norm. For regularly spaced Eulerian grids, omitting edge attributes and node positions does not negatively impact the model. Encoder: In the encoding stage, we employ four distinct MLPs denoted by ℰlag,E_lag,V, ℰlag,EE_lag,E, ℰeu,E_eu,V, and ℰeu,EE_eu,E to embed node and edge features into latent space. Formally, lag h_lag =ℰlag,(lag), =E_lag,V(V_lag), lag e_lag =ℰlag,E(Elag), =E_lag,E(E_lag), (4) eu h_eu =ℰeu,(eu), =E_eu,V(V_eu), eu e_eu =ℰeu,E(Eeu). =E_eu,E(E_eu). (5) Here, lag∈ℝNlag×dlatent h_lag ^N_lag× d_latent and eu∈ℝNeu×dlatent h_eu ^N_eu× d_latent. Downsampler: Our downsampler consists of M successive blocks, each comprising a SAGPooling (Lee et al., 2019) operation followed by message passing. Given a node feature ∈ℝN×dlatent h ^N× d_latent as well as edge index E with latent features e. SAGPooling first computes a scalar importance score for each node via a self-attention GNN: =GNN(,) s=GNN( h, e). Then, it selects the top-k nodes with the highest importance scores: TopK(,k)TopK(s,k), where k=⌈λN⌉k= λ N and λ is the downsmpling rate, and retain only the corresponding node features and the induced subgraph. Next, we apply a message passing GNN with the update rule i=∑(j→i)∈Eϕm(i,j,i,j)m_i= _(j→ i)∈ E _m(h_i,h_j,e_i,j) and i←ϕh(i,i)h_i← _h(h_i,m_i), where ϕm _m and ϕh _h are parametrized by MLPs. We apply this downsampling procedure independently to both the Lagrangian and Eulerian graphs, yielding latent node embeddings lag∈ℝN~lag×dlatentq_lag N_lag× d_latent and eu∈ℝN~eu×dlatentq_eu N_eu× d_latent, with associated particle or node coordinates y and ~ y, respectively. By construction, N~eu≪N~lag N_eu N_lag. To reduce the number of hyperparameters, we use a shared pooling rate λ for both the Lagrangian and Eulerian graphs. Empirically, this choice is sufficient to achieve strong performance. Nevertheless, if computational resources permit, we recommend assigning separate rates to the two graphs and performing hyperparameter tuning. The downsampler plays a critical role by reducing the redundancy inherent in particle-based fluid representations. In regions where the flow is smooth or nearly uniform, such as laminar zones, the dynamics of individual particles become indistinguishable. They move together or remain nearly static, offering little new information if treated separately. The downsampler condenses these redundant particles into a smaller set of representative super-nodes, capturing their shared properties without expending computation on unnecessary features. This reduction not only improves efficiency but also ensures that the model allocates its capacity to regions with complex, highly variable dynamics, where detailed interactions matter. The processor, introduced in the next section, then operates on this compressed representation, solving the governing equation at a global level. Without this downsampler, the processor would be forced to carry the full burden of redundant particle interactions, wasting computation while diluting attention away from the regions that matter most. Processor: We define a learnable aggregation kernel :ℝd×ℝd→ℝ>0K:R^d×R^d _>0, parametrized by MLPS, that maps a Lagrangian particle position and an Eulerian node position to a positive weight. For each particle i, K assigns a weight to each Eulerian node j, satisfying the normalization condition ∑j=1N~eu(i,~j)=1 _j=1 N_euK( y_i, y_j)=1. Let ∈ℝN~lag×N~eu W N_lag× N_eu denote the weight matrix produced by this kernel. Using W, we aggregate the downsampled Lagrangian particles onto the downsampled Eulerian nodes: lag=⊤lag∈ℝN~eu×dlatent. _lag= W q_lag N_eu× d_latent. We then apply multi-head self-attention to the aggregated Lagrangian feature, =MHA(Query=lag,Key=lag,Value=lag), A=MHA\! (Query=Z_lag,\ Key=Z_lag,\ Value=Z_lag ), (6) producing coarse predictions based solely on the aggregated Lagrangian representation. Next, we perform cross-attention between the aggregated Lagrangian feature and the downsampled Eulerian feature, =MHA(Query=lag,Key=eu,Value=eu), C=MHA\! (Query= Z_lag,\ Key= q_eu,\ Value= q_eu ), (7) allowing the model to jointly consider Lagrangian and Eulerian features. The outputs are then combined and scattered back to the downsampled Lagrangian nodes: =+α,^lag=∈ℝN~lag×dlatent, F= A+α\, C, q_lag= W\, F N_lag× d_latent, (8) where α is a learnable parameter. The processor is inspired by PIC (Francis Harvey Harlow, 1955; Dawson, 1983; Tskhakaya, 2008), adapted to the hierarchical latent space. In classical PIC, particles are projected onto a grid using voxelization to solve field equations. The learnable aggregation kernel K plays the role of this projection operator, while the self-attention on aggregated Lagrangian feature functions as a latent field solver on a coarse level. The Lagrangian formulation naturally tracks the material derivative along particle trajectories, providing an accurate representation of advection-dominated dynamics, while the Eulerian formulation encodes the spatial structure of the flow field on a fixed stencil, facilitating the computation of differential operators such as divergence and vorticity that govern incompressibility. Self-attention on aggregated Lagrangian representation alone can only mix Lagrangian-derived information, whereas cross-attention is the channel by which each Lagrangian query token retrieves the field-derived quantities it cannot reconstruct internally. This complementarity is also what makes the design stable over long autoregressive rollouts. Pure-Lagrangian neural surrogates are known to drift, since small per-step errors compound into unphysical particle distributions over hundreds of steps. Because Eulerian nodes never move and their features are tied to fixed spatial locations, the Eulerian branch acts as a spatial anchor. Cross-attention re-grounds the Lagrangian update at every step against this stable, field-based reference, suppressing error accumulation and improving long-horizon fidelity. We repeat this layer for L steps. Note that the downsampled Lagrangian particles and Eulerian nodes are fairly small. The calculation of the weight matrix and aggregation does not incur significant computational overhead. Upsampler: With ^lag q_lag, we reconstruct the full particle latent representation ^lag∈ℝNlag×dlatent h_lag ^N_lag× d_latent by executing M successive upsampling operations, each coupled with message-passing steps, while integrating a residual link from the downsampling stage. Decoder: The latent state ^lag h_lag is passed through a learnable MLP decoder to produce the predicted position, velocity, or acceleration. If the network output is velocity or acceleration, position is computed with Euler integration. Remark: While NeuralMPM (Rochman-Sharabi et al., 2025) shares a conceptual link with our work in that it also aggregates particles to Eulerian nodes, our work is not simply an extension of this method with additional components. The voxelization approach used in this prior work is highly restrictive, limiting it exclusively to regular grids, and they only perform a single aggregation step at the beginning of the process before scattering back at the end. In contrast, our method is highly adaptive. We utilize a completely different aggregation strategy that performs distinct particle aggregations at each layer, which we believe significantly enhances the model’s expressivity. Furthermore, because our Eulerian representation is defined on a general mesh graph rather than a fixed grid, it naturally extends to arbitrary domains and adapts to complex geometries and deformable domains, a capability we demonstrate in subsequent experiments. By combining the hierarchical, layer-wise aggregation with a cross-attention mechanism, our framework achieves geometric generality, efficiency, and fidelity in a manner entirely separate from existing methods. Metric GNS SEGNN GraphUnet GraphTransformer AdvDIFFormer NeuralMPM PhysicsNFP Ours 2D TGV MSEMSE 7.920×10−37.920× 10^-3 1.260×10−21.260× 10^-2 5.356×10−35.356× 10^-3 2.963×10−22.963× 10^-2 1.586×10−21.586× 10^-2 3.051×10−23.051× 10^-2 7.140×10−37.140× 10^-3 4.551×10−34.551× 10^-3↑15.02% KEKE 1.795×10−31.795× 10^-3 2.854×10−32.854× 10^-3 1.246×10−31.246× 10^-3 4.729×10−34.729× 10^-3 2.339×10−32.339× 10^-3 6.402×10−36.402× 10^-3 1.652×10−31.652× 10^-3 1.082×10−31.082× 10^-3↑13.20% SinkhornSinkhorn 4.150×10−54.150× 10^-5 3.722×10−43.722× 10^-4 3.767×10−53.767× 10^-5 3.099×10−33.099× 10^-3 1.744×10−41.744× 10^-4 9.607×10−49.607× 10^-4 3.323×10−53.323× 10^-5 1.422×10−51.422× 10^-5↑57.22% 2D LDC MSEMSE 3.542×10−53.542× 10^-5 6.293×10−66.293× 10^-6 4.413×10−64.413× 10^-6 6.490×10−46.490× 10^-4 3.955×10−53.955× 10^-5 2.308×10−32.308× 10^-3 2.629×10−62.629× 10^-6 1.868×10−61.868× 10^-6↑28.94% KEKE 7.880×10−67.880× 10^-6 1.043×10−51.043× 10^-5 1.163×10−71.163× 10^-7 1.635×10−41.635× 10^-4 1.316×10−51.316× 10^-5 5.724×10−25.724× 10^-2 1.068×10−71.068× 10^-7 7.240×10−87.240× 10^-8↑32.23% SinkhornSinkhorn 8.239×10−68.239× 10^-6 1.388×10−51.388× 10^-5 2.227×10−62.227× 10^-6 1.803×10−51.803× 10^-5 2.214×10−52.214× 10^-5 3.172×10−43.172× 10^-4 8.768×10−78.768× 10^-7 4.675×10−74.675× 10^-7↑46.68% 2D RPF MSEMSE 2.629×10−12.629× 10^-1 5.779×10−25.779× 10^-2 4.067×10−24.067× 10^-2 8.835×10−18.835× 10^-1 2.842×10−12.842× 10^-1 1.044×10−11.044× 10^-1 3.862×10−23.862× 10^-2 2.426×10−22.426× 10^-2↑37.19% KEKE 7.961×10−27.961× 10^-2 6.374×10−26.374× 10^-2 2.471×10−22.471× 10^-2 4.104×1004.104× 10^0 5.848×10−25.848× 10^-2 2.694×10−22.694× 10^-2 3.181×10−23.181× 10^-2 1.595×10−21.595× 10^-2↑35.42% SinkhornSinkhorn 3.627×10−23.627× 10^-2 2.134×10−22.134× 10^-2 7.909×10−37.909× 10^-3 4.120×10−14.120× 10^-1 4.812×10−24.812× 10^-2 2.163×10−22.163× 10^-2 1.684×10−21.684× 10^-2 6.140×10−36.140× 10^-3↑22.37% 2D DAM MSEMSE 2.363×10−42.363× 10^-4 9.245×10−49.245× 10^-4 1.505×10−41.505× 10^-4 4.552×10−34.552× 10^-3 6.360×10−36.360× 10^-3 2.010×10−42.010× 10^-4 1.458×10−41.458× 10^-4 9.171×10−59.171× 10^-5↑37.10% KEKE 9.461×10−79.461× 10^-7 1.817×10−61.817× 10^-6 9.235×10−79.235× 10^-7 3.783×10−43.783× 10^-4 6.238×10−46.238× 10^-4 6.080×10−66.080× 10^-6 7.856×10−77.856× 10^-7 4.709×10−74.709× 10^-7↑40.06% SinkhornSinkhorn 2.131×10−42.131× 10^-4 1.935×10−41.935× 10^-4 3.710×10−43.710× 10^-4 7.751×10−37.751× 10^-3 8.843×10−38.843× 10^-3 2.196×10−42.196× 10^-4 1.886×10−41.886× 10^-4 1.216×10−41.216× 10^-4↑35.55% Table 1: Comparison of MSE, Kinetic Energy (KE) error, and Sinkhorn divergence across five rollout steps on 2D datasets. The relative performance gain against the strongest baseline is computed as Baseline−OursBaseline×100% Baseline-OursBaseline× 100\%. 5 Experiment Dataset We use Lagrangian particle-based datasets (Toshev et al., 2024a) solved based on the weakly-compressible Navier–Stokes equations using SPH (Desbrun & Cani, 1996; Hoover, 2006; Monaghan, 1992) with the quintic kernel. The Taylor–Green Vortex (TGV) contains a periodic-domain simulation initialized from an analytical solution without any external forcing. The Lid-Driven Cavity Flow (LDC) introduces no-slip boundaries, including a moving top lid that drives the internal circulation, creating sharp shear layers and corner vortices. Reverse Poiseuille Flow (RPF) features opposing constant body forces in the upper and lower halves of the domain. Dam Break (DAM) captures the transient free-surface flow following the sudden release of a water column into an open downstream region, resulting in splashing and wave breaking. The datasets are sampled every 100 steps of the ground truth solver. Thus, a five-step rollout corresponds to 500 simulation steps. We refer the readers to Appendix C for additional details on the datasets. Task Setup and Baselines To evaluate our model, we benchmark it against a range of strong baselines. We include GNS (Sanchez-Gonzalez et al., 2020) and SEGNN (Brandstetter et al., 2021), the top-performing models reported in LagrangeBench (Toshev et al., 2024a). Graph-Unet (Gao & Ji, 2019) and GraphTransformer (Dwivedi & Bresson, 2021) are also considered due to their wide adoption and proven effectiveness in graph-based learning tasks. In addition, we evaluate against AdvDIFFormer (Wu et al., 2025) and PhysicsNFP (Jiang et al., 2025), which are recently proposed architectures that have achieved competitive results in diverse physical simulation settings. We also include NeuralMPM (Rochman-Sharabi et al., 2025), as it represents the work most closely related to our approach. We refer to Appendix D for detailed training and implementation of all models. Metric GraphUnet PhysicsNFP Ours 3D TGV MSEMSE 4.512×10−14.512× 10^-1 3.034×1003.034× 10^0 3.649×10−13.649× 10^-1↑19.14% KEKE 3.719×1003.719× 10^0 1.176×1021.176× 10^2 3.215×1003.215× 10^0↑13.56% SinkhornSinkhorn 8.799×10−28.799× 10^-2 7.560×10−17.560× 10^-1 1.094×10−21.094× 10^-2↑87.56% 3D LDC MSEMSE 4.810×10−34.810× 10^-3 2.455×10−22.455× 10^-2 3.059×10−33.059× 10^-3↑36.40% KEKE 1.464×10−41.464× 10^-4 1.807×10−21.807× 10^-2 1.271×10−41.271× 10^-4↑13.17% SinkhornSinkhorn 4.647×10−44.647× 10^-4 3.452×10−23.452× 10^-2 2.436×10−42.436× 10^-4↑47.59% 3D RPF MSEMSE 2.161×10−22.161× 10^-2 6.386×10−26.386× 10^-2 1.810×10−21.810× 10^-2↑16.25% KEKE 1.167×10−21.167× 10^-2 7.377×10−27.377× 10^-2 1.022×10−21.022× 10^-2↑12.38% SinkhornSinkhorn 3.721×10−43.721× 10^-4 8.551×10−38.551× 10^-3 3.018×10−43.018× 10^-4↑18.89% Table 2: Comparison of MSE, Kinetic Energy (KE) error, and sinkhorn divergence across five rollout steps on 3D datasets. The relative performance gain against the strongest baseline is computed as Baseline−OursBaseline×100% Baseline-OursBaseline× 100\%. Model performance is evaluated over five-step rollouts using three metrics. The mean squared error (MSE) measures the accuracy of trajectory predictions. The kinetic energy error (KE) captures global discrepancies in the fluid’s energy. Finally, the Sinkhorn divergence quantifies the distance between predicted and reference particle distributions using optimal transport. All models are configured to directly predict particle positions from unperturbed inputs. We deliberately exclude auxiliary training strategies found in prior literature (Pfaff et al., 2021; Sanchez-Gonzalez et al., 2020) to isolate the specific contributions of our proposed modules. Metric GraphUnet PhysicsNFP Ours 3600 MSEMSE 5.356×10−35.356× 10^-3 7.140×10−37.140× 10^-3 4.551×10−34.551× 10^-3↑15.02% KEKE 1.246×10−31.246× 10^-3 1.652×10−31.652× 10^-3 1.082×10−31.082× 10^-3↑13.20% SinkhornSinkhorn 3.767×10−53.767× 10^-5 3.323×10−53.323× 10^-5 1.422×10−51.422× 10^-5↑57.22% 4489 MSEMSE 6.396×10−36.396× 10^-3 1.623×10−21.623× 10^-2 4.572×10−34.572× 10^-3↑28.52% KEKE 1.530×10−31.530× 10^-3 4.895×10−34.895× 10^-3 1.087×10−31.087× 10^-3↑28.97% SinkhornSinkhorn 2.229×10−52.229× 10^-5 1.397×10−31.397× 10^-3 1.171×10−51.171× 10^-5↑47.47% 10000 MSEMSE 6.798×10−36.798× 10^-3 1.796×10−21.796× 10^-2 5.999×10−35.999× 10^-3↑11.75% KEKE 1.697×10−31.697× 10^-3 6.252×10−36.252× 10^-3 1.468×10−31.468× 10^-3↑13.47% SinkhornSinkhorn 2.444×10−52.444× 10^-5 1.540×10−31.540× 10^-3 9.640×10−69.640× 10^-6↑60.55% 40000 MSEMSE 6.923×10−36.923× 10^-3 1.669×10−21.669× 10^-2 5.063×10−35.063× 10^-3↑26.87% KEKE 1.403×10−31.403× 10^-3 6.631×10−36.631× 10^-3 1.205×10−31.205× 10^-3↑14.10% SinkhornSinkhorn 5.086×10−65.086× 10^-6 2.183×10−32.183× 10^-3 4.295×10−64.295× 10^-6↑15.55% Table 3: Comparison of MSE, Kinetic Energy (KE) error, and Sinkhorn divergence on higher resolution simulations. The relative performance gain against the strongest baseline is computed as Baseline−OursBaseline×100% Baseline-OursBaseline× 100\%. 5.1 Main Results Table 1 summarizes the quantitative results for the 2D experiments, where our proposed model demonstrates consistently superior predictive accuracy. In particular, our method significantly outperforms NeuralMPM, validating the efficacy of our architectural design choices. The most substantial gains are observed in the Sinkhorn divergence metric, with improvements of 57.22%57.22\% and 46.68%46.68\% on the TGV and LDC datasets, respectively. This indicates that our model is far more effective at capturing the underlying distributional geometry and global structure of the data than the baseline methods, rather than merely fitting local features. This success is driven by the downsampler’s ability to prioritize high-importance areas, while the Eulerian cross-attention ensures the predictions remain grounded in a consistent global structure. 3D Datasets To further assess the effectiveness of our approach, we compare the proposed model against the two strongest baselines, GraphUnet and PhysicsNFP, on challenging 3D datasets. The results, summarized in Table 2, demonstrate that our model consistently surpasses both baselines across evaluation metrics. This indicates that the advantages of our design generalize beyond 2D settings and remain robust even in the more complex 3D scenarios. 5.2 Higher Resolution and Scalability Simulation High-resolution simulations are essential for providing a detailed and accurate description of fluid flow, enabling more precise queries of velocity and pressure at arbitrary locations. However, they introduce two compounding challenges: the inherent difficulty of modeling complex micro-scale details and the computational burden of processing massive particle counts. Our approach is specifically designed to address these distinct hurdles. To manage the large particle counts, our downsampler condenses indistinguishable particles, effectively reducing computational overhead. Simultaneously, our hierarchical framework allows for the separate modeling of micro-scale details and coarse-level large structures. We validate this capability using the Taylor Green Vortex, varying particle spacing from 0.02 to 0.005, resulting in 3,600 to 40,000 particles. This experiment serves as a scalability test for both the model’s scalability and its ability to resolve intricate flow features. We refer the readers to Appendix C.1 for dataset details. The results are summarized in Table 3, where we compare our proposed model against the two strongest baselines, GraphUnet and PhysicsNFP, introduced in Section 5.1. Our model consistently outperforms both baselines. In particular, PhysicsNFP shows a noticeable drop in performance as the resolution increases, indicating that hierarchical architectures such as GraphUnet and our model are better suited for simulations with larger particle counts. Metric GraphUnet PhysicsNFP Ours MSEMSE 5.572×10−55.572× 10^-5 1.150×10−21.150× 10^-2 4.955×10−54.955× 10^-5↑11.07% KEKE 1.207×10−71.207× 10^-7 9.560×10−39.560× 10^-3 9.980×10−89.980× 10^-8↑17.31% SinkhornSinkhorn 1.949×10−61.949× 10^-6 2.130×10−22.130× 10^-2 1.682×10−61.682× 10^-6↑13.70% Table 4: Comparison of MSE, Kinetic Energy (KE) error, and Sinkhorn divergence on fluid-solid interaction over an irregular domain. The relative performance gain against the strongest baseline is computed as Baseline−OursBaseline×100% Baseline-OursBaseline× 100\%. 5.3 Fluid-Solid Interaction Simulation The voxelization strategy adopted in Rochman-Sharabi et al. (2025) is restricted to fixed, rectangular grids, which prohibits generalization to irregular domains that frequently arise in fluid–solid interaction problems. We visualize the data in Figure 2. To evaluate our model’s ability to handle such scenarios, we conduct experiments on the Flow Around Cylinder task within an irregular domain, comparing against the two strongest baselines, GraphUnet and PhysicsNFP. We refer the readers to Appendix C.2 for dataset details. The results, reported in Table 4, demonstrate that our proposed model consistently outperforms both baselines. These findings highlight not only the effectiveness of our design but also its broader applicability to complex geometries where voxelization-based approaches face fundamental limitations. Metric GraphUnet PhysicsNFP Ours MSEMSE 3.602×10−53.602× 10^-5 1.379×10−41.379× 10^-4 2.421×10−52.421× 10^-5↑32.79% KEKE 1.469×10−51.469× 10^-5 7.203×10−67.203× 10^-6 4.229×10−64.229× 10^-6↑41.29% SinkhornSinkhorn 5.742×10−55.742× 10^-5 1.646×10−41.646× 10^-4 3.107×10−53.107× 10^-5↑45.89% Table 5: Comparison of MSE, Kinetic Energy (KE) error, and Sinkhorn divergence on deforming domain. The relative performance gain against the strongest baseline is computed as Baseline−OursBaseline×100% Baseline-OursBaseline× 100\%. 5.4 Deformable Domain Simulation Lagrangian simulations excel at modeling deformations and free-surface boundaries. In addition to the dam break scenario evaluated in Section 5.1, we also test our model on a dynamically deforming domain. Specifically, we simulate a 2D circular drop of fluid that is subjected to a velocity field, stretching it into an ellipse (Monaghan, 1994). Because the fluid is incompressible, its total area remains strictly conserved during this deformation. We provide a visualization of this process in Figure 2. Figure 2: Visualization of datasets with irregular and deforming domains. Top row: Periodic flow around a fixed cylinder centered in the domain, creating an irregular flow field. Bottom row: A fluid droplet subjected to a velocity field, demonstrating continuous stretching and domain deformation. The experimental results are presented in Table 5. As in previous experiments, we compare our approach against GraphUnet and PhysicsNFP. The results show that our proposed model achieves superior performance. This is expected, as our architectural design is inherently well-suited for continuously deforming domains. 5.5 Extended Rollout We investigate our model’s performance on extended rollouts using the Water dataset (Sanchez-Gonzalez et al., 2020). This dataset features a volume of water dropped into a box, involving complex splashing dynamics. We compare our model against GraphUnet, as PhysicsNFP is unable to handle the varying particle counts inherent in this dataset. We apply all auxiliary training strategies to both the baseline and our proposed models. These models predict higher-order temporal derivatives and utilize Euler integration to calculate subsequent positions. Additionally, we incorporate noise into the training data to enhance rollout stability. As shown in Table 6, our model maintains superior performance over extended horizons, surpassing the strongest baseline at every step. Notably, the kinetic energy error at the 200-step rollout is orders of magnitude smaller than that of the baseline. We attribute this to the cross-attention with Eulerian features, which helps the model correct particle drift. Steps Model Metric 25 50 75 100 150 200 MSEMSE 1.801×10−41.801× 10^-4 1.080×10−31.080× 10^-3 2.810×10−32.810× 10^-3 6.234×10−36.234× 10^-3 2.111×10−22.111× 10^-2 4.201×10−24.201× 10^-2 KEKE 3.252×10−23.252× 10^-2 8.987×10−18.987× 10^-1 6.726×1006.726× 10^0 3.114×1013.114× 10^1 3.652×1023.652× 10^2 1.545×1031.545× 10^3 GraphUnet SinkhornSinkhorn 1.225×10−41.225× 10^-4 4.905×10−44.905× 10^-4 1.398×10−31.398× 10^-3 3.966×10−33.966× 10^-3 1.967×10−21.967× 10^-2 3.235×10−23.235× 10^-2 MSE 1.612×10−51.612× 10^-5 8.718×10−58.718× 10^-5 5.244×10−55.244× 10^-5 1.895×10−31.895× 10^-3 9.250×10−39.250× 10^-3 1.681×10−21.681× 10^-2 KE 1.015×10−51.015× 10^-5 1.237×10−41.237× 10^-4 3.230×10−33.230× 10^-3 4.894×10−24.894× 10^-2 1.212×1001.212× 10^0 3.735×1003.735× 10^0 Ours Sinkhorn 1.397×10−51.397× 10^-5 1.314×10−41.314× 10^-4 9.414×10−49.414× 10^-4 3.680×10−33.680× 10^-3 1.716×10−21.716× 10^-2 3.042×10−23.042× 10^-2 Table 6: Comparison of MSE, Kinetic Energy (KE) error, and Sinkhorn divergence on the Water dataset over 200 rollout steps. We compare against the strongest baseline, GraphUNet. PhysicsNFP is excluded, as it cannot handle a varying number of particles. 6 Additional Experiment Ablation Study and Model Analysis: We conduct extensive ablation studies to evaluate each component of our model, alongside an analysis of key parameters, in Appendix E. We demonstrate that every proposed component contributes positively to model performance across all metrics. Additional Baselines on Eulerian Mesh Simulator: Because Eulerian mesh-based neural simulators operate on graph structures analogous to our Lagrangian framework, they are inherently compatible with our setting. Consequently, we compare our approach against a comprehensive suite of state-of-the-art Eulerian mesh models in Appendix F. Additional Experiments on Diverse Materials: While the primary focus of this work is Lagrangian fluid simulation, other materials such as sand and deformable goop can also adopt a Lagrangian representation. In Appendix G, we demonstrate that our model generalizes well to these diverse material domains and achieves strong performance. Computational Complexity: We report the runtime and memory consumption of the evaluated models in Appendix H. Despite processing both Lagrangian and Eulerian representations, our approach leverages parallelization to maintain highly competitive inference speeds. 7 Conclusion In this work, we investigate the limitations of Lagrangian neural simulators and identify the sources of their performance gap relative to Eulerian approaches. To overcome these challenges, we introduced a Hybrid Lagrangian–Eulerian model that unifies the complementary strengths of both representations through particle downsampling, particle aggregation, and cross-attention. Our experiments demonstrate that this hybrid framework achieves higher accuracy, stronger robustness, and improved stability across diverse fluid regimes. These results highlight the promise of bridging particle- and grid-based formulations to advance the next generation of neural fluid simulators. References Agertz et al. (2007) Agertz, O., Moore, B., Stadel, J., Potter, D., Miniati, F., Read, J., Mayer, L., Gawryszczak, A., Kravtsov, A., Nordlund, A., Pearce, F., Quilis, V., Rudd, D., Springel, V., Stone, J., Tasker, E., Teyssier, R., Wadsley, J., and Walder, R. Fundamental differences between sph and grid methods: Simulating fluids using sph and grid techniques. Monthly Notices of the Royal Astronomical Society, 380(3):963–978, August 2007. ISSN 0035-8711. doi: 10.1111/j.1365-2966.2007.12183.x. URL http://dx.doi.org/10.1111/j.1365-2966.2007.12183.x. Alkin et al. (2024) Alkin, B., Fürst, A., Schmid, S., Gruber, L., Holzleitner, M., and Brandstetter, J. Universal physics transformers. arXiv preprint arXiv:2402.12365, 2024. Behrmann et al. (2019) Behrmann, J., Grathwohl, W., Chen, R. T., Duvenaud, D., and Jacobsen, J.-H. Invertible residual networks. In International Conference on Machine Learning, p. 573–582. PMLR, 2019. Brackbill & Ruppel (1986) Brackbill, J. U. and Ruppel, H. M. FLIP: A method for adaptively zoned, particle-in-cell calculations of fluid flows in two dimensions. Journal of Computational Physics, 65(2):314–343, 1986. doi: 10.1016/0021-9991(86)90211-1. Brandstetter et al. (2021) Brandstetter, J., Hesselink, R., van der Pol, E., Bekkers, E., and Welling, M. Geometric and physical quantities improve e(3) equivariant message passing. 2021. Brandstetter et al. (2022) Brandstetter, J., Worrall, D. E., and Welling, M. Message passing neural PDE solvers. In International Conference on Learning Representations, 2022. URL https://openreview.net/forum?id=vSix3HPYKSU. Cao (2021) Cao, S. Choose a transformer: Fourier or galerkin, 2021. URL https://arxiv.org/abs/2105.14995. Cao et al. (2023) Cao, Y., Chai, M., Li, M., and Jiang, C. Efficient learning of mesh-based physical simulation with bi-stride multi-scale graph neural network. In International Conference on Machine Learning, 2023. URL https://openreview.net/forum?id=2Mbo7IEtZW. Dawson (1983) Dawson, J. M. Particle simulation of plasmas. Rev. Mod. Phys., 55:403–447, Apr 1983. doi: 10.1103/RevModPhys.55.403. URL https://link.aps.org/doi/10.1103/RevModPhys.55.403. Dehnen & Aly (2012) Dehnen, W. and Aly, H. Improving convergence in smoothed particle hydrodynamics simulations without pairing instability: Sph without pairing instability. Monthly Notices of the Royal Astronomical Society, 425(2):1068–1082, August 2012. ISSN 0035-8711. doi: 10.1111/j.1365-2966.2012.21439.x. URL http://dx.doi.org/10.1111/j.1365-2966.2012.21439.x. Desbrun & Cani (1996) Desbrun, M. and Cani, M.-P. Smoothed particles: a new paradigm for animating highly deformable bodies. In Proceedings of Eurographics Workshop on Computer Animation and Simulation, Poitiers, France, August 1996. Dwivedi & Bresson (2021) Dwivedi, V. P. and Bresson, X. A generalization of transformer networks to graphs, 2021. URL https://arxiv.org/abs/2012.09699. Elrefaie et al. (2024) Elrefaie, M., Morar, F., Dai, A., and Ahmed, F. Drivaernet++: A large-scale multimodal car dataset with computational fluid dynamics simulations and deep learning benchmarks. In Globerson, A., Mackey, L., Belgrave, D., Fan, A., Paquet, U., Tomczak, J., and Zhang, C. (eds.), Advances in Neural Information Processing Systems, volume 37, p. 499–536. Curran Associates, Inc., 2024. Elrefaie et al. (2025) Elrefaie, M., Dai, A., and Ahmed, F. Drivaernet: A parametric car dataset for data-driven aerodynamic design and prediction. Journal of Mechanical Design, 147(4), March 2025. ISSN 1528-9001. doi: 10.1115/1.4068104. URL http://dx.doi.org/10.1115/1.4068104. Francis Harvey Harlow (1955) Francis Harvey Harlow, Martha Evans, R. D. R. A machine calculation method for hydrodynamic problems. Report, Los Alamos Scientific Laboratory of the University of California, 1955. Gao & Ji (2019) Gao, H. and Ji, S. Graph u-nets. In International Conference on Machine Learning, p. 2083–2092, 2019. Gingold & Monaghan (1977) Gingold, R. A. and Monaghan, J. J. Smoothed particle hydrodynamics: Theory and application to non-spherical stars. Mon. Not. Roy. Astron. Soc., 181:375, 1977. Greenfeld et al. (2019) Greenfeld, D., Galun, M., Basri, R., Yavneh, I., and Kimmel, R. Learning to optimize multigrid PDE solvers. In International Conference on Machine Learning, p. 2415–2423. PMLR, 2019. Hao et al. (2023) Hao, Z., Ying, C., Wang, Z., Su, H., Dong, Y., Liu, S., Cheng, Z., Zhu, J., and Song, J. Gnot: A general neural operator transformer for operator learning. arXiv preprint arXiv:2302.14376, 2023. Hao et al. (2024) Hao, Z., Su, C., Liu, S., Berner, J., Ying, C., Su, H., Anandkumar, A., Song, J., and Zhu, J. Dpot: Auto-regressive denoising operator transformer for large-scale pde pre-training. arXiv preprint arXiv:2403.03542, 2024. He et al. (2022) He, X., Wang, Y., and Li, J. Flow completion network: Inferring the fluid dynamics from incomplete flow information using graph neural networks. Physics of Fluids, 34(8), August 2022. ISSN 1089-7666. doi: 10.1063/5.0097688. URL http://dx.doi.org/10.1063/5.0097688. Herde et al. (2024) Herde, M., Raonic, B., Rohner, T., Käppeli, R., Molinaro, R., de Bezenac, E., and Mishra, S. Poseidon: Efficient foundation models for PDEs. In The Thirty-eighth Annual Conference on Neural Information Processing Systems, 2024. URL https://openreview.net/forum?id=JC1VKK3UXk. Hoover (2006) Hoover, W. G. Smooth Particle Applied Mechanics: The State of the Art. World Scientific, Singapore, 2006. Huang et al. (2024) Huang, J., Yang, G., Wang, Z., and Park, J. J. Diffusionpde: Generative pde-solving under partial observation, 2024. URL https://arxiv.org/abs/2406.17763. Jiang et al. (2025) Jiang, H., Wang, J., Zhu, X., and He, Y. Topology-aware neural flux prediction guided by physics. In Forty-second International Conference on Machine Learning, 2025. URL https://openreview.net/forum?id=8Um2YotdbD. Jing et al. (2024) Jing, G., Wang, H., Li, X., Wang, G., and Yang, Y. An airflow velocity field reconstruction method with sparse or incomplete data using physics-informed neural network. Journal of Building Engineering, 88:109231, July 2024. Published: 1 July 2024. Johnson (1996) Johnson, N. L. The legacy and future of CFD at Los Alamos. In Proceedings of the 1996 Canadian CFD Conference. Office of Scientific and Technical Information (OSTI), 1996. OSTI 244662. Kingma & Ba (2017) Kingma, D. P. and Ba, J. Adam: A method for stochastic optimization, 2017. URL https://arxiv.org/abs/1412.6980. Kruse et al. (2021) Kruse, J., Ardizzone, L., Rother, C., and Köthe, U. Benchmarking invertible architectures on inverse problems. arXiv preprint, 2021. Lee et al. (2019) Lee, J., Lee, I., and Kang, J. Self-attention graph pooling. In Proceedings of the 36th International Conference on Machine Learning, 09–15 Jun 2019. Li et al. (2025) Li, R., Wan, G., Huang, Z., Liu, Z., Wang, H., Luo, X., Wang, W., and Sun, Y. Flow field reconstruction with sensor placement policy learning. In The Thirty-ninth Annual Conference on Neural Information Processing Systems, 2025. URL https://openreview.net/forum?id=1esFEjMUBS. Li et al. (2020) Li, Z., Kovachki, N., Azizzadenesheli, K., Liu, B., Bhattacharya, K., Stuart, A., and Anandkumar, A. Neural operator: Graph kernel network for partial differential equations, 2020. URL https://arxiv.org/abs/2003.03485. Li et al. (2021) Li, Z., Kovachki, N. B., Azizzadenesheli, K., liu, B., Bhattacharya, K., Stuart, A., and Anandkumar, A. Fourier neural operator for parametric partial differential equations. In International Conference on Learning Representations, 2021. URL https://openreview.net/forum?id=c8P9NQVtmnO. Li et al. (2023a) Li, Z., Kovachki, N. B., Choy, C., Li, B., Kossaifi, J., Otta, S. P., Nabian, M. A., Stadler, M., Hundt, C., Azizzadenesheli, K., and Anandkumar, A. Geometry-informed neural operator for large-scale 3d PDEs. In Thirty-seventh Conference on Neural Information Processing Systems, 2023a. URL https://openreview.net/forum?id=86dXbqT5Ua. Li et al. (2023b) Li, Z., Zheng, H., Kovachki, N., Jin, D., Chen, H., Liu, B., Azizzadenesheli, K., and Anandkumar, A. Physics-informed neural operator for learning partial differential equations, 2023b. URL https://arxiv.org/abs/2111.03794. List et al. (2022) List, B., Chen, L.-W., and Thuerey, N. Learned turbulence modelling with differentiable fluid solvers: physics-based loss functions and optimisation horizons. Journal of Fluid Mechanics, 949:A25, 2022. Lu et al. (2021) Lu, L., Jin, P., Pang, G., Zhang, Z., and Karniadakis, G. E. Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators. Nature Machine Intelligence, 3(3):218–229, 2021. Luo et al. (2025) Luo, H., Wu, H., Zhou, H., Xing, L., Di, Y., Wang, J., and Long, M. Transolver++: An accurate neural solver for PDEs on million-scale geometries. In Forty-second International Conference on Machine Learning, 2025. URL https://openreview.net/forum?id=AM7iAh0krx. McCabe et al. (2024) McCabe, M., Blancard, B. R.-S., Parker, L. H., Ohana, R., Cranmer, M., Bietti, A., Eickenberg, M., Golkar, S., Krawezik, G., Lanusse, F., Pettee, M., Tesileanu, T., Cho, K., and Ho, S. Multiple physics pretraining for spatiotemporal surrogate models. In The Thirty-eighth Annual Conference on Neural Information Processing Systems, 2024. URL https://openreview.net/forum?id=DKSI3bULiZ. Mo & Magri (2024) Mo, Y. and Magri, L. Reconstructing unsteady flows from sparse, noisy measurements with a physics-constrained convolutional neural network, 2024. URL https://arxiv.org/abs/2409.00260. Monaghan (1994) Monaghan, J. Simulating free surface flows with sph. Journal of Computational Physics, 110(2):399–406, 1994. ISSN 0021-9991. doi: https://doi.org/10.1006/jcph.1994.1034. URL https://w.sciencedirect.com/science/article/pii/S0021999184710345. Monaghan (1992) Monaghan, J. J. Smoothed particle hydrodynamics. Annual Review of Astronomy and Astrophysics, 30:543–574, 1992. doi: 10.1146/annurev.a.30.090192.002551. Nairn (2003) Nairn, J. A. Material point method calculations with explicit cracks. Computer Modeling in Engineering & Sciences, 4(6):649–664, 2003. doi: 10.3970/cmes.2003.004.649. Paszke et al. (2019) Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Köpf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., Bai, J., and Chintala, S. Pytorch: An imperative style, high-performance deep learning library, 2019. URL https://arxiv.org/abs/1912.01703. Pfaff et al. (2021) Pfaff, T., Fortunato, M., Sanchez-Gonzalez, A., and Battaglia, P. W. Learning mesh-based simulation with graph networks. In International Conference on Learning Representations, 2021. Prantl et al. (2022) Prantl, L., Ummenhofer, B., Koltun, V., and Thuerey, N. Guaranteed conservation of momentum for learning particle-based fluid dynamics. In Oh, A. H., Agarwal, A., Belgrave, D., and Cho, K. (eds.), Advances in Neural Information Processing Systems, 2022. URL https://openreview.net/forum?id=6niwHlzh10U. Price (2012) Price, D. J. Smoothed particle hydrodynamics and magnetohydrodynamics. Journal of Computational Physics, 231(3):759–794, February 2012. ISSN 0021-9991. doi: 10.1016/j.jcp.2010.12.011. URL http://dx.doi.org/10.1016/j.jcp.2010.12.011. Ramachandran et al. (2021) Ramachandran, P., Bhosale, A., Puri, K., Negi, P., Muta, A., Dinesh, A., Menon, D., Govind, R., Sanka, S., Sebastian, A. S., Sen, A., Kaushik, R., Kumar, A., Kurapati, V., Patil, M., Tavker, D., Pandey, P., Kaushik, C., Dutt, A., and Agarwal, A. PySPH: A Python-based Framework for Smoothed Particle Hydrodynamics. ACM Transactions on Mathematical Software, 47(4):1–38, December 2021. ISSN 0098-3500, 1557-7295. doi: 10.1145/3460773. Rochman-Sharabi et al. (2025) Rochman-Sharabi, O., Lewin, S., and Louppe, G. A neural material point method for particle-based emulation. Transactions on Machine Learning Research, 2025. ISSN 2835-8856. URL https://openreview.net/forum?id=zSK81A2hxQ. Sanchez-Gonzalez et al. (2020) Sanchez-Gonzalez, A., Godwin, J., Pfaff, T., Ying, R., Leskovec, J., and Battaglia, P. W. Learning to simulate complex physics with graph networks, 2020. URL https://arxiv.org/abs/2002.09405. Sappl et al. (2019) Sappl, J., Seiler, L., Harders, M., and Rauch, W. Deep learning of preconditioners for conjugate gradient solvers in urban water related problems. arXiv preprint, 2019. Shen et al. (2024) Shen, J., Marwah, T., and Talwalkar, A. UPS: Efficiently building foundation models for PDE solving via cross-modal adaptation. Transactions on Machine Learning Research, 2024. ISSN 2835-8856. URL https://openreview.net/forum?id=0r9mhjRv1E. Sulsky et al. (1994) Sulsky, D., Chen, Z., and Schreyer, H. L. A particle method for history-dependent materials. Computer Methods in Applied Mechanics and Engineering, 118(1):179–196, 1994. doi: 10.1016/0045-7825(94)90112-0. Sun et al. (2023) Sun, Z., Yang, Y., and Yoo, S. A neural PDE solver with temporal stencil modeling. arXiv preprint, 2023. Taccari et al. (2022) Taccari, M. L., Nuttall, J., Chen, X., Wang, H., Minnema, B., and Jimack, P. K. Attention u-net as a surrogate model for groundwater prediction. Advances in Water Resources, 163:104169, 2022. ISSN 0309-1708. doi: https://doi.org/10.1016/j.advwatres.2022.104169. URL https://w.sciencedirect.com/science/article/pii/S0309170822000458. Teng & Choromanska (2019) Teng, Y. and Choromanska, A. Invertible autoencoder for domain adaptation. Computation, 7(2):20, 2019. Toshev et al. (2024a) Toshev, A., Galletti, G., Fritz, F., Adami, S., and Adams, N. Lagrangebench: A lagrangian fluid mechanics benchmarking suite. Advances in Neural Information Processing Systems, 36, 2024a. Toshev et al. (2023) Toshev, A. P., Galletti, G., Brandstetter, J., Adami, S., and Adams, N. A. Learning lagrangian fluid mechanics with e(33)-equivariant graph neural networks, 2023. Toshev et al. (2024b) Toshev, A. P., Erbesdobler, J. A., Adams, N. A., and Brandstetter, J. Neural sph: Improved neural modeling of lagrangian fluid dynamics, 2024b. URL https://arxiv.org/abs/2402.06275. Tran et al. (2023) Tran, A., Mathews, A., Xie, L., and Ong, C. S. Factorized fourier neural operators. In The Eleventh International Conference on Learning Representations, 2023. URL https://openreview.net/forum?id=tmIiMPl4IPa. Tskhakaya (2008) Tskhakaya, D. Chapter 6: The particle-in-cell method. In Fehske, H., Schneider, R., and Weiße, A. (eds.), Computational Many-Particle Physics, volume 739 of Lecture Notes in Physics, p. 271–?,. Springer, Berlin, Heidelberg, 2008. ISBN 978-3-540-74685-0. doi: 10.1007/978-3-540-74686-7. Viswanath et al. (2024) Viswanath, H., Chang, Y., Berner, J., Chen, P. Y., and Bera, A. Reduced-order neural operators: Learning lagrangian dynamics on highly sparse graphs. arXiv preprint arXiv:2407.03925, 2024. Wang et al. (2024) Wang, H., Cao, Y., Huang, Z., Liu, Y., Hu, P., Luo, X., Song, Z., Zhao, W., Liu, J., Sun, J., Zhang, S., Wei, L., Wang, Y., Wu, T., Ma, Z.-M., and Sun, Y. Recent advances on machine learning for computational fluid dynamics: A survey, 2024. URL https://arxiv.org/abs/2408.12171. Wang et al. (2025) Wang, H., Li, R., Xu, F., Sun, F., Han, K., Huang, Z., Wan, G., Chang, C., Luo, X., Wang, W., and Sun, Y. Fd-bench: A modular and fair benchmark for data-driven fluid simulation, 2025. URL https://arxiv.org/abs/2505.20349. Wang & Wang (2024) Wang, T. and Wang, C. Latent neural operator for solving forward and inverse pde problems, 2024. URL https://arxiv.org/abs/2406.03923. Williams et al. (2023) Williams, C., Falck, F., Deligiannidis, G., Holmes, C. C., Doucet, A., and Syed, S. A unified framework for u-net design and analysis. In Thirty-seventh Conference on Neural Information Processing Systems, 2023. URL https://openreview.net/forum?id=43ruO2fMjq. Wu et al. (2023) Wu, H., Hu, T., Luo, H., Wang, J., and Long, M. Solving high-dimensional pdes with latent spectral models. In International Conference on Machine Learning, 2023. Wu et al. (2024) Wu, H., Luo, H., Wang, H., Wang, J., and Long, M. Transolver: A fast transformer solver for pdes on general geometries. In International Conference on Machine Learning, 2024. Wu et al. (2025) Wu, Q., Yang, C., Zeng, K., and Bronstein, M. M. Supercharging graph transformers with advective diffusion. In Forty-second International Conference on Machine Learning, 2025. URL https://openreview.net/forum?id=MaOYl3P84E. Xu et al. (2025) Xu, J., Huang, H., Zou, C., Savva, M., Wei, Y., and Chen, W. Hybrid neural-mpm for interactive fluid simulations in real-time, 2025. URL https://arxiv.org/abs/2505.18926. Yadav et al. (2025) Yadav, V., Casel, M., and Ghani, A. Rf-pinns: Reactive flow physics-informed neural networks for field reconstruction of laminar and turbulent flames using sparse data. Journal of Computational Physics, 524:113698, March 2025. Published: 1 March 2025. Ye et al. (2024) Ye, Z., Huang, X., Chen, L., Liu, H., Wang, Z., and Dong, B. PDEformer: Towards a foundation model for one-dimensional partial differential equations. In ICLR 2024 Workshop on AI4DifferentialEquations In Science, 2024. URL https://openreview.net/forum?id=GLDMCwdhTK. Zhong et al. (2023) Zhong, Y., Fukami, K., An, B., et al. Sparse sensor reconstruction of vortex-impinged airfoil wake with machine learning. Theor. Comput. Fluid Dyn., 37:269–287, 2023. doi: 10.1007/s00162-023-00657-y. (74) Zhou, H., Ma, Y., Wu, H., Wang, H., and Long, M. Unisolver: Pde-conditional transformers towards universal neural pde solvers. In Forty-second International Conference on Machine Learning. Appendix A Impact Statement Our neural simulator has the potential to transform computational fluid dynamics (CFD) and the broader sciences by enabling faster, more scalable, and more adaptable modeling tools. The development of approaches that combine data-driven learning with physically grounded representations opens new opportunities for fast simulations in areas such as engineering design, climate modeling, and multiphase flow analysis. By improving the efficiency and robustness of CFD, these methods can accelerate scientific discovery and expand access to large-scale simulations that were previously computationally prohibitive. We do not identify any adverse societal or ethical impacts associated with this paper. Appendix B Limitation and Future Work While this framework significantly enhances the rollout performance of Lagrangian neural simulators, a primary limitation remains the divergence of trajectories over extreme temporal horizons, such as those spanning 100 physical hours. This drift is not unique to our architecture. Rather, it stems from the irreducible local truncation errors inherent in both neural simulators and numerical solvers. As these microscopic errors compound over millions of integration steps, the global truncation error eventually accumulates, leading to a potential breakdown in physical fidelity. This phenomenon represents a fundamental, systemic challenge in numerical analysis and the simulation of dynamical systems. Because error accumulation is an intrinsic property of discrete solvers, achieving stability over infinite horizons remains an open problem for the community. Future research could investigate specialized neural architectures that explicitly regularize long-term stability or incorporate error-correction manifolds. By moving beyond the imitation of classical solvers, there is a distinct opportunity to develop neural simulators that eventually surpass the stability and accuracy of traditional numerical methods. Resolving these foundational issues regarding long-term numerical stability will require a sustained, collective effort from the broader physics-ML community. Fully addressing these systemic error accumulation patterns is beyond the scope of the current work. Appendix C Datasets We provide a brief description of the main benchmark dataset used in this work and refer the readers to Toshev et al. (2024a) for additional details. Taylor Green Vortex The dataset represents a decaying flow system characterized by periodic boundary conditions and a specific initial velocity field that leads to kinetic energy decay through viscous interactions. Dataset generation employs a Lagrangian SPH scheme, solving the weakly-compressible Navier-Stokes equations coupled with a barotropic equation of state. The solver setup involves randomly drawing particle positions, which are then relaxed under periodic boundary conditions via 1000 steps of SPH relaxation to create the initial state. The 2D TGV case is defined on a 1×11× 1 spatial domain with 2,500 particles at a Reynolds number of 100, using a particle spacing (Δx x) of 20×10−320× 10^-3 and a time step (Δt t) of 40×10−340× 10^-3. Conversely, the 3D TGV case utilizes a 2π×2π×2π2π× 2π× 2π domain populated by 8,000 particles at a Reynolds number of 50, with a Δx x of 314.16×10−3314.16× 10^-3 and a Δt t of 500×10−3500× 10^-3. The 2D version follows an analytical solution of exponentially decaying velocity. Lid-Driven Cavity Flow The dataset simulates a flow within a confined domain driven by a moving top wall. This dataset is generated using the Lagrangian SPH solver, which solves the weakly-compressible Navier-Stokes equations. The ground truth setup initializes the fluid velocity to zero and runs the simulation until the system reaches a statistically stationary equilibrium state, at which point data collection begins. While the solver uses a generalized wall boundary condition with multiple layers of dummy particles to enforce impermeability and no-slip conditions, the dataset retains only the innermost layer. In the 2D case, the simulation comprises 2,708 particles within a 1.12×1.121.12× 1.12 spatial domain at a Reynolds number of 100, utilizing a particle spacing (Δx x) of 20×10−320× 10^-3 and a time step (Δt t) of 40×10−340× 10^-355. The 3D case involves 8,160 particles in a 1.25×1.25×0.51.25× 1.25× 0.5 domain, also at a Reynolds number of 100, with a Δx x of 41.667×10−341.667× 10^-3 and a Δt t of 90×10−390× 10^-3. Notably, the 3D version applies periodic boundary conditions in the z-direction to recover the dynamics of the 2D solution. Reverse Poiseuille Flow The dataset simulates a fully periodic flow driven by a spatially varying external force field. This dataset is generated using the Lagrangian SPH solver, which solves the weakly-compressible Navier-Stokes equations. The ground truth setup initializes with zero fluid velocity and applies a force of magnitude 1 in the lower half of the domain and -1 in the upper half, running the simulation until it reaches a statistically stationary equilibrium state before data collection begins. In the 2D case, the system consists of 3,200 particles within a 1×21× 2 spatial domain at a Reynolds number of 10, utilizing a particle spacing (Δx x) of 25×10−325× 10^-3 and a time step (Δt t) of 40×10−340× 10^-3. The 3D counterpart involves 8,000 particles in a 1×2×0.51× 2× 0.5 domain, also at a Reynolds number of 10, with a Δx x of 50×10−350× 10^-3 and a Δt t of 100×10−3100× 10^-35. Although theoretically similar to laminar channel flow, this specific setup at Re=10Re=10 is not fully laminar and exhibits mixing, making the dynamics more diverse for learning tasks. Dam Break The dataset is designed to benchmark the simulation of free-surface flows involving significant topological changes. Generated using the Lagrangian SPH solver to solve the weakly-compressible Navier-Stokes equations, this dataset distinctly calculates density via evolution equations rather than summation to accurately account for the lack of full kernel support at the free surface. The ground truth setup initializes the system by randomly drawing particle positions and letting them relax under non-periodic boundary conditions for 1000 steps. To manage the solid boundaries, the solver employs generalized wall boundary conditions with dummy particles to enforce no-slip and impermeability, though only the innermost particle layer is retained in the final dataset. The 2D simulation is defined on a spatial domain of 5.486×2.125.486× 2.12 and contains 5,740 particles. While this problem is typically modeled as inviscid, this specific implementation includes a small amount of physical viscosity to minimize wall artifacts, resulting in a Reynolds number of 40,000, with a particle spacing (Δx x) of 20×10−320× 10^-3 and a time step (Δt t) of 30×10−330× 10^-3. C.1 Higher Resolution Simulation Data To generate higher-resolution simulation data for the Taylor-Green Vortex, we systematically reduce the initial particle spacing, Δx x, starting from 20×10−320× 10^-3 and progressively decreasing it to 15×10−315× 10^-3, 10×10−310× 10^-3, and 5×10−35× 10^-3. These configurations yield simulations with increasing particle counts of 3,600, 4,489, 10,000, and 40,000, respectively. Reducing the inter-particle spacing enhances the spatial resolution of the Lagrangian discretization, which is critical for numerical accuracy. A denser particle distribution enables more precise approximation of spatial derivatives, resulting in more accurate gradient computations and improved fidelity when querying field variables such as velocity and pressure. Furthermore, this setup serves as a scalability experiment for the neural surrogate models. The significant increase in particle number challenges the model to generalize to larger, more computationally intensive graph structures while maintaining predictive accuracy. C.2 Fluid-Solid Interaction Data The simulation takes place in a rectangular fluid domain with physical dimensions of 18.0×18.018.0× 18.0 units, extending from x=0x=0 to 18.018.0 and y=−9.0y=-9.0 to 9.09.0. A cylinder with a diameter of 1.21.2 units is positioned at coordinates (6.0,0.0)(6.0,0.0) to obstruct the flow. The simulation uses a fixed time step (dtdt) of approximately 0.00270.0027 seconds, calculated initially based on the minimum of CFL and viscous stability criteria, without utilizing adaptive timestepping. The system saves the simulation state every 10 time steps, effectively capturing the dynamics of the flow at a Reynolds number of 200. The dataset is generated using Ramachandran et al. (2021). C.3 Deformable Domain Data The data is generated by simulating the Navier-Stokes equations for an incompressible fluid, modeling the evolution of a circular patch of fluid deformed into an ellipse by an initial linear velocity field. This process perfectly conserves the area of the fluid. The simulation utilizes approximately 5,026 particles at standard resolution and dynamically calculates an adaptive time step of roughly Δt=5.27×10−6 t=5.27× 10^-6 seconds using a CFL constraint. We generate 10000 frames for training, 1000 for validation, and 1000 for testing. The dataset is generated using Ramachandran et al. (2021). Appendix D Model Training and Implementation Details All models are implemented using the PyTorch framework (Paszke et al., 2019) and trained with the Adam optimizer (Kingma & Ba, 2017). All models are configured to directly predict particle positions from unperturbed inputs. We deliberately exclude these auxiliary training strategies found in prior literature (Pfaff et al., 2021; Sanchez-Gonzalez et al., 2020) to isolate the specific contributions of our proposed modules. GNS The model employs a message-passing GNN architecture (Sanchez-Gonzalez et al., 2020), adhering to the hyperparameter configuration detailed in Toshev et al. (2024a). Both the encoder and decoder are composed of 4-layer MLPs with ReLU activations, while the processor utilizes 10 message-passing layers with a latent dimension of 128. SEGNN We implement the geometric architecture proposed by Brandstetter et al. (2021), adopting the hyperparameter guidelines from Toshev et al. (2024a). The network comprises 10 message-passing layers with a latent dimension of 64, using 2-layer MLP blocks throughout. We restrict the maximum spherical harmonic degree (LmaxL_max) to 1 for both hidden representations and attributes. The input features are treated as isotropic. GraphUnet We adopt the hierarchical architecture proposed by Gao & Ji (2019), which features three downsampling and upsampling stages with a pooling ratio of 0.60.6. The encoder and decoder modules are composed of 4-layer MLPs with ReLU activations, mapping to a latent dimension of 128128. At the coarsest resolution (the bottom level), the processor employs a stack of 4 message-passing GNN layers. GraphTransformer We adopt the model architecture proposed in Dwivedi & Bresson (2021). We use 66 layers with a latent dimension of 128128. The encoder and decoder modules are composed of 4-layer MLPs with ReLU activations. AdvDIFFormer "We adopt the model architecture proposed in Wu et al. (2025), configured with a hidden channel size of 128. We set the approximation order K to 1 and utilize 2 attention heads. NeuralMPM We adopt the architecture proposed by Rochman-Sharabi et al. (2025), which operates on a fixed 32×3232× 32 Eulerian grid. The central processor employs a U-Net structure with a base channel dimension of 64. PhysicsNFP We adopt the topology-aware architecture proposed by Jiang et al. (2025). The model comprises 10 layers with a hidden channel dimension of 128128. Ours Our architecture incorporates M=2M=2 downsampling and upsampling stages with a ratio of 0.80.8. The model operates with a latent dimension of 128128 and utilizes a central processor comprising L=4L=4 steps. Both the encoder and decoder modules are constructed as 4-layer MLPs with ReLU activations. The model is trained using a learning rate schedule that decays from 10−410^-4 to 10−610^-6. Appendix E Model Analysis We conduct extensive analysis on the Lid-Driven Cavity Flow dataset to isolate the contribution of each module. Ablation Study on Downsampler: We perform an ablation study to evaluate the effect of the downsampler and upsampler. Results are reported in Figure 3. Across all metrics, the model equipped with these components consistently outperforms the variant without them over the five-step rollout. This observation aligns with our design rationale. In the Lid-Driven Cavity flow dataset, many particles in the lower-left region remain nearly stationary, while significant dynamics occur elsewhere. The downsampler and upsampler encourage the model to focus its capacity on regions with strong flow variation. As a result, the model with these modules achieves notably lower MSE and kinetic energy error. We attribute this to the difference in the processor, as the processor of the ablated model must explicitly model interactions among all particles, whereas our hierarchical framework delegates the coarse forward PDE dynamics to the processor and leaves the fine-scale refinements to the downsampler and upsampler. Figure 3: Five-step rollout results of the ablation study on the downsampler and upsampler, reported on the Lid-Driven Cavity flow dataset. Ours denotes the proposed model, while Ablation refers to the variant without the downsampler and upsampler. Ablation Study on Combined Self–Cross Attention Feature: We conduct an ablation study to examine the role of the combined self–cross attention feature. Specifically, we remove the cross-attention between Lagrangian and Eulerian features and retain only the Lagrangian self-attention, replacing Equation 8 with Equation 9. Results are shown in Figure 4. =,^lag=∈ℝN~lag×dlatent. F= A, q_lag= W\, F N_lag× d_latent. (9) The model that integrates both Lagrangian self-attention and Lagrangian–Eulerian cross-attention consistently outperforms the ablated variant. This highlights the importance of cross-attention in enhancing overall performance by leveraging information from both representations. Figure 4: Five-step rollout results of the ablation study on the combined self-cross attention feature, reported on the Lid-Driven Cavity flow dataset. Ours denotes the proposed model, while Ablation refers to the variant without cross attention. Ablation Study on Aggregating Lagrangian Nodes to Eulerian Nodes: We also examine the role of Eulerian features in defining aggregation points for Lagrangian features. In the full model, downsampled Eulerian nodes serve as anchors for aggregating Lagrangian representations. To assess their contribution, we design an ablation study where these nodes are replaced by a learnable alternative: an MLP that takes the positions of downsampled Lagrangian particles as input and outputs logits indicating the assignment of each particle to aggregation points. The number of aggregation points is kept equal to the number of downsampled Eulerian nodes, and the logits determine the contribution of each particle to its assigned aggregation node. Results are reported in Figure 5. The full model that aggregates Lagrangian features onto downsampled Eulerian nodes consistently achieves stronger performance across all metrics. The full model that aggregates Lagrangian features onto downsampled Eulerian nodes consistently achieves stronger performance across all metrics. We attribute these gains to two factors. First, the downsampled Eulerian nodes naturally share a spatial distribution similar to that of the downsampled Lagrangian particles. Aggregating Lagrangian features onto this Eulerian structure produces representations that remain dynamically consistent with the Eulerian features, which is particularly beneficial during cross-attention, where the Eulerian branch provides a global perspective while the Lagrangian branch supplies fine-scale details. Second, the MLP-based aggregation points discard explicit spatial structure, weakening the self-attention’s ability to resolve PDE dynamics at a coarse level and diminishing the effectiveness of cross-attention in aligning complementary representations. Figure 5: Five-step rollout results of the ablation study on the role of Eulerian feature in Lagrangian aggregation, reported on the Lid-Driven Cavity flow dataset. Ours denotes the proposed model, while Ablation refers to the variant that uses an MLP to determine aggregation points. Analysis on Eulerian Discretization: The Eulerian features are represented on a discretized grid with NeuN_eu nodes. We investigate how the discretization resolution affects model performance using the Lid-Driven Cavity flow dataset, defined on a bounded 1.12×1.121.12× 1.12 domain. Each direction is discretized with s nodes, and the results are shown in Figure 6. Among the tested resolutions, s=32s=32 yields the best performance. Extremely coarse grids such as s=8s=8, or overly fine grids such as s=64s=64, perform the worst. With too few Eulerian nodes, the representation lacks sufficient information about the global fluid field. In addition, aggregating Lagrangian features onto such a small set of nodes leads to significant information loss, as particles from diverse locations and flow behaviors are forced to map to the same averaged values. Figure 6: Five-step rollout results of varying the discretization resolution of Eulerian nodes, reported on the Lid-Driven Cavity flow dataset. Analysis on Downsampling Ratio: We study the impact of the downsampling ratio, λ, on model performance, with results shown in Figure 7. The best accuracy is obtained at λ=0.8λ=0.8, while more aggressive choices such as λ=0.2λ=0.2 or λ=0.4λ=0.4 lead to a sharp decline in performance. This behavior is expected, as the purpose of downsampling is to merge particles that move coherently or remain nearly static, since these contribute little additional information when considered separately. However, if the ratio is too aggressive, the model also discards particles carrying essential information about the fluid dynamics, which in turn degrades the solution. Figure 7: Five-step rollout results of varying the downsampling ratio, reported on the Lid-Driven Cavity flow dataset. Analysis on Number of Downsampling Layers: We analyze how the number of downsampling layers, M, influences model performance, with results shown in Figure 8. The best results are obtained with M=2M=2, followed by M=1M=1, while deeper configurations with M=3M=3 or M=4M=4 lead to a significant drop in performance. These findings are consistent with the results from the study on the downsampling ratio. Increasing the number of downsampling layers progressively reduces the particle set, which can cause the loss of particles carrying critical information about the fluid dynamics. Interestingly, two layers outperform a single layer. We attribute this to the hierarchical structure. Applying downsampling twice enables the model to better identify and preserve important particles, while additional rounds of message passing improve the transfer of information from distant particles to the aggregation nodes. Figure 8: Five-step rollout results of varying the number of downsampling layers, reported on the Lid-Driven Cavity flow dataset. Appendix F Additional Baselines on Eulerian Mesh Simulator Since Eulerian mesh-based neural simulators operate on graph structures analogous to our Lagrangian framework, they are inherently compatible with our setting. Consequently, we evaluate our approach against a broad suite of state-of-the-art Eulerian mesh models, including the Graph Kernel Operator (GKO) (Li et al., 2020), General Neural Operator Transformer (GNOT) Hao et al. (2023), Directional Transport-Aware Graph Neural Network (DTA-GNN) (Li et al., 2025), Graph Interaction Operator for Reduced-Order Modeling (GIOROM) (Viswanath et al., 2024), Geometry Informed Neural Operator (GINO) (Li et al., 2023a), Transolver (Wu et al., 2024), and Transolver++ (Luo et al., 2025). The experiments are conducted on the Lid-Driven Cavity Flow dataset to assess performance. Metric GKO GNOT DTA-GNN GIOROM GINO Transolver Transolver++ Ours MSEMSE 5.205×10−45.205× 10^-4 1.011×10−51.011× 10^-5 3.505×10−63.505× 10^-6 8.335×10−68.335× 10^-6 3.601×10−43.601× 10^-4 1.308×10−51.308× 10^-5 2.665×10−52.665× 10^-5 1.868×10−61.868× 10^-6 KEKE 1.161×10−41.161× 10^-4 9.113×10−79.113× 10^-7 5.127×10−75.127× 10^-7 1.765×10−61.765× 10^-6 8.656×10−58.656× 10^-5 1.376×10−61.376× 10^-6 1.417×10−51.417× 10^-5 7.240×10−87.240× 10^-8 SinkhornSinkhorn 6.609×10−66.609× 10^-6 6.470×10−76.470× 10^-7 1.560×10−61.560× 10^-6 1.245×10−61.245× 10^-6 7.321×10−67.321× 10^-6 2.001×10−62.001× 10^-6 2.693×10−62.693× 10^-6 4.675×10−74.675× 10^-7 Table 7: Comparison of MSE, Kinetic Energy (KE) error, and Sinkhorn divergence on Lid-Driven Cavity flow dataset. The relative performance gain against the strongest baseline is computed as Baseline−OursBaseline×100% Baseline-OursBaseline× 100\%. While several baselines incorporate aggregation or downsampling techniques, our method introduces fundamental improvements in spatial awareness and domain flexibility. For instance, GINO projects particles onto fixed grids to employ Fourier Neural Operators (FNO), a reliance that renders it inapplicable to irregular domains where our method excels. Furthermore, while Transolver and Transolver++ aggregate points into physical tokens based on latent states, our approach leverages physical coordinates to explicitly preserve the underlying spatial structure. Similarly, although GIOROM generates sparse representations of the input graph, it relies on static heuristics such as random sampling or farthest-point selection. In contrast, our downsampler dynamically identifies and selects salient features based on the data, ensuring the preservation of critical information and the production of more meaningful sparse representations. The results are summarized in Table 7. We observe that our proposed model achieves significantly superior performance compared to the Eulerian mesh baseline simulators. This highlights the effectiveness of our design choices, particularly regarding sparse representation construction and particle aggregation techniques. Appendix G Additional Experiments on Diverse Materials Although this work focuses primarily on Lagrangian fluid simulation, we also evaluate our model’s performance on diverse materials, including sand, goop, and multi-material interactions. Using datasets from Sanchez-Gonzalez et al. (2020), we report performance across 200 rollout steps. The results, summarized in Table 8, demonstrate that our model generalizes effectively to these alternative materials. We exclude PhysicsNFP from this comparison as it is unable to accommodate varying particle counts. Steps Data Model Metric 25 50 75 100 150 200 MSEMSE 1.668×10−41.668× 10^-4 1.616×10−31.616× 10^-3 5.883×10−35.883× 10^-3 2.173×10−22.173× 10^-2 2.421×10−22.421× 10^-2 3.240×10−23.240× 10^-2 KEKE 6.536×10−36.536× 10^-3 2.414×10−12.414× 10^-1 8.072×1008.072× 10^0 4.850×1014.850× 10^1 3.271×1013.271× 10^1 7.792×1017.792× 10^1 GraphUnet SinkhornSinkhorn 3.095×10−43.095× 10^-4 3.849×10−33.849× 10^-3 9.788×10−39.788× 10^-3 3.947×10−23.947× 10^-2 3.078×10−23.078× 10^-2 3.099×10−23.099× 10^-2 MSE 1.321×10−41.321× 10^-4 9.191×10−49.191× 10^-4 3.458×10−33.458× 10^-3 9.409×10−39.409× 10^-3 9.422×10−39.422× 10^-3 1.367×10−21.367× 10^-2 KE 3.784×10−33.784× 10^-3 2.153×10−12.153× 10^-1 4.183×1004.183× 10^0 3.370×1013.370× 10^1 2.775×1012.775× 10^1 4.773×1014.773× 10^1 Goop Ours Sinkhorn 2.599×10−42.599× 10^-4 1.847×10−31.847× 10^-3 6.752×10−36.752× 10^-3 1.838×10−21.838× 10^-2 9.632×10−39.632× 10^-3 1.938×10−21.938× 10^-2 MSEMSE 5.578×10−55.578× 10^-5 3.152×10−33.152× 10^-3 1.378×10−21.378× 10^-2 2.729×10−22.729× 10^-2 4.462×10−24.462× 10^-2 5.060×10−25.060× 10^-2 KEKE 7.982×10−47.982× 10^-4 1.145×1001.145× 10^0 2.373×1012.373× 10^1 9.424×1019.424× 10^1 2.675×1022.675× 10^2 3.590×1023.590× 10^2 GraphUnet SinkhornSinkhorn 9.766×10−59.766× 10^-5 5.228×10−35.228× 10^-3 2.567×10−22.567× 10^-2 5.222×10−25.222× 10^-2 8.543×10−28.543× 10^-2 9.783×10−29.783× 10^-2 MSE 4.545×10−54.545× 10^-5 2.913×10−32.913× 10^-3 1.346×10−21.346× 10^-2 2.429×10−22.429× 10^-2 2.348×10−22.348× 10^-2 2.05×10−22.05× 10^-2 KE 6.196×10−46.196× 10^-4 9.243×10−19.243× 10^-1 2.153×1012.153× 10^1 7.226×1017.226× 10^1 7.446×1017.446× 10^1 5.995×1015.995× 10^1 Sand Ours Sinkhorn 8.544×10−58.544× 10^-5 4.312×10−34.312× 10^-3 2.536×10−22.536× 10^-2 4.572×10−24.572× 10^-2 4.102×10−24.102× 10^-2 3.477×10−23.477× 10^-2 MSEMSE 1.909×10−41.909× 10^-4 9.311×10−49.311× 10^-4 2.244×10−32.244× 10^-3 3.854×10−33.854× 10^-3 8.534×10−38.534× 10^-3 1.770×10−21.770× 10^-2 KEKE 2.270×10−22.270× 10^-2 7.012×10−17.012× 10^-1 4.394×1004.394× 10^0 1.283×1011.283× 10^1 3.196×1013.196× 10^1 1.369×1021.369× 10^2 GraphUnet SinkhornSinkhorn - - - - - - MSE 1.773×10−41.773× 10^-4 7.528×10−47.528× 10^-4 1.571×10−31.571× 10^-3 2.490×10−32.490× 10^-3 7.056×10−37.056× 10^-3 1.259×10−21.259× 10^-2 KE 1.923×10−21.923× 10^-2 4.449×10−14.449× 10^-1 1.981×1001.981× 10^0 4.214×1004.214× 10^0 2.013×1012.013× 10^1 6.836×1016.836× 10^1 Multi Material Ours Sinkhorn - - - - - - Table 8: Comparison of MSE, Kinetic Energy (KE) error, and Sinkhorn divergence across datasets containing diverse materials. Sinkhorn divergence is omitted for the multi-material dataset because it does not account for material type. Figure 9: Analysis of the runtime and memory usage of models. Appendix H Computational Complexity We evaluate the runtime and memory consumption of different models. Although our approach processes both Lagrangian and Eulerian representations, it is naturally amenable to parallelization at inference time. In particular, the encoding and downsampling stages for the two representations can be executed concurrently, which reduces computational time. As illustrated in Figure 9, this design results in only a marginal increase in inference time compared to the baselines. We acknowledge that our model requires additional memory during inference, which could pose a limitation in certain scenarios. Nonetheless, we argue that the overhead remains modest, well within the capacity of widely available commercial GPUs. Appendix I Hardware Specification We implement models in PyTorch (Paszke et al., 2019). All experiments can be run on servers/workstations with the following configuration: • 80 CPUs, 503G Mem, 8 x NVIDIA V100 GPUs. • 48 CPUs, 220G Mem, 8 x NVIDIA TITAN XP GPUs. • 96 CPUs, 1.0T Mem, 8 x NVIDIA A100 GPUs. • 64 CPUs, 1.0T Mem, 8 x NVIDIA RTX A6000 GPUs. • 224 CPUs, 1.5T Mem, 8 x NVIDIA L40S GPUs. • 128 CPUs, 480G Mem, 8 × NVIDIA RTX 4090 GPUs.