E($3$) Equivariant Graph Neural Networks for Particle-Based Fluid
  Mechanics
We contribute to the vastly growing field of machine learning for engineering
systems by demonstrating that equivariant graph neural networks have the
potential to learn more accurate dynamic-interaction models than their
non-equivariant counterparts. We benchmark two well-studied fluid flow systems,
namely the 3D decaying Taylor-Green vortex and the 3D reverse Poiseuille flow,
and compare equivariant graph neural networks to their non-equivariant
counterparts on different performance measures, such as kinetic energy or
Sinkhorn distance. Such measures are typically used in engineering to validate
numerical solvers. Our main findings are that while being rather slow to train
and evaluate, equivariant models learn more physically accurate interactions.
This indicates opportunities for future work towards coarse-grained models for
turbulent flows, and generalization across system dynamics and parameters.
Artur P. Toshev, Gianluca Galletti, Johannes Brandstetter, Stefan Adami, Nikolaus A. Adams
# E(3) Equivariant Graph Neural Networks for Particle-Based Fluid Mechanics

###### Abstract

We contribute to the vastly growing field of machine learning for engineering systems by demonstrating that equivariant graph neural networks have the potential to learn more accurate dynamic-interaction models than their non-equivariant counterparts. We benchmark two well-studied fluid flow systems, namely the 3D decaying Taylor-Green vortex and the 3D reverse Poiseuille flow, and compare equivariant graph neural networks to their non-equivariant counterparts on different performance measures, such as kinetic energy or Sinkhorn distance. Such measures are typically used in engineering to validate numerical solvers. Our main findings are that while being rather slow to train and evaluate, equivariant models learn more physically accurate interactions. This indicates opportunities for future work towards coarse-grained models for turbulent flows, and generalization across system dynamics and parameters.

## 1 Particle-based fluid mechanics

Navier-Stokes equations (NSE) are omnipresent in fluid mechanics, hydrodynamics or weather modeling. However, for the majority of problems, solutions are analytically intractable, and obtaining accurate predictions necessitates falling back to numerical solution schemes. Those can be split into two categories: grid/mesh-based (Eulerian description) and particle-based (Lagrangian description).

**Smoothed Particle Hydrodynamics.** In this work, we investigate Lagrangian methods and more precisely the Smoothed Particle Hydrodynamics (SPH) approach, which was independently developed by Gingold & Monaghan (1977) and Lucy (1977) to simulate astrophysical systems. Since then, SPH has established as the preferred approach in various applications ranging from free surfaces such as ocean waves (Violeau & Rogers, 2016), through fluid-structure interaction systems (Zhang et al., 2021), to selective laser melting in additive manufacturing (Weiraher et al., 2019).

Figure 1: Velocity magnitude of Taylor-Green vortex (a) and x-velocity of reverse Poiseuille (b).
The main idea behind SPH is to represent the fluid properties at discrete points in space and to use truncated radial interpolation kernel functions to approximate them at any arbitrary location. The kernel functions are used to estimate state statistics which define continuum-scale interactions between particles. The justification for truncating kernel support is the assumption of local interactions between particles. The resulting discretized equations are then integrated in time using numerical integration techniques like symplectic Euler by which the particle positions are updated.

To generate training data, we implemented our own SPH solver based on the transport velocity formulation by Adamir et al. (2013), which promises a homogeneous particle distribution over the domain. We then selected two flow cases, both of which are well-known in the fluid mechanics community: the 3D laminar Taylor-Green Vortex and the 3D reverse Poiseuille Flow. We are planning to open-source the datasets in the near future.

Taylor-Green Vortex.The Taylor-Green vortex system (TGV, see Figure 1 (a)) with Reynolds number of \(\text{Re}=100\) is neither laminar nor turbulent, i.e. there is no layering of the flow (typical for laminar flows), but also the small scales caused by vortex stretching do not lead to a fully developed energy cascade (typical for turbulent flows) Brachet et al. (1984). The TGV has been extensively studied starting with Taylor & Green (1937) and continuing all the way to Sharma & Sengupta (2019). The TGV system is typically initialized with a velocity field given by

\[u=-\cos(kx)\cos(ky)\cos(kz)\;,\qquad v=\sin(kx)\cos(ky)\cos(kz)\;,\qquad w=0\;, \tag{1}\]

where \(k\) is an integer multiple of \(2\pi\). The TGV datasets used in this work consist of 8/2/2 trajectories for training/validation/testing, where each trajectory comprises 8000 particles. Each trajectory spans 1s physical time and was simulated with \(dt=0.001\) resulting in 1000 time steps per trajectory. The ultimate goal would be to learn the dynamics over much larger time steps than those taken by the numerical solver, but with this dataset we just want to demonstrate the applicability of learned approaches to reproducing numerical solver results.

Reverse Poiseuille Flow.The Poiseuille flow, i.e. laminar channel flow, is another well-studied flow case in fluid mechanics. However, channel flow requires the treatment of wall-boundary conditions, which is beyond the focus of this work. In this work, we therefore consider data obtained by reverse Poiseuille flow (RPF, see Figure 1 (b)) (Fedosov et al., 2008), which essentially consists of two opposing streams in a fully periodic domain. Those flows are exposed to opposite force fields, i.e., the upper and lower half are accelerated in negative \(x\) direction and positive \(x\) direction, respectively. Due to the fact that the flow is statistically stationary (the vertical velocity profile has a time-independent mean value), the RPF dataset consists of one long trajectory spanning 120s. The flow field is discretized by 8000 particles and simulated with \(dt=0.001\), followed by sub-sampling at every 10th step. Learning to directly predict every 10th state is what we call temporal coarse-graining. The resulting number of training/validation/testing instances is the same as for TGV, namely 8000/2000/2000.

## 2 (Equivariant) graph network-based simulators

We first formalize the task of autoregressively predicting the next state of a Lagrangian fluid mechanics simulation based on the notation from Sanchez-Gonzalez et al. (2020). Let \(X^{t}\) denote the state of a particle system at time \(t\). One full trajectory of \(K+1\) steps can be written as \(\mathbf{X}^{t_{0:K}}=(X^{t_{0}},\dots,X^{t_{K}})\). Each state \(\mathbf{X}^{t}\) is made up of \(N\) particles, namely \(\mathbf{X}^{t}=(\mathbf{x}^{t}_{1},\mathbf{x}^{t}_{2},\dots\mathbf{x}^{t}_{N})\), where each \(\mathbf{x}_{i}\) is the state vector of the \(i\)-th particle. However, the inputs to the learned simulator can span multiple time instances. Each node \(\mathbf{x}^{t}_{i}\) can contain node-level information like the current position \(\mathbf{p}^{t}_{i}\) and a time sequence of \(H\) previous velocity vectors \(\hat{\mathbf{p}}^{t_{i+k-H:k}}\), as well as global features like the external force field \(\mathbf{f}_{i}\) in the reverse Poiseuille flow. To build the connectivity graph, we use an interaction radius of \(\sim 1.5\) times the average interparticle distance. This results in around 10-20 one-hop neighbors.

Graph Network-based Simulator.The Graph Network-based Simulator (GNS) framework (Sanchez-Gonzalez et al., 2020) is one of the most popular learned surrogates for engineering particle-based simulations. The main idea of the GNS model is to use the established encoder-processor-decoder architecture (Battaglia et al., 2018) with a processor that stacks several message passing layers (Gilmer et al., 2017). One major strength of the GNS model lies in its simplicity given that all its building blocks are simple MLPs. However, the performance of GNS when predicting
long trajectories strongly depends on choosing the right amount of Gaussian noise to perturb input data. Additionally, GNS and other non-equivariant models are less data-efficient (Batzner et al., 2022). For these reasons, we implement and tune GNS as a comparison baseline, and use it as an inspiration for which setup, features, and hyperparameters to use for equivariant models.

Steerable E(3)-equivariant Graph Neural Network.Steerable E(3)-equivariant Graph Neural Networks (SEGNNs) (Brandstetter et al., 2022) are an instance of E(3)-equivariant GNNs, i.e., GNNs that are equivariant with respect to isometries of the Euclidean space (rotations, translations, and reflections). Most E(3)-equivariant GNNs that are tailored towards molecular property prediction tasks, (Batzner et al., 2022; Batatia et al., 2022) restrict the parametrization of the Clebsch-Gordan tensor products to an MLP-parameterized embedding of pairwise distances. In contrast, SEGNNs use general steerable node and edge attributes which can incorporate any kind of physical quantity, and directly learn the weights of the Clebsch-Gordan tensor product. Indeed, extensions of methods such as NequIP (Batzner et al., 2022) towards general physical features would results in something akin to SEGNN.

Steerable attributes strongly impact the Clebsch-Gordan tensor products, and thus finding physically meaningful edge and node attributes is crucial for good performance. In particular, we chose edge attributes \(\hat{a}_{ij}=V(\mathbf{p}_{ij})\), where \(V(\cdot)\) is the spherical harmonic embedding and \(\mathbf{p}_{ij}=\mathbf{p}_{i}-\mathbf{p}_{j}\) are the pairwise distances. We further choose node attributes \(\hat{a}_{i}=V(\tilde{\mathbf{p}}_{i})+\sum_{k\in\mathcal{N}(i)}\hat{a}_{ik}\), where \(\tilde{\mathbf{p}}_{i}\) are averaged historical velocities and \(\mathcal{N}(i)\) is the \(i\)-neighborhood. As for node and edge features, we found that concatenated historical velocities for the nodes and pairwise displacements for the edges capture best the Navier-Stokes dynamics.

For training SEGNNs, we verified that adding Gaussian noise to the inputs (Sanchez-Gonzalez et al., 2020) indeed significantly improves performance. We further found that explicitly concatenating the external force vector \(\mathbf{f}_{i}\) to the node features boosts performance in the RPF case. However, adding \(\mathbf{f}_{i}\) to the node attributes \(\hat{a}_{i}t=V(\mathbf{f}_{i})+V(\tilde{\mathbf{p}}_{i})+\sum_{k\in\mathcal{N }(i)}\hat{a}_{ik}\) does not improve performance.

Other models, like EGNN by Satorras et al. (2021), achieve equivariance by working with invariant messages, but it does not allow the same flexibility in terms of features. On a slightly more distant note, there has been a rapid raise in physics-informed machine learning (Raissi et al., 2019) and operator learning (Li et al., 2021), where functions or surrogates are learned in an Eulerian (grid-based) way. SEGNN is a sound choice for Lagrangian fluid mechanics problems since it is designed to work directly with vectorial information and particles.

## 3 Results

The task we train on is the autoregressive prediction of accelerations \(\tilde{\mathbf{p}}\) given the current position \(\mathbf{p}_{i}\) and \(H=5\) past velocities of the particles. We measured the performance of the GNS and the SEGNN models in four aspects when evaluating on the test dataset: (i) _Mean-squared error_ (MSE) of particle positions \(\text{MSE}_{p}\) when rolling out a trajectory over 100 time steps (1 physical second for both flow cases). This is also the validation loss during training. (ii) _Sinkhorn distance_, as an optimal transport distance measure between particle distributions. Lower values indicate that the particle distribution is closer to the reference one. (iii) _Kinetic energy_\(E_{kin}\) (\(=0.5mv^{2}\)) as a global measure of physical behavior. Performance comparisons are summarized in Table 1. GNS and SEGNN models have roughly the same number of parameters for Taylor-Green (both have 5 layers and 128-dim features), whereas for reverse Poiseuille SEGNN has three times less parameters than GNS (SEGNN has 64-dim features). Looking at the timing in Table 1, equivariant models of similar size are one order of magnitude slower than non-equivariant ones. This is a known result and is related to the constraint of how the Clebsch-Gordan tensor product can be implemented on accelerators like GPUs.

Taylor-Green Vortex.One of the major challenges of the Taylor-Green dataset are the varying input and output scales throughout a trajectory, by up to one order of magnitude. Consequently, this results in the larger importance of first-time steps in the loss even after data normalization. Figure 2 (a) summarizes the most important performance properties of the Taylor-Green vortex experiment. In general, both models match the ground truth kinetic energy well, but GNS drifts away from the reference SPH curve earlier. Both learned solvers, seem to preserve larger system velocities resulting in higher \(E_{kin}\). The rollout MSE for this case matches the behavior seen in \(E_{kin}\).
**Reverse Poiseuille Flow.** The challenge of the reverse Poiseuille case lies in the different velocity scales between the main flow direction (\(x\)-axis) and the \(y\) and \(z\) components of the velocity. Although such unbalanced velocities are used as inputs, target accelerations in \(x\)-, \(y\)-, and \(z\)-direction all underlie similar distributions. This, combined with temporal coarsening makes the problem sensitive to input deviations. Figure 2 (b) shows that SEGNN reproduces the particle distribution almost perfectly, whereas GNS shows signs of particle clustering, resulting in a larger Sinkhorn distance. Interestingly, the shear layers in-between the inverted flows (around planes \(y=\{0,1,2\}\)) seem to have the largest deviation from the ground truth, which could be source of clusters, see Figure 3.

## 4 Future Work

In this work, we demonstrate that equivariant models are well suited to capture underlying physics properties of particle-based fluid mechanics systems. Natural future steps are enforcing physical behaviors such as homogeneous particle distributions, and including recent developments for neural PDE training into the training procedure of Sanchez-Gonzalez et al. (2020). The latter include e.g., the push-forward trick and temporal bundling (Brandstetter et al., 2022). One major weakness of recursively applied solvers, which these strategies aim to mitigate, is error accumulation, which in most cases leads to out-of-distribution states, and consequently unphysical behavior after several rollout steps. We conjecture that together with such extensions equivariant models offer a promising direction to tackle some of the long-standing problems in fluid mechanics, such as the learning of coarse-grained representations of turbulent flow problems, e.g. Taylor-Green (Brachet et al., 1984), or learning the multi-resolution dynamics of NSE problems (Hu et al., 2017).

\begin{table}
\begin{tabular}{c|c c c|c c c}  & \multicolumn{2}{c|}{Taylor-Green vortex} & \multicolumn{4}{c}{Reverse Poiseuille flow} \\ \cline{2-7}  & SEGNN & GNS & SPH & SEGNN & GNS & SPH \\ \hline \(\text{MSE}_{\mathbf{p}}\) & 7.7e-5 & 1.3e-4 & - & 7.7e-3 & 8.0e-3 & - \\ \(\text{MSE}_{E_{kin}}\) & 5.3e-5 & 1.3e-4 & - & 2.8e-1 & 3.0e-1 & - \\ \hline Sinkhorn & 1.3e-7 & 1.1e-7 & - & 7.8e-8 & 1.9e-6 & - \\ Time [ms] & 290 & 32 & 9.7 & 180 & 33 & 110 \\ \(\#\) params & 720k & 630k & - & 180k & 630k & - \\ \end{tabular}
\end{table}
Table 1: Performance measures on the Taylor-Green vortex and reverse Poiseuille flow. The Sinkhorn distance is averaged over test rollouts, the inference time is obtained for one rollout step of 8000 particles.

Figure 2: Taylor-Green vortex (a) and reverse Poiseuille (b) performance evolution.