Title: LAGRANGIAN FLUID SIMULATION WITH CONTINUOUS CONVOLUTIONS

Abstract: We present an approach to Lagrangian fluid simulation with a new type of convolutional network. Our networks process sets of moving particles, which describe fluids in space and time. Unlike previous approaches, we do not build an explicit graph structure to connect the particles but use spatial convolutions as the main differentiable operation that relates particles to their neighbors. To this end we present a simple, novel, and effective extension of N-D convolutions to the continuous domain. We show that our network architecture can simulate different materials, generalizes to arbitrary collision geometries, and can be used for inverse problems. In addition, we demonstrate that our continuous convolutions outperform prior formulations in terms of accuracy and speed.

1. INTRODUCTION:

Understanding physics can help reasoning about our environment and interacting with it. Neural networks have emerged as a particularly promising approach to capture the complexity of natural phenomena from data (Ling et al., 2016; Tompson et al., 2017; Morton et al., 2018 ). An important aspect of learning physics with neural networks is the choice of representation. Lagrangian representations based on particles are particularly popular and have supported recent results with rigid bodies, deformable solids, and fluids (Battaglia et al., 2016; Mrowca et al., 2018; Li et al., 2019) . Many of these approaches use graph structures to define interactions; the existence of an edge determines in a binary fashion whether two particles interact.

However, a wide range of physical processes such as fluid mechanics are described by continuous partial differential equations rather than discrete graph structures. The continuous, volumetric, and tightly coupled nature of these processes causes inherent difficulties for graph-based approaches, such as a large number of edge connections that must be established, tracked, and disengaged as the particles move. In this work, instead of using graphs as the underlying representation, we adopt a continuous viewpoint. We propose to use convolutional networks (ConvNets) with continuous convolutions on particles for learning fluid mechanics. We treat fluids as spatially continuous functions sampled at a finite set of (continuously evolving) positions and process them with a novel continuous convolution layer. This matches the continuous nature of the problem more closely and simplifies the definition of neural networks by abstracting the underlying particle representation.

We extend the grid-based filter representation commonly used for discrete convolutions to the continuous domain by simple linear interpolation. Linear interpolation of the filters allows efficient lookup of spatially varying filter values at arbitrary positions while retaining the compactness and efficiency of the grid representation. In addition, we use a window function to define the support of the filters and a ball-to-cube mapping to support spherical receptive fields. We show that our convolutions, despite their simplicity, perform better than more sophisticated representations (Wang et al., 2018; Schenck & Fox, 2018) .

With the presented continuous convolution layer, we develop an efficient ConvNet architecture for learning fluid mechanics. The network processes sets of particles. We use dynamic particles to represent the fluid and static particles to describe the boundary of the scene. Modeling the scene boundary with particles makes it easy to apply our network to new scenes and allows the network to learn collision handling in a unified framework. Our network generalizes to arbitrary obstacle configurations and can simulate a range of material behavior. To demonstrate the usefulness of a learned -and hence differentiable -fluid simulator, we show that material properties can be estimated from observed simulation data. Experimental results indicate that the presented approach outperforms a state-of-the-art graph-based framework (Li et al., 2019) .

2. RELATED WORK:

Fluids encompass a range of materials that are important in everyday life and throughout science and engineering (Wilcox, 2006) . A prominent set of methods, known as Smoothed Particle Hydrodynamics (SPH), employs a Lagrangian viewpoint to simulate these phenomena. SPH originated in particle-based models for astrophysics (Gingold & Monaghan, 1977) and has become extremely popular for simulating complex flows in many fields of science (Monaghan, 1988) . Among others, it has been highly successful for modeling interface flows (Colagrossi & Landrini, 2003) , complex multi-phase phenomena (Hu & Adams, 2006) , and even magnetohydrodynamics (Price, 2012).

SPH and its variants have been widely used to model complex real-world phenomena for visual effects. Following the earlier development of physics-based fluid simulation for special effects by Foster & Metaxas (1996) , the introduction of SPH (Müller et al., 2003) has led to a large class of powerful algorithms (Solenthaler & Pajarola, 2009; Bender & Koschier, 2015) . These applications often involve complex geometries at large scales, and extensions such as FleX (Macklin et al., 2014) and pressure-aware rigid-body coupling (Gissler et al., 2018) broaden the framework to encompass many physical phenomena.

Lagrangian flows have also been considered in machine learning. The pioneering work of Ladický et al. (2015) demonstrated that flow representations can be learned with regression forests. CNNs were used by Tompson et al. (2017) to accelerate the expensive pressure correction step of gridbased solvers, while other works have focused on super-resolution (Xie et al., 2018) , learning time evolution via the Koopman operator (Morton et al., 2018) , and learning reduced representations (Wiewel et al., 2019; Kim et al., 2019) . Differentiable SPH solvers were proposed to solve control tasks for robotic applications (Schenck & Fox, 2018) . Generic physics simulations for Lagrangian rigid and deformable bodies were considered in a series of works that developed graph-based representations (Battaglia et al., 2016; Sanchez-Gonzalez et al., 2018) . Such graph-based representations were recently applied directly to fluid simulation (Mrowca et al., 2018; Li et al., 2019) . We share with these recent works the goal of modeling Lagrangian fluids with differentiable neural networks, but take a different tack: rather than using graph-based representations, we work with point clouds and continuous convolutions over the spatial domain.

From a technical perspective, our approach is also related to existing works that apply convolutions on point clouds. A number of methods transferred the convolution concept to point clouds in the context of semantic classification and segmentation of 3D objects (Hua et al., 2018; Atzmon et al., 2018; Hermosilla et al., 2018; Li et al., 2018; Su et al., 2018; Wu et al., 2019; Xu et al., 2018; Lei et al., 2019) . Particularly notable in our context is the work of Wang et al. (2018) , who used continuous convolutions to compute the scene flow between two point clouds, and the aforementioned work of Schenck & Fox (2018) , who used convolutions to implement a differentiable version of position-based fluids (Macklin & Müller, 2013) . Both works define continuous convolution operators that can support regression tasks and we compare to them directly in Section 6. Most closely related to our filter representation is the work of Fey et al. (2018) . They use B-splines to define continuous filters and propose to use spherical coordinates to implement spherical receptive fields. Spherical coordinates are problematic due to singularities and require special treatment, which we avoid with a ball-to-cube mapping. Furthermore, the output of the operator of Fey et al. (2018) can be discontinuous. We show in our ablation study in Section 6 that applying a window function to guarantee a continuous output is advantageous for our task.

3. BASICS:

Fluids have been studied for centuries, and the Navier-Stokes equations for incompressible fluids are well established (Batchelor, 1967) :

EQUATION

A common approach to solve these partial differential equations is to approximate the fluid with a set of smooth particles (Monaghan, 1988) . Each particle corresponds to a continuous blob of matter and carries the local properties of the fluid, such as velocity and density, which move with the flow. This is motivated by the fact that a function A(x) can be represented by an integral interpolation

EQUATION

where δ(x) denotes the Dirac delta function. This equation can be discretized as

EQUATION

where V i is the volume at the given point in space and

EQUATION

Here W (x, h) is a smooth kernel or convolution with radius h, usually in the form of a Gaussian distribution, but more complex functions can also be used. In practice, the kernel is finite and yields localized neighborhoods of particles that interact via interactions weighted by the kernel of its derivatives. In this way, the continuous description for fluids from Equation 1 can be discretized in a Lagrangian fashion and solved numerically. Typically, internal forces are calculated based on the local pressure, viscosity, and surface tension, which give an update for the position of each particle. Below, we adopt the position-based fluids (PBF) method (Macklin & Müller, 2013; Macklin et al., 2014) , which likewise is based on SPH, but reformulates the updates as constraints on the positions.

4. CONTINUOUS CONVOLUTIONS:

The discrete convolution operator as commonly used in ConvNets is defined as

EQUATION

where f and g are the input and the filter function, x is the position, τ is the shift vector, and Ω is the set of shift vectors that defines the support of the filter function. On regular data such as images, the positions x range over a regular grid and the shift vectors τ are integer-valued, i.e. x, τ ∈ Z d for some dimensionality d. Analogously in the continuous domain, this convolution is defined as

EQUATION

where x and τ are real-valued vectors, i.e. x, τ ∈ R d .

We now adapt this definition to unstructured point clouds. In this setting we have a finite number of points that sample the function f but do not lie on a grid. For a point cloud with i = 1, .., N points with values f i at positions x i , we define the convolution at position x as

EQUATION

N (x, R) is the set of points within a radius R around x. a is a scalar function that can be used for density normalization specific to the points x i and x as in Hermosilla et al. (2018) . In the simplest case, a can be constant: a = 1. In our case we want to ensure a smooth response of our convolution under varying particle neighborhoods, therefore we define a as a window function:

EQUATION

A similar function has been used by Müller et al. (2003) in the SPH framework. ψ is another scalar function for normalization, which can be set in our implementation as either We use ψ(x) = 1, since changes in the density of particles are an important feature for simulating fluids.

EQUATION

For the filter function g we simply use a regular grid to store the filter values but use linear interpolation to make g a continuous function. In addition, we use a mapping Λ(r) of a unit ball to a unit cube to implement spherical filters as shown in Figure 1 . We use the mapping described by Griepentrog et al. ( 2008) and give the detailed function in the appendix. The intermediate coordinate mapping Λ provides the flexibility to implement different spatial shapes while keeping the advantages of a regular grid for the storage and lookup of filter values.

Note that Equation 7 uses a similar approximation as in the SPH framework (Equation 3). Assuming that each point represents the same volume, V i is a constant factor in Equation 3, which we drop in our definition.

5. LEARNING FLUID MECHANICS WITH CONVOLUTIONAL NETWORKS:

Our goal is to learn fluid mechanics from observing the motion of particles. The input to our ConvNet is a set of particles with corresponding features. Since position itself is not a feature but simply defines the particle's position in space, we must assign a feature vector to each particle.

The feature vector we use is a constant scalar 1 accompanied by the velocity v and the viscosity ν.

A particle p n i at timestep n with its position and input feature vector is thus a tuple (x n i , [1, v n i , ν i ]). Defining the velocity explicitly as an input feature allows us to compute intermediate velocities and positions as in Ladický et al. (2015) and to apply external forces and pass this information to the network. We compute the intermediate positions x n * i and velocities v n * i beginning with timestep n with Heun's method as

EQUATION

)

The vector a ext is an acceleration through which we can apply external forces to control the fluid or to simply apply gravity. The intermediate positions and velocities lack any interactions between particles or the scene, which we are going to implement with a ConvNet. To enable the network to handle collisions with the environment we define a second set of static particles s j . We sample particles on the boundaries of the scene with normals n j as the feature vectors, i.e. s j = (x j , [n j ]).

Our network implements the function

EQUATION

) which uses convolutions to combine features from both particle sets. ∆x is a correction of the position which accounts for all particle interactions including the collision handling with the scene. Finally, we apply the correction to update positions and velocities for n + 1 as

EQUATION

)

Note that the updated position x n+1 i depends on the output vector ∆x i and allows us to directly define our learning objective on the particle positions.

Figure 2 : Schematic of our network with a depth of four. In the first depth level we compute convolutions at each dynamic particle location with the set of static particles that defines the environment as well as the dynamic particle set. We also directly process the features of each particle via a fully-connected stream. In the following levels, we compute convolutions only on the dynamic particles. At each level we use addition to aggregate the features computed by convolutions and fully-connected layers. Between the second and third level we also include a residual connection. The final level generates the position correction ∆x. Operations annotated with a * are followed by the ReLU activation function. All CConv and FC operations use an additive bias.

5.1. NETWORK ARCHITECTURE:

We use a simple convolutional architecture with an effective depth of four. An overview of the network is shown in Figure 2 . Since we want to compute the correction for all dynamic particles in our scene, we compute convolutions for the intermediate positions defined in Equation 11. Our network is a sequence of continuous convolutions (CConv), which are defined by an input particle set, the positions at which we want to evaluate the convolution, its filters G, and the radius R. For instance, to describe a convolution on the static particles s i at intermediate positions

EQUATION

where f i are the computed output features for each position x n * i . G is a 5D array storing all filters in the layout [width, height, depth, ch in , ch out ]. In contrast to discrete convolutions on a regular grid, the spatial filter dimensions here do not define the receptive field but the resolution of the filters. The receptive field depends only on the radius R, which specifies the spatial extent. Throughout our network we use filters with a spatial resolution of [4, 4, 4] and a radius of 4.5 times the particle radius.

For convolutions within the dynamic particles we exclude the particle at which we evaluate the convolution and instead process the particle's own features in a stream of fully-connected layers. After each depth level we then combine the result from the convolutions and the fully-connected layers by addition. This can be interpreted as a convolution with a spatial resolution of 4×4×4 + 1. We found that this design improves accuracy and allows us to use smaller filters with even sizes (see Table 2 ).

5.2. TRAINING PROCEDURE:

We train our fluid simulation network in supervised fashion based on particle trajectories produced by classic ("ground-truth") physics simulation. Our loss is defined as follows:

EQUATION

)

The ground-truth position at timestep n + 1 is denoted by xn+1 (Bender & Koschier, 2015) , a high-fidelity solver that was executed with small timesteps (down to 0.001s). Despite using a much larger timestep (0.02s), our convolutional network produces results of comparable visual fidelity. Note that our particles are falling slightly more slowly due to differences in the integration of positions and the much larger timestep. See the supplementary video. neighbors across our experiments. Particles with few neighbors are close to the surface or interact with the scene boundary. Both cases are important for fluid simulation because particles near the surface define the liquid-air interface, which is particularly salient, and particles near the scene boundary require collision handling. The parameter γ = 0.5 makes our loss function more sensitive to small particle motions, which is important for increasing the accuracy and visual fidelity for small fluid flows. During training we predict particle positions for two future timesteps, namely n + 1 and n + 2. The combined loss L is

EQUATION

We found that optimizing a loss defined over two frames improves the overall quality of the simulation. (Optimization for three frames did not result in further improvements.) We optimize L over 50,000 iterations with Adam (Kingma & Ba, 2015) and a learning rate decay with multiple steps, starting with a learning rate of 0.001 and stopping with 1.56 • 10 -5 .

5.3. DATASETS:

We have trained our network on multiple datasets. For quantitative comparisons with prior work we trained our network on the dam break data from Li et al. (2019) . The scene simulates the behavior of a randomly placed fluid block in a static box. We generate 2000 scenes for training and 300 for testing. The data was generated with FleX, which is a position-based simulator that targets real-time applications (Macklin et al., 2014) .

We also trained our network on more challenging data generated with DFSPH (Bender & Koschier, 2015) , which prioritizes simulation fidelity over runtime. DFSPH can generate accurate fluid flows with very low volume compression: a desired property. We generate ground-truth data by randomly placing multiple bodies of fluid in 10 different box-like scenes and simulating them for 16 seconds each with an adaptive timestep of up to 1 kHz. The time resolution of the generated data is 50 Hz. We show a qualitative comparison of our method to the ground truth in Figure 3 . We generate 200 scenes for training and 20 scenes for the test set. To train networks that can deal with multiple materials, we additionally generate 200 scenes with fluids of varying viscosity. For estimating material properties, we generate 7 test scenes that only differ in the viscosity parameter.

6. EVALUATION:

Baselines. We compare our method to DPI-Nets (Li et al., 2019) , which were previously shown to significantly outperform prior formulations such as hierarchical relation networks (Mrowca et al., 2018) . Figure 4 provides a qualitative comparison. Quantitative results are reported in Table 1 . To analyze the accuracy of the forward step for each method, we compute the average error of the particle positions with respect to the ground truth. We use every 5 th frame for initialization and compute the deviation from the ground truth for two subsequent frames, denoted by n + 1 and n + 2. In addition, we report the average distance from the ground-truth particles to the closest particle in the prediction for the whole sequence, to measure long-term similarity. We compute the distance for frame n as

EQUATION

where X n is the set of predicted particle positions for frame n, xn i is the ground-truth position for particle i, and N is the number of particles.

We also compare our continuous convolution formulation to continuous convolution representations used in SPNets (Schenck & Fox, 2018) , PCNN (Wang et al., 2018) , KPConv (Thomas et al., 2019) , and SplineCNN (Fey et al., 2018) . To this end, we plug the respective convolution operators into our network architecture as shown in Figure 2 . This facilitates controlled comparisons in which the overall architecture and simulation setup are fixed and only the convolution operators are varied. As shown in Table 1 , our method outperforms all baselines with respect to both accuracy and inference time.

We have trained and tested all methods on the dam break dataset from Li et al. (2019) as well as our data (generated with a high-fidelity simulator (Bender & Koschier, 2015) ). To facilitate the training on our data for all methods, we generated a simplified version with a constant number of particles (6,000) and a single box environment (as shown in Figure 4 ). We train DPI-Nets for 5 epochs and all other networks for 50,000 iterations, which corresponds to 2.7 epochs for the dam break data and 5 epochs for our data. Training takes about a day for our method with our convolutions on an NVIDIA RTX 2080Ti. Training with PCNN convolutions (Wang et al., 2018) takes about 2 days on 4 GPUs. Note that our reimplementation of PCNN convolutions uses Tensorflow's built-in functions, which consume a lot of memory and necessitate multi-GPU training. The DPI-Nets model trains in about one day. Training with SplineCNN Convs takes 3 to 4 days on a single GPU. We got the best results for this method with spherical kernel coordinates and closed splines. For KPConvs we used a Quadro RTX 6000 with 24 GB of RAM due to the higher memory requirements. Training took about 1 day with 15 kernel points. For SPNets convolutions, we estimated a training time of more than 29 days with 3×3×3 filters by extrapolating from timing of a smaller number of iterations. We thus only report inference time for this method. Note that the convolutions of SPNets were designed to implement the position-based fluids algorithm, while we use them here in a more general network architecture with a much larger number of channels, which explains the very long runtime.

Ablation study. We perform an ablation study to evaluate our decisions in the design of the continuous convolution operator and the network. We study the importance of the interpolation, the window Table 1 : Accuracy and runtime analysis. We compare the average error between the ground-truth particle positions and two predicted future frames on the test set. Additionally, we report the average distance from the ground truth to the prediction over the whole sequence. In this test mode some methods become unstable after a few frames. In the last column we report the average inference time per frame.

function, the loss design, and the architecture choice. The results are reported in Table 2 . This table reports error measures (averaged over the test sequences) that were used in Table 1 , and also reports the errors for two predicted frames initialized with the frames at the end of each sequence to measure the errors for small flow velocities. Large errors for small velocities can yield perceptually salient artifacts: rather than being still, fluid particles jitter or churn. Table 2 : Ablation study. We compare the average error between the ground-truth particle positions and two predicted future frames on the test set, evaluated over whole test sequences (left) and just on the final frames of each test sequence (middle; this focuses on frames with small motion). The rightmost column shows the average distance from the ground truth to the predicted point set over whole test sequences. Ours w/o interpolation uses nearest-neighbor instead of trilinear interpolation for the convolution filters. Ours w/o window uses a(x i , x) = 1 in Equation 7. Ours w/ naïve loss uses Euclidean distance as the loss, i.e. we set γ = 1 and φ i = 1 in Equation 16. Ours w/o FC uses only convolutions and includes the central particle in the convolution (rather than separately processing its features via an FC layer).

Generalization. In Figure 5 we show that our network generalizes well to scenes with drastically different geometry than seen during training. (The training set uses only box-like containers. See the appendix for visualization.) These scenarios demonstrate that we can emit particles during simulation, which can be costly for methods that build and maintain explicit graph structures. We compare generalization performance quantitatively to DFPSH on a complex scene in Figure 6 . Figure 7 demonstrates generalization along a different dimension. Here we show that we can set the viscosity of the fluid at test time to a value not seen during training. The fluid shape used in this example was also not seen during training. Material estimation. In this experiment we apply our network to an inverse problem: material estimation from observation. We use our network to estimate the viscosity parameter of a fluid from the particle movement. 

7. CONCLUSION:

We have developed continuous convolutional networks for Lagrangian fluid simulation. We have introduced a simple formulation for continuous convolutions and demonstrated its accuracy and speed. Our model captures a wide range of complex material behavior, offers long-term stability, and generalizes to new situations such as varying particle counts, domain geometries, and material properties. There are numerous directions for future work, such as extending the framework to incorporate rigid and deformable solids. We will release the code to facilitate such development.

Our continuous convolution implementation will be made available as part of Open3D (Zhou et al., 2018) .

A APPENDIX A.1 IMPLEMENTATION DETAILS:

To accelerate the computation of Equation 7 we use existing general matrix multiplication primitives. Similar to standard convolutions in deep learning frameworks we build a matrix with patches that we then multiply with the filter matrix. We can account for the linear interpolation by applying the interpolation to the patch matrix instead of the filters. For 3D point clouds a single point contributes to up to 8 voxels in a patch.

Another crucial part of our implementation is the nearest neighbor search. Since the positions of the fluid particles update with each timestep, we have to rebuild the neighborhood information for every frame. We implement the neighborhood search with spatial hashing (we use the hash function proposed in Teschner et al. (2003) ). We explicitly store all particle neighbors in a compact list, which allows us to reuse the information for multiple convolutions operating on the same point sets. Table 4 compares the average frame runtimes of our method with the baselines.

A.2 DATASET GENERATION:

To enable generalization to new environments we use 10 different containers (see Figure 8 ) in our data generation process. During scene generation we randomly sample an environment and place up to 3 fluid bodies with different shapes and sizes (see Figure 9 ) and random initial velocities in the scene. We simulate each generated scene with DFSPH using the SPlisHSPlasH framework1 for 4 : Runtime analysis. We compare the average per frame inference time and the time used for the nearest neighbor search (NNS). Our convolution achieves the shortest inference times in comparison even if NNS times would be excluded. Irrespective of that, the methods used for finding neighbors can have a significant contribution to the total runtime. Since fluid particles are moving, acceleration structures for the neighbor search have to be rebuilt each frame. This affects the KD-Tree methods as well as the methods using spatial hashing. Note that we use the same NNS for PCNN and Ours (denoted with *). For all methods except for SPNets the inference time is shorter on the smaller DPI DamBreak data (3456 particles compared to the 6000 particles of our data). We attribute this to a higher number of neighbors on the DPI DamBreak, which is about 49 on average compared to the average 40 neighbors on our datasets. All runtimes were measured on a system with an Intel Core i9-7960 and an NVIDIA RTX 2080Ti.

16 seconds to ensure that each scene contains frames with small particle velocities. To create the particle representation of the box-like containers, we do Poisson-Disc sampling on the mesh surface with the tools provided by DFPSH and add surface normals to each particle. Müller et al. (2003) . Ours triangular window uses a triangular window function.

Ours w/o window does not use a window function. This case can also be interpreted as a rectangular window since we only consider points within a radius R.

A.3 TRAINING DETAILS:

We use the Tensorflow framework for implementing the training procedure. We use Adam as optimizer and train with a batch size of 16 and an initial learning rate of 0.001. We half the learning rate at steps 20000, 25000, . . . , 45000. For the convolutions we use the random uniform initializer with range [-0.05, 0.05]. All other weights are initialized with the respective default initializers of Tensorflow version 1.12. The output of the network is scaled with 1 128 to roughly adjust the output range to the ground truth position correction of the training data.

The unit of length of the training data is meter. The particle radius used in DFSPH is h = 0.025 m. For our convolutions, we use spherical filters with an empirically determined radius of R = 4.5h.

A.4 WINDOW FUNCTION:

We compare 3 different choices for the window function a from Equation 8 in Table 5 . The results show that enforcing a continuous output with a window function yields better results. Further, using a window function similar to kernels used in SPH codes gives better results than a simple triangular window. Learning the window function is therefore a possible direction to extend our framework to further improve the accuracy.

A.5 COORDINATE MAPPING FUNCTION:

We use the ball to cube mapping described in (Griepentrog et al., 2008) to map a position within a spherical region to the filter values stored in a regular grid. We give here the functions used for the 3-D case as used in our implementation. For more details see (Griepentrog et al., 2008) .

:

https://github.com/InteractiveComputerGraphics/SPlisHSPlasH

Figures and Tables:

Figure 1: We use spherical filter shapes for our continuous convolutions but use regular grids to store the filter values. The left part of the figure shows a spherical region with radius R and a point with relative coordinates r with respect to the center. We transform r via a mapping Λ to normalized coordinates in a regular grid. The thin dotted lines illustrate the distortion of the mapping. To look up the final filter value we use trilinear interpolation in the regular grid.

Figure 3: Comparison to ground-truth physics simulation. Two fluid bodies collide. Top: simulation by our trained network. Bottom: simulation of the same scenario by DFSPH(Bender & Koschier, 2015), a high-fidelity solver that was executed with small timesteps (down to 0.001s). Despite using a much larger timestep (0.02s), our convolutional network produces results of comparable visual fidelity. Note that our particles are falling slightly more slowly due to differences in the integration of positions and the much larger timestep. See the supplementary video.

Figure 4: Qualitative comparison with DPI-Nets on a test sequence from our dataset. Two fluid bodies collide inside a container. DPI-Nets works well on data with little variance but has problems with more complex scenes and high particle velocities. The DPI-Nets simulation becomes unstable immediately after the fluid hits the box. The fluid behavior predicted by our network matches the ground truth more closely and remains stable for the whole sequence. The two networks have been trained on the same data. Test sequences are distinct from training sequences. See the supplementary video.

Figure 5: Generalization to environments with drastically different geometry than seen during training. Top: we use an emitter to fill up a virtual river with fluid particles, demonstrating generalization with respect to scene geometry and the number of particles. Bottom: a waterfall scene showing the fluid particles and the particle representation of the environment. See the supplementary video.

Figure 8: We sample from 10 different box-like containers during data generation. For the simplified version of our dataset used in the quantitative comparison with the baselines we only use the first container (leftmost container in the first row).

Figure 9: We randomly place fluid bodies of different initial shapes in the scene during data generation. We sample from 5 different shapes and vary the size, the orientation and the initial particle velocity. All particles from the same fluid body start with the same initial velocity. The image shows the particles generated from each shape for a specific size and orientation.

The training data for this experiment contains 200 sequences with random viscosity parameters between 0.01 and 0.3. For testing we generate 5 new sequences with 100 frames and random viscosity within the same range as used during training and use an initial fluid shape not present in the training data. To test generalization we generate two sequences with viscosity values outside the training range, namely 0.35 and 0.4. To estimate the viscosity we backpropagate through the trained network and optimize ν with gradient descent. Table3reports the results, which indicate that our network can be used to estimate material properties from observed data.

Table content: <table><tr><td>GT viscosity</td><td>0.044</td><td>0.127</td><td>0.174</td><td>0.233</td><td>0.269</td><td>0.350</td><td>0.400</td><td>Mean</td></tr><tr><td>Avg. estimated viscosity</td><td>0.027</td><td>0.150</td><td>0.202</td><td>0.255</td><td>0.277</td><td>0.322</td><td>0.336</td><td/></tr><tr><td>Avg. relative error (%)</td><td colspan="3">38.377 18.542 15.604</td><td>9.568</td><td>3.235</td><td colspan="2">7.954 16.016</td><td>15.614</td></tr></table>

Application to an inverse problem: material estimation from observed fluid motion. To estimate the viscosity we backpropagate through our network and optimize ν with gradient descent. For each scene, we run the procedure 10 times, each time with random initialization, and report the average. Viscosity values 0.350 and 0.400 are outside the range that was used during training and are used to test generalization.

Table content: <table/>



Table content: <table/>

Comparison of different window functions. Ours uses a window function similar to the poly6 kernel used in

Table content: <table/>

Peer review (iclr2020):
This paper applies 3D convolutions to the problem of Lagrangian fluid simulation. The primary difficulty in this is that, unlike for Eulerian fluid simulation, which represents the fluid as a grid and adapts nicely to 3D convolutions, in Lagrangian simulations the fluid is represented as an unordered set of particles. It is not straight-forward to apply 3D convolutions on such a data structure, however this paper proposes a method to apply the same regular-grid kernels used in grid-based convolutions to the particle structure. To do this, several points around the kernel are evaluated by first using trilinear interpolation between the particles to get feature values at those points and then convolving those values with the kernel weights. This results in a new particle in the next layer up with those features. In the paper, this method is used to train the weights of the network to reproduce fluid dynamics generated by a simulator. The results show that the proposed method was able to model fluid dynamics over 2 timesteps more accurately than other methods and can do so quickly. While I have some reservations about this paper (detailed below), on the whole I think it is a quality contribution and should be accepted. This paper contributes a novel method for performing 3D convolutions on unordered particle sets, and it shows that the learned fluid dynamics generalize to novel situations. One major hurdle to applying modern convolutional learning techniques to Lagrangian methods is the mismatch between the layout of the data (unordered particles) and the layout of the kernels (regular grid). This paper presents a novel way of bridging that divide, and it shows that the proposed method actually works by applying it to the problem of fluid dynamics and successfully learning it. However, one major concern I had was that it seems all of the training data was generated in box-like environments, which could easily lead to overfitting. This was alleviated by the results showing that although the network was trained only in boxes, it generalized to environments with channels and waterfalls (as seen in the video). This is a powerful result and shows that this method really did learn fluid dynamics and not just a shortcut that only works in boxes. I do think this paper can be improved in a few aspects however. The biggest issue is that the quantitative analysis of the core functionality (reproducing fluid physics) is lacking. The paper only reports results for error after at most 2 timesteps, which is not nearly long enough to determine if the output is accurate. Furthermore, the results are only reported for the box scenes, not the generalization scenes mentioned above. Qualitatively, from the videos, it is clear that the output does at least somewhat model fluid dynamics, but it would be much better to have hard numbers to back that up. I suspect the authors discovered that Lagrangian systems are sufficiently chaotic that after only a few timesteps the particle positions have diverged significantly. This is not a bug but a feature of such systems. In Lagrangian fluids, the particles are but an approximation of the fluid, and unlike Eulerian systems*, multiple different sets of particles can approximate the same fluid. This makes particle position only useful as a measure of error if the trained model can perfectly reproduce the fluid dynamics. But of course it can*t (trained networks aren*t ever perfect in practice), and so small errors quickly compound into large particle position disparities. So even though the trained network models the fluid well overall, the particles end up in completely different locations. Instead a better error metric would be something like measuring the difference between the surface of the fluids, or the velocities or densities at various locations. These are agnostic to the particular particle positions, but still measure how well two different sets of particles represent the same fluid. Using a metric like this, it would be nice to see error graphs over time for both the box and generalization scenes. A couple other smaller points. The chaotic divergence behavior of DPI-Nets seems inconsistent with that paper. Is this possibly a bug in the way it was implemented here? Additionally, the paper states that the convolutions of SPNets were *specifically designed to implement the position-based fluids algorithm* but that it was used in the paper with a much larger number of channels. If it was designed only to work for that one algorithm, how were the number of channels increased? That is unclear. Also, the average error for SPNets is not shown in Table 1 and it is not stated why. *Assuming same grid shape, size, and position.