Year: 2020

Abstract:
We present a data-driven model for fluid simulation under Lagrangian representation. Our model uses graphs to describe the fluid field, where physical quantities are encoded as node and edge features. Instead of directly predicting the acceleration or position correction given the current state, we decompose the simulation scheme into separate parts advection, collision, and pressure projection. For these different reasoning tasks, we propose two kinds of graph neural network structures, node-focused networks, and edge-focused networks. By introducing physics prior knowledge, our model can be efficient in terms of training and inference. Our tests show that the learned model can produce accurate results and remain stable in scenarios with a large number of particles and different geometries. Unlike many previous works, further tests demonstrate that our model is able to retain many important physical properties of incompressible fluids, such as minor divergence and reasonable pressure distribution. Additionally, our model can adopt a range of time step sizes different from ones using in the training set, which indicates its robust generalization capability.

heading: 1 INTRODUCTION
text: For many science and engineering problems, fluids are an essential integral part. How to simulate fluid dynamics accurately has long been studied by researchers and a large class of numerical models have been developed. However, computing high-quality fluid simulation is still computationally expensive despite the advances in computing power. Also, the time of calculation usually increases drastically when the resolution of the simulating scene scales up. A common way to alleviate computing costs is using a data-driven model. Recent progress in the machine learning domain opens up the possibility of employing learning algorithms to learn and model fluid dynamics.
In this paper, we propose a graph-based data-driven fluid dynamics model (Fluid Graph Networks, FGN), which consists of simple multi-layer perceptron and graph inductive architectures. Our model predicts and integrates forward the movement of incompressible fluids based on observations. Compared to previous works in this domain (Ummenhofer et al., 2020; Sanchez-Gonzalez et al., 2020), our model enjoys traceability of physical properties of the system, like low velocity-divergence and constant particle density, and it can predict reasonable pressure distribution. Experiments demonstrate that our model can remain stable and accurate in long-term simulation. Although our model is entailed and customized for fluid simulation, it can be extended to simulation of other dynamics under the Lagrangian framework, as it takes universal features (positions, velocities, particle density) under the Lagrangian framework as input.

heading: 2 RELATED WORKS
text: Our model is built upon the Lagrangian representation of fluid, where continuous fluids are discretized and approximated by a set of particles. The most prominent advantage of the Lagrangian method is that the particle boundary is the material interface, which makes boundary conditions easy to impose, especially when the material interface is large and changing violently. A well-known Lagrangian method is Smooth Particle Hydrodynamics (SPH)(Monaghan, 1988). SPH and its variants are widely used in the numerical physic simulation, especially fluid dynamics under various environments. Particle-based fluid simulation (Müller et al., 2003) introduces SPH model to simulate fluids and
generate realistic visual effects. Moving particle semi-implicit method (MPS) (Koshizuka and Oka, 1996) markedly improves the accuracy and stability of incompressible fluid simulation by introducing a pressure projection procedure that emulates Eulerian grid-based methods. Weakly compressible SPH (WSPH)(Becker and Teschner, 2007) introduces equation of state to model the pressure during the simulation. Predictive-corrective incompressible SPH (Solenthaler and Pajarola, 2009) and divergence-free SPH (Bender and Koschier, 2015) use iterative method to improve the accuracy of pressure calculation in incompressible flow simulation.
Modeling fluid dynamics in a data-driven way has been explored and studied by many researchers. With advances in machine learning algorithms, many data-driven models employing machine learning algorithms have been built. Ladický et al. (2015) reformulate the Navier-Stokes equation as a regression problem and build a regressor using random forest, which significantly improves the calculation efficiency. Tompson et al. (2016), Xiao et al. (2020) learn the pressure projection under the Eulerian framework with a convolutional neural network, which accelerates the fluid simulation. Wiewel et al. (2018) bring significant speed-up by learning a reduced-order representation and predicting the pressure field with an LSTM-based model. Morton et al. (2018) learn the dynamics of airflow around a cylinder based on Koopman theory. de Avila Belbute-Peres et al. (2020) predict fluid flow by combining grid-based method with graph convolutional neural networks.
Learning and reasoning particle dynamics under graph representation has the following benefits and conveniences. First, particle-based methods model physics phenomena as interactions between particles within a local area. This imposes an inductive bias for learning under the Lagrangian framework: dynamics have a strong locality. The locality of unstructured data under Lagrangian representation can be captured by aggregation operation on graphs, such as GCN and other variants (Kipf and Welling, 2016; Hamilton et al., 2017). Second, unlike Eulerian grid-based methods, Lagrangian particle-based methods do not have explicit and structured grid, which makes standard Convolutional Neural Network (CNN) cannot be directly applied to particles without feature processing (Wang et al., 2018; Ummenhofer et al., 2020). Third, many dynamics are based on pairwise relation between particles, like collision, which can be easily interpreted as edge attributes of a graph. Given these factors, recently there have been a rich class of works that use graph neural networks (Scarselli et al., 2009) to learn and reason about underlying physics of interacting objects and particles. (Battaglia et al., 2016; Chang et al., 2016; Sanchez-Gonzalez et al., 2018; Li et al., 2018; Mrowca et al., 2018)

heading: 3 MODEL

heading: 3.1 FLUID DYNAMICS
text: The governing equation for incompressible fluids is the Navier-Stokes equation and the continuity equation as follows (Batchelor, 2000):
Du Dt = −∇p ρ + ν∇2u + g, (1)
∇ · u = 0. (2) To describe the fluid field, there are two kinds of systems, Eulerian and Lagrangian ones. In this work, we adopt a Lagrangian system. A common method to solve the Navier-Stokes equation and discretize fluids under the Lagrangian framework is Smooth Particle Hydrodynamics (SPH) method (Monaghan, 1988), where physical quantities at an arbitrary point in the space are approximated by the states of nearby particles.
In SPH, an arbitrary scalar (or vector) field A (r) at location r can be represented by a convolution:
A (r) = ∫ A (r′)W (|r− r′| , h) dV (r′) , (3)
where W is weighting function or smooth kernel as defined in SPH, h is the smoothing length, which defines the range of particles to be considered and V (r′) is the volume at r. Numerically, the interpolation can be approximated by replacing the integration with a summation.
Based on this model, equation equation 1 and equation 2 can be discretized. The discrete equation system is usually solved under a predictor-corrector scheme, prediction based on advection and correction based on physical properties (such as divergence-free constraint).

heading: 3.2 MODEL
text: Fluids are time-dependent dynamical systems, where location of particles, r, is described by equation of form: dr/dt = f(r). When building a data-driven model to learn and solve this system, we assume the system is Markovian, that is, the state of the system at time step n+ 1 depends only on the state of the previous time step n. The update mechanism in our model can be represented as:
{xn+1,vn+1} = Gθ ({xn,vn}) . (4)
Here {xn,vn} denotes the positional information and velocity of fluid field at time step n. Datadriven model Gθ, parameterized by θ , maps the state of time step n to time step n+ 1.
In order to build a robust and accurate data-driven model, the structure of our model is physicinformed, which enables the model to give interpretable output without losing many physical properties of the system. In general, our model mimics the predictor-corrector scheme and includes three parts, advection net, collision net, and pressure net. They can be divided into two types of graph networks (GN) according to the network structure (Battaglia et al., 2018). Specifically, advection net and pressure net are node-focused graph networks, while collision net is edge-focused networks. As each of these networks has a specific task and different output, they are trained on different data separately.
Node-focused Graph Network Advection net is responsible for the prediction of advection effect and pressure net is responsible for pressure projection. Considering a particle i, the node-focused graph network first aggregates node features from neighbor particles {vj |∀j ∈ N (i)} and output node embedding fi. The embedding fi will then be passed to a processor gR. gR will predict the desirable physical quantities oi (i.e. acceleration a in advection net and pressure p in pressure net).
The whole message passing procedure can be defined as:
oi = gR (gA (Vi)) , Vi = {vj}∀j∈N (i). (5)
Edge-focused Graph Network To prevent particle penetration and increase model stability, we propose a graph network model that is responsible for predicting the effect of collision. As the relative position and relative velocity will have different signs with different observation perspective (i.e. relative velocity vij = −vji), thus the graph in collision net is directed. In collision net, relative features, eij (relative positions, relative velocities between particle i and j), are passed to processor fR as edge features. The processor will output the edge embedding rij between each pair of nodes. Lastly, edge embedding rij is aggregated via aggregator fA, to gather the influence from all nearby particles and predict an overall effect oi on the center particle i. The whole process is defined as:
oi = fA (fR (Ei)) , Ei = {eij}∀j∈N (i). (6)
The advantage of using relative position and velocity instead of global ones as input features is that this explicitly imposes a spatial invariance to the network, given that collision between two particles is invariant to the global positions they are at.

heading: 4 IMPLEMENTATION
text: We adopt the numerical model in SPH to evaluate the physical quantities like particle density and differential operators like gradient. To construct graph representation for particles, we establish edges between particles within the control radius.

heading: 4.1 UPDATE SCHEME
text: In general, given the state (position xn and velocity vn) of the current time step n, we derive the state of next time step n + 1 by passing the state information through advection net, collision net, and pressure net sequentially. The input features to advection net are positions and velocities of particles, [xn,vn], along with g, which indicates the external body force per mass of fluid, and viscosity parameter ν, which denotes the magnitude of fluid viscosity. The advection net predicts acceleration of particles:
aadv = Gadv (x n,vn,g, ν) , (7)
and updates the state of fluid particles to an intermediate state [x∗,v∗].
v∗ = vn + aadv∆t, (8) x∗ = xn + v∗∆t. (9)
Where aadv = [aadv1 , ...,a adv N ], x n = [xn1 , ...,v n N ], v n = [vn1 , ...,v n N ] for particles i ∈ {1, .., N}. We will use the same notation throughout illustration.
The collision net takes relative positions and velocities between particles, [x∗i −x∗j ,v∗i −v∗j ] as input, and predicts correction to the velocity,
∆v = Gcol(x ∗ r ,v ∗ r), (10)
where [x∗r ,v ∗ r ] denotes the relative position and velocity in intermediate state. The velocity is then updated with predicted correction: v∗∗ = v∗ + ∆v. (11)
The updated intermediate position and velocity are taken as input by the pressure net, along with particle number density ρ.
p̂ = Gpres(x ∗,v∗∗,ρ). (12)
The state of fluid field is then updated to next time step n+ 1,
vn+1 = v∗∗ − ∇p̂ ρc ∆t, (13) xn+1 = x∗ + vn+1∆t. (14)
Where p̂ = [p̂1, ..., p̂N ], ρ = [ρ1, ..., ρN ], ρc is the density parameter of fluid. Predicting pressure of fluid field using particle density and velocity is based on the observation that advection will incur a temporary compression on fluid body, which means fluid density has changed. Therefore the goal of pressure net is to impose a pressure projection to mitigate these deviations.
During the above calculation, the global positional information is only used to construct graph on fluid particles and will not be passed into aggregator and processor as features. The relative position and particle density are normalized before input.

heading: 4.2 NETWORK ARCHITECTURES
text: For node-focused graph network, to derive a smooth response of the field with respect to spatial location, the aggregation from layer l − 1 to l is defined as:
a (l) i =
∑ f
(l−1) j W (|ri − rj | , h)∑ W (|ri − rj | , h) + f (l−1) i ,∀j ∈ N (i), (15)
f (l) i = σ ( W · a(l)i ) , (16)
where the aggregator sums up the features f (l−1)j from neighbor vertices {vj |j ∈ N (i)} using smooth kernel as weight function, and here self-connection is added to every vertex. Linear transformation W and non-linear transformation σ are then applied to the aggregated features. In practice we found that two layer of aggregations is enough for the model to produce reasonably accurate output (Adding more aggregation layers does not bring in significant improvements).
As for the edge-focused network, the aggregation is simply defined as: ai = ∑
j∈N (i)
W · rji, (17)
where rji is edge feature processed by processor, W is a linear transformation matrix. The aggregation in the edge-focused network is at the last layer, so no non-linearity is included here.
The processors in both networks are implemented as shared MLP, where they are shared among nodes or edges depending on network types (e.g. in the node-focused network, the processor is an MLP shared among each node). In a node-focused network, the processor has three hidden layers, with the input node embedding f of size 128. In an edge-focused network, the processor has four hidden layers, with input edge attributes [xr,vr] of size 6.

heading: 5 EXPERIMENTS

heading: 5.1 TRAINING
text: Dataset We use the Moving Particle Semi-implicit method (MPS) (Koshizuka and Oka, 1996) with an improved pressure solver (Lee et al., 2011) to generate high-fidelity simulation data of incompressible flow. MPS is a numerical method based on SPH which prioritizes accuracy over calculation speed. It enforces the incompressibility of the fluid field by solving the pressure Poisson equation. We created 20 scenes by randomly placing fluid blocks, solid obstacles, initializing fluid particles with random velocity (See A.2 for full detail of training dataset settings and training strategy).
Loss Function and Optimization We train three networks, advection net, collision net, and pressure net separately. Each network is trained in a supervised way by optimizing the mean squared error between prediction ŷ and ground truth y.
L = 1
N N∑ i=1 ||yi − ŷi||22 (18)
We normalize particle density before inputting into the pressure net, which accelerates and stabilizes training. In the processor, we add LayerNorm (Ba et al., 2016) after activation to each layer (except
for the output layer). The parameters of the model are optimized with Adam (Kingma and Ba, 2014) optimizer. We implement the model in PyTorch. All the training and experiments are mainly carried out on NVIDIA GTX 1080Ti GPU.

heading: 5.2 EVALUATION
text: Baselines Besides the comparison against ground truth data, we also compare our model to two recent works that use data-driven approaches to simulate fluids under the Lagrangian framework. Ummenhofer et al. (2020) use the continuous convolutional kernel to learn fluid dynamics and they reported that their model has outperformed other works in this domain. Sanchez-Gonzalez et al. (2020) propose a graph network-based simulator (GNS) as a general-purpose physic simulator under Lagrangian representation. Our model and GNS both transform and pass messages of fluid field via graph structures, but GNS consists of a far larger and deeper network with multiple sub-blocks and thus contains much more parameters than ours. For all baseline models, we adopt the same training strategy from original papers but train them on our dataset (See A.5 for full detail).
Metrics To conduct quantitative analysis, we evaluate model performance based on several metrics. We report the asymmetric version of Chamfer Distance between the simulated results of different models and ground truth sequence. The asymmetric Chamfer distance for two collections of particles X,Y (from X to Y ) is defined as:
L (X,Y ) = 1
N ∑ x∈X min y∈Y d (x,y) , (19)
where N is the total particle number of point cloud collection X , and distance function d(x,y) is evaluated using L2-norm ‖x−y‖2. We investigate two essential physical quantities for incompressible fluid simulation - velocity divergence and particle density deviation of fluid field. In addition, we use normalized mean absolute error (MAE) and relative tolerance to evaluate the error of advection net and pressure net on single frame inference respectively. The normalized MAE from prediction ŷ to ground truth y is defined as:
LMAE = 1
N ∑ |ŷ − y| |y| . (20)
The relative tolerance of the numerical solution x̂ to a linear system Ax = b can be defined as:
tol = ||Ax̂− b||2 ||b||2 . (21)

heading: 6 RESULTS
text: Performance In order to measure the performance of our FGN model as a physic simulator, we performed simulations on several different test cases. In the first two cases, dam collapse and water
fall, we qualitatively measure the results of different models via visualization1 of fluids in Figure 4. We report the physical property of simulation results and its Chamfer distance to ground truth data. Quantitative results over the whole simulation sequence are listed in Table 1 (See Figure 9 in A.3 for error trend figures). The results show that our model FGN gives the best accuracy in retaining physical properties and position prediction. In addition, we study the performance of each sub-network in our model as stand-alone solvers for sub-dynamical systems and report their relative error in Table 2. For the advection net, we challenge it by applying a different set of material parameters (i.e. different gravity and viscosity parameters). We report the normalized mean absolute error (MAE) between prediction and ground truth. For pressure net, we evaluate the relative tolerance of its predicted solution p̂ to the discretized pressure poisson equation (i.e. Ap = b). The relative low error demonstrates the capability of our sub-networks in learning and predicting physics.
1Video link: https://sites.google.com/view/fluid-graph-network-video/home
Generalization To test out the model’s capability of generalization. We apply our model to test cases with conditions that are beyond training distributions. In the first case, we study how our model will predict the pressure distribution of circular flow around a cylinder. This scene contains an inflow on the left side which keeps emitting particles during the simulation and an outflow on the right side. We challenge our model’s robustness by applying different time step sizes to it. Figure 4 shows qualitative comparison of model output under different time step sizes. In addition, we simulate on two scenes which contain much more complex geometries. Visualizations of two test scenarios containing complex geometries are shown in Figure 5. Although our model is trained with only a fairly small dataset, it remains accurate under several different conditions. This demonstrates it capability of generalization and robustness. (More details on quantitative analysis are in A.3. Figure 10 shows position error trend under different time step sizes. Figure 11 shows position error trend under complex scenes.)
Ablation Study We test the performance of different types of aggregators used in pressure net, as pressure net has the largest impact on overall prediction accuracy. 2 We compare our aggregator against graph convolution networks (GCN) from Kipf and Welling (2016), Hamilton et al. (2017)’s graph SAGE using mean aggregator , and MLP w/o any graph aggregation operation. In general, all aggregators give a similar performance on overall position prediction (similar Chamfer distance error), yet our model significantly improves the quality of predictions in terms of maintaining a constant density (See Table 4 in A.4 for full detail).
In addition, we report the model size (parameter number) and runtime benchmark in Table 3. Our model is very efficient in terms of training and inference, as it has far less trainable parameters than others.

heading: 7 CONCLUSION
text: In this paper, we present a data-driven Lagrangian fluid model for incompressible fluid simulation by decomposing simulation scheme as separate reasoning tasks based on Navier-Stokes equation. It can preserve many essential physical properties of the fluid field such as low volume compression, and predict reasonable pressure distribution. Our model also has generalization capability, where it can remain stable when extrapolating to a wide range of different geometries and adopting different time step sizes. In general, our work is an advance in learning on unstructured data with graph neural networks, and enriches the paradigm of combining learning-based methods with physical models as well.
2During the preliminary experiment, we find that aggregators in advection net do not have significant impact on overall model performance when viscosity parameter of fluids is small, therefore ablation results on it are not listed here; Aggregator in collision net is just summing up all the predicted edge attributes and conducting linear transformation. Using more advanced structure does not bring in further improvement.

heading: APPENDIX A
text: A.1 IMPLEMENTATION DETAILS
Graph construction We build the graph by establishing edges between particles within limit radius. We perform the neighborhood searching on GPU by using cell sort algorithm. All the edge attributes are stored in sparse matrices. Although the limit radius does not have significant impact on training loss and larger limit radius will increase the computing cost drastically, we found that small limit radius can influence the long-term stability of the model. For advection net and pressure projection net, we select the limit radius to be three times the particle diameter D (D = 0.050m), and 0.9D for collision net.
Numerical model In general, we calculate all the gradient operators and smooth kernel based on numerical model of SPH.
In SPH model, particle density is defined as: ni = ∑ i 6=j W (|ri − rj | , h),∀j ∈ N (i). (22)
For scalar quantity φi at location ri, we approximate its gradient by:
∇φi = d
n∗ ∑ j (φj − φi)∇Wij , (23)
where n∗ is the constant particle number density derived by calculating the maximum particle number density at the initial frame, d is the dimension of the problem, Wij denotes the smooth kernel function value between particle i and j. Similarly velocity divergence is defined as:
∇ · vi = d
n∗ ∑ j (vj − vi) · ∇Wij , (24)
We adopt the same smooth kernel function from Koshizuka and Oka (1996), which is very simple to evaluate.
Wij =
{ h
||ri−rj ||2 − 1 if ||ri − rj ||2 < h 0 if ||ri − rj ||2 ≥ h
(25)
A.2 TRAINING
Dataset Generation We place one or two of the following basic obstacles (as shown in Figure 6) in training scenes. In addition to obstacles, each scene is a cubic box (80x80x40) containing a fluid block (25x25x10) (as shown in Figure 7). We place the fluid block at random place in the box and initialize its velocity of one direction by uniformly sampling from U (0, 0.1). We generated 20 scenes adopting above settings and simulate each training scene to 1000 time steps with step size dt = 0.002.
Strategy In each time step of the ground truth simulator, there are mainly three steps. First, advect fluids with body force and viscosity:
v∗ = vn + aadv∆t, (26) x∗ = xn + v∗∆t, (27)
and then solve the pressure poisson equation,
∇2p = 1 ∆t ∇ · v∗, (28)
lastly,
vn+1 = v∗ − ∇p ρc ∆t, (29) xn+1 = x∗ + vn+1∆t. (30)
Here, we use [x∗,v∗,g, ν] as the training features and a∗ (i.e. (v∗ − vn)/∆t) as the training label for advection net. [x∗,v∗] and particle density ρ (evaluated based on x∗) are the training features for pressure net, with pressure p as label. To train collision net, we simulate another particle system that is updated only based on elastic collision rule and applies no other dynamics. We use relative velocity and position before collision as inputs, and use velocity difference (i.e. ∆v) as output target.
We train each network for 100,000 iterations of gradient updates and decay learning rate from 0.001 to 0.0000625. For the training of three sub-networks (advection net, collision net, pressure net), the batch size of each network is 16, 4 and 32 respectively. To allow the mini-batching of different graph, we mini-batch the adjacency matrix of different scenes by creating a large sparse matrix and stacking adjacency matrix on the diagonals.
A.3 METRICS EVALUATION DETAILS
We show the error trend of dam collapse and water fall scene under different evaluation metrics in Figure 9. The position loss analysis of different time step sizes is shown in Figure 10. Position error of complex scenarios is shown in Figure 11.
A.4 ABLATION STUDY DETAILS
Performance of different graph message aggregate structures is listed in Table 4.
A.5 BASELINES IMPLEMENTATION
Continuous Convolution We use the open-source implementation from Ummenhofer et al. (2020) 3. To give a fair benchmark result we train their network with our dataset. Note our training dataset is much smaller than theirs and does not include complex geometries. As in our work, we model
3https://github.com/intel-isl/DeepLagrangianFluids
solid obstacles and wall as virtual particles, so we transform these virtual particles into surface and corresponding normals before inputting them into CConv network. The original CConv uses a time step size of 0.02s, but given such a large time step size, the qualitative results can be distinct from ground truth and other model with smaller time step size. Hence during training and comparison we adopt a time step size of 0.002s for CConv model.
Graph Network-based Simulator We implemented GNS following the description in SanchezGonzalez et al. (2020). We build GNS with 10 unshared GN blocks, conditioned on 5 previous velocities and input relative positions as edge features. We chose the connectivity radius to be 2.1D, so that the number of neighbors is around 20. Sanchez-Gonzalez et al. (2020) use finite difference to calculate acceleration and velocity, but in our implementation we explicitly maintain an array to store the velocities of all particles. Additionally, we do not use learned embedding but simple zero and one to indicate particle material type, as in our testing domain there are only two kinds of particles - solid and fluid. The loss function and training procedure are implemented as described in Sanchez-Gonzalez et al. (2020), including noise injection and similar normalization techniques.
For the implementation of GN block, as Sanchez-Gonzalez et al. (2020) states "We use GNs without global features or global updates (similar to an interaction network)", so we implement the GN block update mechanism following the description in Battaglia et al. (2018).
e′k = φ e(ek,vrk ,vsk ,u) (31) v′i = φ v(ê′k,vi,u) (32) ê′i = ρ e→v (E′i) (33)
where φ is the update function and implemented as MLP here, ρ is aggregation function which aggregates all the edge attributes to its center vertex, u is the global feature and here we append them to the nodal features as input feature, rk denotes receiver vertices and sk denotes sender vertices.
In the testing stage, as authors did not state how to set initial velocities, so we just warm start the simulation by calculating the first 5 frames using MPS method and apply GNS to the rest frames.

Review summaries:
The authors propose a learned model specialized on learning Lagrangian fluid dynamics for incompressible fluids. The model is a hybrid between a simulator with explicit advection, collision and pressure correction stages, and a learned model, trained by supervising each of those stages. The authors demonstrate improved stability/conservation of physical properties for a model, and some flexibility to the time-step being changed at test time.

The paper deals with the prediction of 3D Lagrangian Fluid Simulations. Therefore the problem is divided into 3 subproblems, oriented on numerical simulations. An advection part, where the acceleration of the particles is calculated, a collision step, where the boundary effects are included, and a pressure prediction part, where the pressure for maintaining the volume is determined. A graph-based network is used for each part, which is either node or edge-based according to the requirements.