Symmetric Basis Convolutions for Learning Lagrangian Fluid Mechanics
Authors: Rene Winchenbach
Nils Thuerey
Venue: ICLR 2024 poster (2024)
Keywords: CConv
GNN
SPH
Lagrangian Fluids
Fourier Methods
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)
TL;DR: Using Fourier Terms we build a Continuous Convolution Network that outperforms prior work by leveraging even and odd symmetries.

Abstract: Learning physical simulations has been an essential and central aspect of many recent research efforts in machine learning, particularly for Navier-Stokes-based fluid mechanics. Classic numerical solvers have traditionally been computationally expensive and challenging to use in inverse problems, whereas Neural solvers aim to address both concerns through machine learning. We propose a general formulation for continuous convolutions using separable basis functions as a superset of existing methods and evaluate a large set of basis functions in the context of (a) a compressible 1D SPH simulation, (b) a weakly compressible 2D SPH simulation, and (c) an incompressible 2D SPH Simulation. We demonstrate that even and odd symmetries included in the basis functions are key aspects of stability and accuracy.
Our broad evaluation shows that Fourier-based continuous convolutions outperform all other architectures regarding accuracy and generalization. Finally, using these Fourier-based networks, we show that prior inductive biases, such as window functions, are no longer necessary. An implementation of our approach, as well as complete datasets and solver implementations, is available at https://github.com/orgs/tum-pbs/SFBC.

Published as a conference paper at ICLR 2024
SYMMETRIC BASIS CONVOLUTIONS FOR LEARNING LAGRANGIAN FLUID MECHANICS
Rene Winchenbach & Nils Thuerey
Physics-based Simulation Group
Technical University Munich
(rene. winchenbach, ni 1s. thuerey)@tum.de
ABSTRACT
Learning physical simulations has been essential and central aspect of many recent research efforts in machine learning, particularly for Navier-Stokes. based fluid mechanies. Classic numerical solvers have traditionally been computation- ally expensive and challenging to use in inverse problems, whereas Neural solvers aim to address both concerns through machine learning. We propose general formulation for continuous convolutions using separable basis functions as a su- perset of existing methods and evaluate large set of basis functions in the context of (a) compressible 1D SPH simulation, (b) weakly compressible 2D SPH sim- ulation, and (c) an incompressible 2D SPH Simulation. We demonstrate that even and odd symmetries included in the basis functions are key aspects of stability and accuracy. Our broad evaluation shows that Fourier-based continuous convo- lutions outperform all other architectures regarding accuracy and generalization. Finally, using these Fourier-basee networks, we show that prior inductive biases, such window functions, approach, as well as complete datasets and solver implementations, is available at https:: / wwg.titthbh... com/ ppcltpppps/ae pbs/SFBC SFBC
longer necessary. An implementation of our
INTRODUCTION
Physical simulations play an essential role in many fields of science, e.g., Computational Fluid Dynamics (CFD) (Zhang et al. 2021). with long history of classical numerical research aimed at improving the performance and accuracy of solvers (Vacondio et al., 2021). Nonetheless, numerical solvers typically rely on many insights and intuitions, e.g., which solver and time-stepping scheme to use and which conservation laws to enforce explicitly and implicitly (Popel 2000). While neural
solvers traditionally focus on learning from data (Ummenhofer et al. 2019 Morton al. 2018)
there is often a clear benefit to the integration of inductive biases to enforce dilliculltto-leam aspects,
e.g., rotational invariance (Satorras al. 2021 Lino al. 2022, Thomas et al. 2018), momentum
conservation (Prantl et al. 2022 and symmetry (Wang et al. 2021).
Within CFD, the underlying Partial Differential Equations (PDEs) are either solved on grids with Fulerian approaches or using particles with Lagrangian methods. For Fulerian approaches on struc- tured grids, Convolutional Neural Networks (CNNs) have found broad adoption (Tompson et al.]
2017 Bar-Sinai et al., 2019 Thuerey et al.) 2020); however, applying CNNs to unstructured grids is not directly possible. Graph Neural Networks (GNNs) (Pfaff et al., 2021 Zhou et al., 2020 approaches have been popular for unstructured grids, e.g., GNS (Sanchez-Gonzate et al. 2020).
where cell centers are treated as vertices and cell faces as graph edges. These GNNs have found broad application for various PDEs, e.g., Burger's equation, for unstructured grids and particl-- particle-based
systems [Brandstetter et al., 2022b Allen et al.) 2022). Continuous Convolutions (CConvs) (Wang al. 2018) a subset of GNNs, where interactions
based solely the relative position of cells or particles. These continuous approaches are especially suited for Lagrangian simulations, as particle positions change continuously over time. CConvs often use convolutional filters built using interpolation approaches instead of using Multilayer-Perceptron inherently useful properties, such as smoothness (Fey et al. 2018) and antisymmetry (Prantl al. 2022). Consequently, CConv approaches learn more accurate solutions than general GNNs for an equivalent number parameters in physical problems (Prantl et al. 2022).
(MLPss, as these filters have

---

Published as a conference paper at ICLR 2024
In this context, we propose (a) generalized for- mulation of CConvs using separable basis functions, (b) a Fourier-based architecture with even and odd symmetry for improved inference accuracy and (c) novel dataset consisting of four test cases designed to be challenging while exhibiting quantifiable phys- ical behavior. We perform an exhaustive set of ab- lation studies to investigate which classes of basis functions and inductive biases are beneficial for accurately learning Lagrangian flow physics.
caep for pamile
cordi coreo d01
ataiiiiie Lokk
Figure I: Visual overview of SFBC.
BACKGROUND AND RELATED WORK
followed by an overview of clas- sical solvers in Section 2.2 and prior work regarding machine learning for simulations in Section 2.3
We will first outline our definition of neural solvers in Section 2.1
2.1 NEURAL NETWORK TIME INTEGRATION
At a given time i with particles at position x' E R4, with d spatial dimensions, and corresponding velocities v² E Rd, the goal of physics simulation is to compute a new set of positions and veloc- ities at a new time-point (", i.e., x"" and v"", subject to PDE ō'u F(t, x, u, d² n..... where u is being solved for. Based on the PDE, solver now computes an update for the particle positions, Ax", and particle velocities, Av", for given discrete time-step At, which yields:
x++ X Ax
1+ v¹ + Av'.
Existing neural solvers (Ummenhofer et al., 2019, Prantl et al., 2022) often simplify this problem by inferring the velocity update from the position update, i.e., Av* Ax' Ax//At. /tt Including initial velocity v and external constant forces a', such as gravity, yields our general learning task as
++ x + G(x', O); x² x""xx" Alv' + AlPa',
where G is a Neural Network with parameters θ and feature vector f. Note that the initial velocity v", due to the inference of velocity from position updates, not the actual instantancous velocity of the particles but is defined via v = x xx-x'-r The minimization task is then given
min L(x G(x f O), y*
(3)
where y are known, ground-truth positions from a given dataset and L being a loss function, commonly chosen as L In the following, we refer to solvers of this form as Neural Network Time Integrators (NNTIs). Note that the time-step At and associated time-stepping scheme do not need to be the same for ground-truft simulation and NNTI. For example, it possible to predict multiple ground auth-simulation steps simultaneously using multi-stepping or hofer et al. 2019; Brandstetter et al. 20226).
(Ummen-)
2.2 CLASSICAL SOLVERS
We utilize the Smoothed Particle Hydrodynamics (SPH) method the basis for test cases and to inform some of our inductive biases. SPH is a Lagrangian simulation technique initially intro- duced in astrophysics context (Monaghan 1994) but has found broad application in various fields, e.g., in CFD (Monaghan) 2005) and Computer Graphics (Thmsen et al.] 2013 Macklin & Müller 2013). the heart of SPH are interpolation operations (A) of quantities A carried by particles, with positions x, mass m and density p. using a Gaussian-like kernel function W as
(A,)
Σ Aj mh m
W(x; - -x,,,xx,,) h,,)),
JEN.
where h is the support radius and Nf being the neighbors of i, including itself, which are all particles closer than h, see Koschier al. (2019) for broader overview. This can be seen message- passing step using the edge lengths and features of the adjacent vertices with a summation message- gathering operation. Within SPH many solvers exist for a broad variety of problems and choosing
2

---

Published as a conference paper at ICLR 2024
the correct one for a respective problem can be challenging they are vastly different. As our test cases involve compressible, compeey-om and incompressible SPH solvers to generate the data, understanding their background and requirements is valuable.
For compressible simulations, SPH utilizes either an explicit pressure formulation using a com- pressible Equation of State, do we, using Riemannian solvers (Puri & Ramachandran 2014). which would be an important direction for future research due to their complexity. For weakly com- pressible simulations the most commonly utilized technique is the 8-SPH method (Marrone et al.) 2011) using explicit pressure forces combined with diffusion models for velocity and density terms using very small timesteps, with the 8P-SPH (Sun et al., 2018 variant that we use for our data generation further expanding this approach. Finally, for incompressible SPH (Ihmsen et al.) [013) Benderr& & Koschier, 2013) the simulation revolves around solving a Pressure Poisson Equation us- ing an implicit solver using many iterations per timestep. Overall, the SPH solvers we used vary from straightforward explicit integration over numerically challenging explicit integration with very small timesteps to requiring large linear systems to be solved per timestep.
2.3 NEURAL NETWORKS
Many methods have been proposed to solve PDEs using machine learning techniques BBttaglia et al. (Mrton
(201)); |Morton et al. (2018); Um et al. (2020), where Physics-Informed Neural Networks (PINNs) are among the most prominent approaches to solving such PDEs (Cai et al.) 2021]. PINNS are coordinate networks that learn the solution of a PDE without first requiring discretization using continuous residuals. While powerful, these methods have many drawbacks, including difficulties in training and cannot be applied to discrete simulation data. Another prominent machine learning technique is point-cloud classification and segmentation, e.g., PointNer (Qi et al.] 2017). However, these approaches are not well-suited to Lagrangian simulations as particles
much more inter-
connected than point clouds. Neural Operators have recently also found research interest in solving PDEs (Li al. 2021b] Guibas al., |021 [i et al. 2021a). However, applying these approaches to irregular grids or particle-based simulation is not easily possible. On the other hand, Graph Neural Networks (GNNs) (Sanchez- Gonzalez al. 2020) naturally map to simulations as most simulations can be interpreted as form of message passing on graph.
Graph Neural Networks can be directly applied SPH simulations, where each particle is a ver- tex, and the particle neighhorhoods describe the graph connectivity. GNNs use the graph to generate messages using just the vertex information (Qi et al. 2017), edge information (Wang et al. 2018), edge and vertex (Sanalez-Gonzzlez et al.) 2020), or edge, vertex, and additional feed-through fea- tures collected on vertices using pooling operations (Brandstetter et al., 2022a). These collected messages are then either used directly, as new features on the vertices, combined with existing vertex features. Further operations, such as input encoding and output decoding, have also been proposed (Sanalez-Gonzalez Gonzalez et al., 2020). The message processing performed using MLPs with relatively shallow but broad hidden architectures, e.g., 2 hidden layers with 128 neurons.
Continuous Convolutions are a subset of GNNs and utilize only coordinate distances as inputs to the filter functions, similar to SPH kernel functions (Ummenhofer et al. 2019). These filter functions are then combined with the features of adjacent vertices to form the messages. While this imposes an inductive bias, it can make learning the problem more manageable. The coordinate distances can then be processed using an MLP (Wang et al. 2018) or other function interpolation techniques, e.g., linear interpolation (Ummenhofer et al. 2019) or spline- based interpolation (Fey et al.) 2018). Several extensions of these approaches have been proposed for physical simulations, e.g., using antisymmetry (Prantl et al. 2022) to conserve particle momentum.
Fourier and Chebyshev Methodes: Fourier Neural Operators (FNOs) (Li et al. 202Tb) learn phys- ical simulations by transforming given spatial discretization into spectral representation using an FFT. By applying the learning task in the spectral domain, these operators can learn spatially in- variant behavior but are limited to regularly sampled data due to their reliance on FFTs. Fourier encodings have also been applied in image classification and segmentation to make networks less spatially dependent (Li et al. 202 Chebyshev basis polynomials have also been used in Graph classification tasks using CConvs as higher-order interpolants (Defferrard et al. 2016 He et al. 2022), as well other polynomial bases, e.g., Bernstein polynomials (He et al. 2021
3

---

Published as a conference paper at ICLR 2024
3 METHOD
Our work builds on a generalized formulation of continuous convolutions that acts as a superset of existing methods, such as LinCConv (Ummenhofer et al., 2019) and SplineConv (Fey al. 2018). Based on this model, we describe parameterizations of convolutional filter functions using separable basis functions and explain how prior work fits into this concept. Furthermore, we will describe how symmetries can be built into this model and construct a Fourier-basee convolutional architecture that uses both even and odd symmetry. Finally. we will discuss window functions, an important inductive bias in many existing CConv approaches.
Formulation: In general, convolution a mathematical operation that combines an input function, f: IR- + R. and filter function, g R + R. through an operation defined as (Wang et al. 2018)
(O / sa - T) g(T)dd..
(5)
We then limit g to be compact, i.e., g()) V/T 2 h, with h being cutoff distance, also referred to as support radius within SPH contexts. We can then sample I on the positions of vertices, x, and base 7 on the coordinate distances between connected vertices to discretize the convolution as
(f (fog) (
f(x)) JENS
X; X,
(6)
where: and j are the indices of two vertices, see Appendix A.1 This formulation is then expanded by the inclusion of a normalization term o(x), typically only used in classification tasks, window function w similar in shape to an SPH kernel function, see Appendix (A.2 and coordinate mapping function A Rd -> Rdd see Appendix A.3 Denoting the trainable weights of y by θ this yields
X: h (A(<)=*)))
Xx; [x;=xs1]
x51
(f (fxg)(x)= ((x,)
f(1)) ge JEN.
h
The machine learning task then is to find set of weights & such that (f ge)((x) yxx),, where y represents the supervised ground truth result. We now propose to parametrize ge using a set of one-dimensional basis functions bz(r), with i E [0. n), such that for one-dimensional convolution we get (g)g) = (b(4), O) where (.-)) denotes the inner product. A direct choice for b; would be piece-wise constant function that results in a Nearest Neighbor interpolation or pisce-wise linear function that results in the LinCCony approach; see Appendix [A.4] For a two-dimensional convolution, construct a matrix of basis terms as the outer product of set of separable basis terms in r and y. i.e., by (qx), with E [0, u), and by (qyy) with j E [0, v), such that
(8) While the imposed restriction of the basis terms being separable limits the potential choices for basis terms, it is also in line with prior work, e.g., by Fey et al. (2018). A key benefit of such a separable be computed straightforwardly through transpose of the weight
(g) = (bb((.)) @ ())))),
formulation is that gradients matrix instead of requiring more expensive steps or even involving non-linear optimization, as is required for traditional RBF Networks as proposed by Broomhead & Lowe (9988). It is important to keep in mind that, as continuous convolutions are matrix multiplications, updating + at learning time is still linear operation even if the basis terms are non-lineat see Appendix A.
Several basis terms b have already been considered in prior work, and we will summarize some of them here in ww-dimensional context. A simple choice is using a bi-linear interpolation, e.g., as in the LinCConv approach (Ummenhofer et al. 2019). Here, each entry of the weight matrix O only has a very local influence, and the shape the learned convolution operator go a combination of picce-wise linear functions, and thus also picce-wise linear. While such formulation is straight- forward to implement, it gives no guarantees regarding symmetry, the learned function cannot be smooth, and all learned convolutional filters in larger network will have discontinuities located at the same relative positions. Using cubic B-Splines (Fey et al. 2018) results in smoother learned filters without discontinuities but still does not guarantee any symmetries.
Incorporating Symmetry: The DMCF approach (Prantl et al. 2022) used an amtisymmetric basis to enforce conservation of momentum, which was implemented through explicit mirrrinn of weights.
4

---

Published as a conference paper at ICLR 2024
Figure 2: Quantitative evaluation of the mne-dimensional kernel function (left) and kernel gradient (right) toy problem. Evaluations are based on the L2 error related to the number of base terms.
This had the side effect of reducing the effective number of parameters by factor of two, and symmetry was neglected during backpropagation. Our framework allows us to directly implement the antisymmetry constraint, which enables us to instead modify the basis formulation from Eq. such that it can be applied to any basis function and results in an antisymmetric basis
le
asymm (q) = ([sg((g:) ( )b (2)qzl - 1) @ b" (sgn(qq))
qy]], 9),
(9)
which can be modified to lead to a symmetric formulation by excluding the sgn term:
symm (g) ([b" qa 1) @ by (sgn(q.)) .))]. 9
(10)
While ensuring that all basis functions are antisymmetric is useful for conserving momentum, not all target functions are necessarily antisymmetric, e.g., a density interpolation in SPH is rotationally invariant and can not be learned with such basis. Furthermore, the resulting filter function is not smooth, and the resolution along different coordinate axes is not identical. To resolve these short- comings, we propose set of smooth basis functions with either basis terms influence the outcome for any value q. Our primary choice for such a basis is a Fourier series; see Appendix A.4 for visualization and definition of the two-dimensional Fourier series. Note that this fundamentally differs from applying a Fourier transform, e.g., as done in the FNO approach (Li et al.) 2021b). Architectures like FNO transform an input signal into frequency space through an explicit Fourier transform and learn in the frequency domain. We instead use the Fourier basis to construct a convolutional kernel. The input signal keeps its spatial representation, but the learning task becomes finding the best possible coefficients for a Fourier series. Finally, we also in- clude Chebyshev polynomials, which are popular for approximating non-perinnn non-periodic compact functions and have inherent symmetries. Additional basis constructions are discussed in Appendix A.4
odd symmetry, where all
Window Functions are an important inductive bias in many prior works regarding CConvs, espe- cially within learning physical simulations. The general motivation behind a window function is that interactions should be (a) compact, (b) smooth, and (c) behave like SPH interpolations. This induc- tive bias is primarily informed by SPH methodology and, consequently, prior work generally used SPH kernel functions as window functions such as the Müller (Ummenhofer et al. 2019) or Spiky kernel (Prantl et al. 2022). We evaluate several other functions; details given in Appendix A.2
4 RESULTS
We now propose the SFBC (Symmetric Fourier Basis Convolution) approach using a Fourier basis with no window function and an identity coordinate mapping. We first compare SFBC against other CCony approaches in a toy problem in one and three dimensions to compare the capabilities of the interpolation bases; see Section 4. We then compare SFBC against various baselines in compressible one-dimensional problem; see Section 4.2 Next, we perform an in-depth evaluation of a two-dimensional closed domain simulation focused on inference stability; see Section 4.3 Finally, we evaluate a fluid blob collision scenario to investigate how different basis terms perform with partially occupied support domains; see Section [.4] For more details on the training and setup, Appendix C We also performed a runtime analysis of the most relevant hyperparameters and found that SFBC on average only incurs an increase of 0.5%, 8.5% and 12% over LinCConv in one, two, and three dimensions, respectively, see Appendix C.6 and Figure 38 for details.
4.1 TOY PROBLEMS
To evaluate the capabilities of the different basis functions to learn symmetries, we consider two tasks: (a) symmetric task to learn an SPH kernel interpolation, and (b) an antisymmetric task
5

---

Published as a conference paper at ICLR 2024
18
16
Figure 3: A quantitative evaluation of the relationship between parameter count and test error for basis convolutionsleet),, and a quantitative evaluation of different networks for fixed layout with four message- message-passing steps and 32 features per layer (right). Error bars showing lower to upper 5%.
to learn the SPH gradient, see Appendix B.2 This setup is motivated by the concept that if a machine learning method cannot learn the basic components of an SPH simulation, then learning the overarching SPH simulation is made more difficult. As we want to evaluate the abilities of the different basis functions to act as interpolation functions, we utilize a network with a single message- passing step without activation functions, i.e., the learning task is a linear optimization problem.
Based on the results shown Figure 2 and Appendix C.I.I. we observe clear difference between different basis terms. The SFBC approach performs best but shows decrease in performance with increasing numbers of base terms as this learning task does not require any higher-order harmonics, and learning contribution of zero is potentially challenging. The second best-performing method is the Chebyshev basis, although this method performs several orders of magnitude worse than the Fourier basis. The baseline methods show much worse overall perfor- mance, only getting closer SFBC for a high number terms. The DMCF approach and a Fourier series consisting of odd terms perform the worst these, by construction, can only learn antisym- metric behavior. The inclusion of symmetry significantly improves performance and reduces the lowest achieved L2 error by a factor of 5217 when compared to LinCConv.
Kernel Function:
For this case, shown in Figure 2 and Appendix C.I.2, we see a notably different result. Only few of the methods demonstrate any reasonable performance, and all of these meth- ods have inherently antisymmetric terms. Overall, the Fourier series with only odd symmetry terms performed best, followed by complete Fourier series and DMCF, where the latter was impeded primarily by lack of smoothness for low basis term counts. These results demonstrate that an- tisymmetric learning tasks are challenging for traditional basis functions, and including symmetry and smoothness as an inductive bias significantly improves the overall learning behavior by factor of 7.6 relative to DMCF and 853 relative to LinCConv when compared to the odd Fourier series.
Kernel Gradient:
We expanded this evaluation to three dimensions, Appendix B.5 where we found similar behavior regarding symmetric tasks with SFBC outperforming LinCConv by a factor of two average. As the tensor products of antisymmetric bases not necessarily antisymmetric, we observed much worse performance for the antisymmetric odd Fourier basis. They perform on average three times worse than SFBC (see Appendix C. We also utilized this setup to evaluate how non-ideal learning setups perform and found that SFBC was significantly more resilient against superfluous message-passing steps and input features, see Appendix C.5.2 Overall these findings suggest that including symmetry and smoothness into the basis significantly improves learning per- formance and improves the networks ability to learn in broad range of conditions.
Three Dimensions:
4.2 TEST CASE I: COMPRESSIBLE ID SPH
Based on the results from the toy problems, we now evaluate how these behaviors translate to more holistic learning task where the overall physics update should be learned. The learning setup here a one-dimensional compressible SPH simulation, where the updated velocities per particle should be learned based on the velocity of the current particles. Appendix[B.I [B. provides details of the data generation and underlying simulation model. For the width of the base functions, chose 6, based on the results from the toy problems and a hidden architecture for the MLP-based approaches of 2 deep and 128 wide, in line with
al.] (2020).
Basis Function Methods: We first consider the influence of the network size on the achievable
test error by evaluating large number of hyperparameters, Fig. 3 and Appendix C.2 Our

---

Published as a conference paper at ICLR 2024
Figure 4: A qualitative comparison after 64 inference steps with (f.l..r.) of the ground truth data, our SFBC approach, LinCComv, Fourier basis with window, and Chebyshev with window. The particle data is mapped to a regular grid with color mapping indicating current velocity.
Density
Velecity
Divergence
10
10
19*
Nearest hegrbor
LincCone
Dants
Chetssher Chelysher (2nd knn))
10
SplineCons
STac
Geunaden Kema
Derc
Kuller Karnel
Saley
157
Wend und
Bump RDP
Quarte Sp in
10
Figure 5: The quantitative results are based on the integrated steps regarding four error metrics not included in the training loss.
inference period of 64
results demonstrate that there is clear and noticeable difference between basis functions built upon antisymmetric terms, i.c., Fourier, Fourier (odd) and DMCF, that makes them perform notably better than all other basis functions. Considering the size of the network, made several observations; see Appendix C.2 On the one hand, we found that increasing the number of message-passing steps, e.g., from 3 to 6, only changed the network performance by 100, while increasing the features per layer from 4 to 32, with 3 message-passingg steps, improved performance by an order of magnitude, Fig. 21 On the other hand, also saw clear benefit of the SFBC approach it outperforms DMCF at virtually all sizes, with more significant difference for smaller layouts.
MLP-based Methods: Next we compare to set of MLP-based GNNs using PointNer (Qi al., 2017). GNS Sanchez-Gonzalee et al. 2020), MLSCom (Wang et al.] 2018) and MP-PDE (Brand- stetter et al. 2022b) architectures with message passing steps and features per vertex and mes- sage. Our results show that GNS and MP-PDE perform as well as most convolutional basis functions but not as well as methods with built-in symmetries, see Fig. 3 However, they require significantly more parameters to reach comparable accuracy, i.e., MLSConv, GNS, and MP-PDE require 229K, 393K, and 395K parameters, respectively, compared to 30K for the CConv methods. This highlights the importance of basis convolutions to achieve the same performance with fewer resources. Within the collection MLP baselines, MI.SCovv performs slightly worse but is still comparable other approaches. This indicates that the inductive bias of including a convolution by itself is not sufficient but that the advantages come from constructing them using basis functions with inherently useful properties.
inductive bias that allows the CConv-based networks
Overall, the results indicate that the basis function convolutions perform similarly to MLP-based GNNs while requiring significantly fewer parameters. Furthermore, including a bias of symmetry significantly improves the capabilities of network, e.g., only methods with symmetry were able accurately learn the gradient function. Overall, our SFBC approach shows its robustness by per- forming better than other methods across a wide range of architectural changes, see Appendix
4.3
TEST CASE II: WEAKLY-COMPRESSIN COMPRESSILLE 2D SPH
For the second test case, we focus on the long-term stability of various NNTIs for a closed domain weakly-compressihhe simulation; see Appendix B.3 for details on the data generation. We utilize a network architecture using n 4 with four rounds of message passing, in line with (Ummenhofer ct al. (2019), and train the networks with a maximum rollout of 10 and an evaluation on the test dataset with an inference length of 64. To quantify the performance, we map the particle density and velocity to a regularly sampled grid spanning the closed simulation domain and then assess the L2 difference of density, velocity, and divergence on the grid. We also compute Power Spectral Den- sity (PSD) difference based on the velocity field. The divergence is the central metric for this case it derived metric that is challenging uphold and closely correlated to lack of smoothness,

---

Published as a conference paper at ICLR 2024
e.g., see Fig. Note that naturally require low error in all metrics for truly accurate result. We now evaluate a broad set of choices of basis functions and hyperparameters.
Baseline Comparison:
The proposed Fourier- Fourier-based network performs notably better than other baseline methods, i.., LinCCowv DMCF and SplineConv, and outperforms all other basis functions in almost all metrics, e.g., it outperforms LinCCon by factor of 2.6 regarding divergence enror. The only exception is the nearest neighbor basis regarding density and PSD error but, crucially, not regarding divergence error; see Fig. 5 and Appendix [C3.I] Overall, Chebyshev basis functions also performed well, but only when using a window function, and Spline performs worse than LinCCom, in line with prior observations (Ummenhofer et al. 2019). Finally, for B-spline and other SPH kernels, the higher the order, the better the result. In addition to the quantitative evaluations, we also considered the qualitative results, where the Fourier basis with no window shows the, by far, best prediction, see Fig. In contrast, other methods show superimposed noise on the velocity field, which is combination of both the choice of basis and window function, adding different noise patterns. This superimposition of noise can easily dominate the overall prediction and make the prediction unstable and unusable. We also performed study for more significant number of seeds (32 instead of 4) for the key methods, Appendix C.3.5 and found no significant difference, ie., the divergence error for our method changed from 1.29 to 1.30.
Window Function:
We also evaluated different window functions regarding both the LinCCowv baseline and our proposed basis, Fig. 6 We observed that the choice of window function significantly impacts the network performance, e.g., LinCConv improved in terms of den- sity error by 2.1 times; see Appendix C.3.2 At the same time, the Fourier basis shows 2 times increase for the divergence eтror with the Müller window. This points to fundamental differences stemming from the hasis con- struction. Notably, the Fourier basis has the desirable property of giving a very good performance without requiring a window function. Finally, our eval- uations demonstrate that contrary to intuitions in prior work, using window functions that are not Gaussian-shaped, e.g., a linear window, can outperform existing window functions, e.g., using a linear window showed the lowest density, velocity, and PSD errors for LinCConv. This highlights the importance of choosing the window, as the network's s performance can be fine-tuned for a task.
Divergance
Beste LaiCCrv
Sasie SC
Widdow
Noee
Salky
Guaric Splive
Figure 6: Window Function Study
Fourier Terms: We already observed significantly dif- ferent behavior for different choices of Fourier series terms in the toy problems and now investigate this more closely. We evaluated different variations of Fourier- based networks; see Fig. [ and Appendix C.3.4 where we made several crucial observations. Using only even or odd symmetry basis terms did not lead an overall stable prediction, highlighting that using either symme- try exclusively is not ideal. Furthermore, by changing which harmonies are used for given number of terms, the behavior can be adjusted to be optimal for using window function or not. Replacing the first harmonic cosine term with second harmonic sine term without a window function reduced the divergence error by factor of 1.4.
Diowgence
Sam
Fourte Saries Poureer lord)
Figure 7: Fourier Series Study
Coordinate Mappings: [Ummenhofer et al. (2019) proposed a volume-preserving mapping in the LinCConv approach. To expand the brief evaluations in previous work, used framework to assess the importance of the coordinate mapping. In addition to the volume-presen preserving mapping, mapping via the classic choice of polar coordinates serves baseline, details for which can be found in Appendix CC3.3) Comparing these two variants to networks trained without any coordinate mapping shows no clear advantage as the divergence error increased by up 15%.
Overall, our evaluations show a clear and significant improvement over existing baselines for our SFBC approach. Not only are quantitative metrics improved, e.g., the divergence error is reduced by 30,, but these results were also achieved with fewer inductive biases.
8

---

Published as a conference paper at ICLR 2024
4.4 TEST CASE III: INCOMPRESSIBLE 2D SPH
For our final test case, we focus on the behavior of different basis terms regarding free surfaces for an underlying divergence free solver; see Figure 9 and Appendix C.5 We use the L2 difference of the ground truth and predicted particle positions for training and compute the mean distance of all pre- dicted particle positions the closest particle in the ground truth data for evaluation.
75
Nearest Neighoor
5.
Undtonv
SalineCorv
Chebyshev
0.0
DMCF
Inference
Stop
Point-wise distances for different inference lengths as mean values per step.
Figure
Basis Function Comparison: We use architecture with 6 basis terms, message-passing steps and 32 features to compare the most relevant baselines. The Fourier-based method performs best during the initial inference period, up to times the training rollout length, e.g., 16 inference steps. Still, after that, most methods perform more similar; see Fig. [ and 31 Qualitatively. the Fourier ba- sis performed noticeably better than other methods in preserving the shape of the colliding droplets, as shown in Fig. 9 The window function in this scenario, see Appendix [.4.3] has a significant impact on performance for prior approaches, i.e., LinCConv showed a 3.3x difference by changing the window function, and DMCF was only stable with the Müiler window. In contrast, Fourier and Chebyshev bases exhibited a 20% difference.
Clounciiih
To verify whether the benefits of the Fourier approach hold across varied architectures, we repeated this experiment across different base term counts and network layouts, see Appendix C.4.I and C.4.2 This evaluation shows that the same scaling behavior occurs in test case I, implying that this behavior holds across dif- ferent dimensions and problems. Furthermore, we found the same benefit for our SFBC approach, i.c., there is a notable improvement virtually all network sizes with more pronounced benefit of up to 50% for smaller lay- outs, indicating that SFBC provides robust and accurate basis for learning representations of physical systems.
Varying Network Size:
Kas + Window
Figure 9: Qualitative evaluation of test case starting on frame 10 for 32, 64 and 96 timesteps, velocity color coded.
Overall, our SFBC approach works best in this very chal- lenging scenario and outperforms all baselines. We fur- thermore observed that some baseline methods, especially DMCF, performed vastly differently when using no, simply different, window functions. In contrast, approach remains stable in all cases while performing best without a window function.
CONCLUSIONS
We introduced the Symmetric Fourier Basis Convolution (SFBC) approach as a novel inherently symmetric and smooth continuous convolution approach applied to Lagrangian fluid simulations. Using this technique, we found improved performance our challenging test cases, representing different fluid systems and solvers. Our broad evaluations identified prior inductive biases that are no longer necessary for our Fourier-based approach. At the same time, they confirm the general benefits of learning function-based convolutions unstructured settings.
However, while we did consider a broad set of parameters, we only considered networks with up circa 200 thousand parameters, which is very relevant topic for future work. Our framework allows for easy and flexible explorations of continuous convolutions with new basis functions. We hope our work will inspire future investigations to improve learning methods for unstructured data sets outside of fluid mechanics (Lam et al.) 2022 Reiser et al.) 2022 AAbeet--rissizzball et al.) 2021). rephrased: Furthermore, we would like to expand our formulation to encompass more general non-linear versions of Radial Basis Function networks (Broomhead & Lowe 1988). In addition, algorithms for the initialization of the new types of basis networks, and existing methods such as LinCConv, an where substantial room for improvement exists. Finally, we would like to explore the larger space of basis functions combining different functions for different tasks, e.g., using a different basis for the radial and angular components for spherical coordinate mapping.

---

Published as a conference paper at ICLR 2024
ACKNOWLEDGMENTS
This work was supported by the DFG Individual Research Grant TH 2034/1-2.
ETHICS STATEMENT
Numerical solvers for PDEs can be applied in many fields and play an essential role in design and engineering. Numerical solvers can also be utilized in other less beneficial areas, but our research is connected to any ethically questionable application fields. Predicting the societal impact of more efficient solutions numerical systems is challenging. However, such methods can play positive role either directly by reducing the Carbon footprint running numerical solutions or indirectly by solving inverse problems.
REPRODUCIBILITY STATEMENT
We generated all datasets in our work by ourselves and using custom solvers. Consequently, it is of vital importance that the datasets and solver implementations are freely available, and. accordingly, we have made the dataset with the solver and neural network implementation available open source code at fhttps: /ttpps:/gitee tiggtthu.... wwt.///mmmmsssssseeee /SFBC The open-source code also includes implementations of baselines with configurations for the respective hyperparameters used in our paper.
Outside of the repository, we have also described our evaluations in Section C provided further implementation details in Appendix Section A We have also included detailed descriptions of our data generation setups in all test cases in Appendix Section B We also performed a broad set of hyperparameter studies to ensure that our baseline evaluations are fair and proper, see Section C.3.5
REFERENCES
David
Mohammad Ali Armin, Simon Denman, Clinton Fookes, and Lars Pe- tersson. Graph-based deep learning for medical diagnosis and analysis: past, present and future. Sensors, 21(14):4758, 2021.
Kelsey Allen, Tatiana Sopez-Guevara, Kimberly Stachenfeld, Alvaro Sanzhez- Pe- ter W. Battaglia, Jessica B. Hamrick, and Tobias Pfaff. Inverse design for fluid-structure inter- actions using graph network simulators. In Advances Neural Information Processing Systems, 2022.
Matteo Antuono, Andrea Colagrossi, and Salvatore Marrone. Numerical diffusive terms in weakly- compressible sph schemes. Computer Physics Communications, (83(12):2570-22 25802 2012.
Yohai Bar-Sinai, Stephan Hoyer, Jason Hickey, and Michael P Brenner. Learning data-driven dis- cretizations for partial differential equations. Proceedings of the National Academy of Sciences, 11663))15344415549 2019.
Peter Battaglia, Razvan Pascanu, Matthew Lai, Danilo Jimenez Rezende, et al. Interaction networks for learning about objects, relations and physics. Advances in neural information processing systems, 29, 2016.
Jan Bender and Dan Koschier. Divergence-free smoothed particle hydrodynamics. In Proceedings of the 14th ACM STGGRAPH//Enographis symposium on computer animation, 147-155, 2015.
Johannes Brandstetter, Max Welling, and Daniel E Worrall. Lie point symmetry data augmentation for neural pde solvers. In International Conference on Machine Learning, 2022a.
10

---

Published as a conference paper at ICLR 2024
Johannes Brandstetter, Daniel E. Worrall, and Max Welling. Message passing neural PDE solvers. In The Tenth International Conference on Learning Representations, ICLR. OpenReview.net, 2022b.
Robert Bridson, Jim Houriham, and Marcus Nordenstam. Curl-noise for procedural fluid flow. ACM Trans. Graph., 26(3):46-es, jul 2007. ISSN 0730-0301. doi: 10.1 10.1145/1277377.1777433
David S. Broomheed and David Lowe. Multivariabll functional interpolation and adaptive networks. Complex Syst. 2, 1988.
Shengze Cai, Zhiping Mao, Zhicheng Wang, Minglang Yin, and George Em Kariaaaaii. Physics- informed neural networks (pinns) for fluid mechanies: A review, 2021. URL https://arxi:: org/ hte//bb//a 2105. 09506
Michaël Defferrard, Xavier Bresson, and Pierre Vandergheynst. Convolutional neural networks on graphs with fast localized spectral filtering. Advances in neural information processing systems, 29, 2016.
Walter Dehnen and Hossam Aly. Improving convergence in smoothed particle hydrodynamics sim- ulations without pairing instability. Monthly Notices of the Royal Astronomical Society, 425(2): 1068-1088 1082, 2012.
Matthias Fey, Jan Eric Lenssen, Frank Weichert, and Heinrich Müller. Splinecnn: Fast geomet- ric deep learning with continuous b-spline kernels. In P Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 869-877, 2018.
Jens André Griepentrog, Wolfgang Höppner, Hans-Christoph Kaiser, and Joachim Rehberg. A bi- lipschitz continuuus, volume preserving map from the unit ball onto cube. Note di Matematica, 28( 281)111771 2008.
John Guibas, Morteza Mardani, Zongyi Li, Andrew Tao, Anima Anandkumar, and Bryan Catanzaro. Efficient token mixing for transformers via adaptive fourier neural operators. In International Conference on Learning Representations, 2021.
Mingguo He, Zhewei Wei, Zengfeng Huang, and Hongteng Xu. Bernnet: Learning arbitrary graph spectral filters via bernstein approximation. In Mare' Aurelio Ranzato, Alina Beygelz- imer, Yann N. Dauphin, Percy Liang, and Jennifer Wortman Vaughan (eds.), Advances in Neural Information Processing Systems, pp. 10239-14255 2021.
Mingguo He, Zhewei Wei, and Ji-Rong Wen. Convolutional neural networks on graphs with cheby- shev approximation, revisited. Advances in Neural Information Processing Systems, 35:7264 7276, 2022.
Markus Ihmsen, Jens Cornelis, Barbara Solenthaler, Christopher Horvath, and Matthias Teschner. Implicit incompressible sph. IEEE transactions visualization and computer graphics, 20(3): 426 435, 2013.
Diederik Kingma and Jimmy Ba. Adam: method for stochastic optimization. In Yoshua Bengio and Yann LeCun (eds.), 3rd International Conference on Learning Representations, ICLR, 2015.
Dan Koschier, Jan Bender, Barbara Solenthaler, and Matthias Teschner. Smoothed particle hy- drodynamics techniques for the physics based simulation of fluids and solids. In Wenzel Jakob and Enrico Puppo (eds.), 40th Annual Conference of the European Association for Computer Graphics, Eurographics 2019 Tutorials, pp. 1-41. Eurographics Association, 2019. doi: 10.23 10.23122/et222111033
Remi Lam, Alvaro Sanzhez-Gonzalez, Matthew Willson, Peter Wirnsberger, Meire Fortunato, Alexander Pritzel, Suman Ravuri, Timo Ewalds, Ferran Alet, Zach Eaton-Rosen, Weihua Hu, Alexander Merose, Stephan Hoyer, George Holland, Jacklynn Stott. Oriol Vinyals, Shakir Mo- hamed, and Peter Battaglia. Graphcast: Learning skillful medium-range global weather forecast- ing. 2022.
Yang Li, Si Si, Gang Li, Cho-Jui Hsieh, and Samy Bengio. Learnable fourier features for multi- dimensional spatial positional encoding. Advances in Neural Information Processing Systems, 2021a.
11

---

Published as a conference paper at ICLR 2024
Zongyi Li, Nikola Borislavov Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhat- tacharya, Andrew M. Stuart, and Anima Anandkumar. Fourier neural operator for parametric par- tial differential equations. In 9th International Conference on Learning Representations, 2021b.
Mario Lino, Stathi Fotiadis, Anil A Bharath, and Chris D Cantwell. Multi-scale rotation-equivarian graph neural networks for unsteady eulerian fluid dynamics. Physics of Fluids, 34(8), 2022.
Miles Macklin and Matthias Müller. Position based fluids. ACM Transactions on Graphics (TOG), 32(4):1-12, 2013.
Salvatore Marrone, MAGD Antuono, A Colagrossi, G Colicchio, Le Touzé, and G Graziani. 8- sph model for simulating violent impact flows. Computer Methods in Applied Mechanics and Engineering, 200(13-16):1526-1544 2011.
Joe J Monaghan. Simulating free surface flows with sph. Journal of computational physics, 110(2): 399-406, 1994.
Joe Monaghan. Smoothed particle hydrodynamics. Reports on progress in physics, 68(8):1703, 2005.
Jeremy Morton, Antony Jameson, Mykel J Kochenderfer, and Freddie Witherden. Deep dynamical modeling and control of unsteady fluid flows. In Advances in Neural Information Processing Systems, 2018
Matthias Müller, David Charypar, and Markus Gross. Particle-based fluid simulation for interactive applications. In Proceedings of the 2003 ACM SICGRAPH//Enograhhis symposium on Com- puter animation, pp. 154-159. Citeseer, 2003.
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang. Zachary De Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. Pytorch: An imperative style. sigh-performance deep learning library. In Advances Neural Information Processing Systems 32, pp. 8024 8035. 2019.
Ken Perlin. An image synthesizer. In Proceedings of the 12th Annual Conference on Computer Graphics and Interactive Techniques, SIGGRAPH '85, pp. 287-296, New York, NY, USA, 1985. Association for Computing Machinery. ISBN 0897911660. doi: 10. 10.1145/253343325247.
Tobias Pfaff, Meire Fortunato, Alvaro based simulation with graph networks. In 9th International Conference on Learning Representa- fions, ICLR 2021, 2021.
and Peter W. Battaglia. Learning mesh-
Stephen B. Pope. Turbulent Flows. CB009700111840531.
Cambridge University Press, 2000. doi: 10.1017/
Lukas Prantl, Benjamin Ummenhofer, Vladlen Koltun, and Nils Thuerey. Guaranteed conserva- tion of momentum for learning particle-based fluid dynamics. Advances in Neural Information Processing Systems, 35:6901-6913, 2022.
Daniel Price. Smoothed particle hydrodynamics and magnetohydrooynamiis. Journal of Compu- tational Physics, 231(3):759-794 2012.
Kunal Puri and Prabhu Ramachandran. Approximate riemann solvers for the godunov sph (gsph).
Journal of Computational Physics, 270:432-455 2014.
Charles R Qi, Hao Su, Kaichun Mo, and Leonidas J Guibas. Pointnet: Deep learning on point sets for 3d classification and segmentation. In Proceedings of the IEEE conference computer vision and pattern recognition, pp. 652-660, 2017.
P Rastelli, R Vacondio, JC Marongiu, G Fourtakas, and Benedict D Rogers. Implicit iterative particle shifting for meshless numerical schemes using kernel basis functions. Computer Methods in Applied Mechanics and Engineering, 393: 593:11477 14716, 2022.
12

---

Published as a conference paper at ICLR 2024
Patrick Reiser, Marlen Neubert, André Eberhard, Luca Torresi, Chen Zhou, Chen Shao, Houssam Metni, Clint van Hoesel, Henrik Schopmans, Timo Sommer, et al. Graph neural networks for materials science and chemistry. Communications Materials, 3(1):93, 2022.
Alvaro Sanchez-Gonzales Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovee, and Peter Battaglia. Learning to simulate complex physics with graph networks. In International conference on machine learning, pp. 8459-8468. PMLR, 2020.
Víctor Garcia Satorras, Emiel Hoogeboom, and Max Welling. E(n) equivariant graph neural net- works. In Marina Meila and Tong Zhang (eds.), Proceedings of the 38th International Conference on Machine Learning, pp. 9323-9332, 18-24 Jul 2021.
Barbara Solenthaler and Renato Pajarola. Density Contrast SPH Interfaces. In Eurograph- ics/SIGGRAHH Symposium on Computer Animation. The Eurographics Association, 2008. ISBN 978-3-905674-10-1. doi: 10.2312/SCASSCA0//11-22
PN Sun, Andrea Colagrossi, Salvatore Marrone, Matteo Antuono, and AM Zhang. Multi-resolution delta-plus-sph with tensile instability control: Towards high reynolds number flows. Computer Physics Communications, 224:63-80, 2018.
Vijay Thakkar, Pradeep Ramani, Cris Cecka, Aniket Shivam, Honghao Lu, Ethan Yan, Jack Kosa- ian, Mark Hoemmen, Haicheng Wu, Andrew Kerr, Matt Nicely, Duane Merrill, Dustyn Blasig, Fengqi Qiao, Piotr Majcher, Paul Springer, Markus Hohnerbach. Jin Wang, and Manish Gupta. CUTLASS, January 2023. URL https::/tttttt.t http:/ggttttub.... com CICCONNNVIIIAA.. http://VIIAAA/utlsssssaeesss catlass
Nathaniel Thomas, Tess E. Smidt, Steven Kearnes, Lusann Yang. Li Li, Kai Kohlhoff, and Patrick Riley. Tensor field networks: Rotation- and translation-eqquivrriann neural networks for 3d point clouds. CoRR, abs/ abss/1001 1802.08219, 2018.
abs /ts///110 /w////102.000 08219
Nils Thuerey, Konstantin Weißeenww, Lukas Prantl, and Xiangyu Hu. Deep learning methods for reynolds-averagee navier-stokes simulations of airfoil flows. AIAA Journal, 58(1):25-33) 2020.
doi: 10.2514/1 J055829..
Jonathan Tompson, Kristofer Schlachter, Pablo Sprechmann, and Ken Perlin. Accelerating eulerian fluid simulation with convolutional networks. Proceedings of Machine Learning Research, pp-
2017.
Kiwon Um, Robert Brand, Yun (Raymond) Fei, Philipp Holl, and Nils Thuerey. Solver-in-the-loop Learning from differentiable physics to interact with iterative pde-solvers. In Advances Neural Information Processing Systems 33, 2020.
Benjamin Ummenhofer, Lukas Prantl, Nils Thuerey, and Vladlen Koltun. Lagrangian fluid simu- lation with continuous convolutions. In International Conference on Learning Representations, 2019.
Renato Vacondio, Corrado Altomare, Matthieu De Leffe, Xiangyu Hu, David Le Touzé. Steven Lind, Jean-Christophe Marongiu, Salvatore Marrone, Benedict D. Rogers, and Antonio Souto-Iglesias. Grand challenges for Smoothed Particle Hydrodynamics numerical schemes. Computational Par- ticle Mechanics, 8(3):575-58) 5588, May 2021. doi: 10. 10.1007/s40571-00 -20-00034-111
Rui Wang, Robin Walters, and Rose Yu. Incorporating symmetry into deep dynamics models for improved generalization. In 9th International Conference on Learning Representations, ICLR 2021, Virtual Event, Austria, May 3-7, 2021. OpenReview.net, 2021.
Shenlong Wang, Simon Suo, Wei-Chiu Ma, Andrei Pokrovsky, and Raquel Urtasun. Deep paramet- ric continuous convolutional neural networks. In Proceedings the IEEE conference on computer vision and pattern recognition, pp. 2589-2599 22999 2018.
Chi Zhang, Massoud Rezavand, and Xiangyu Hu. A multi-resolution SPH method for structure interactions. J. Comput. Phys., 429: 10028, 2021. doi: 10.1016/f..c. 10028.
Jie Zhou, Ganqu Cui, Shengding Hu, Zhengyan Zhang, Cheng Yang, Zhiyuan Liu, Lifeng Wang, Changeheng Li, and Maosong Sun. Graph neural networks: A review of methods and applica- tions. AI Open, 1:57-81, 2020. doi: 10. (0.1016/j.aiopen.222 2021.000000...00000..
13

---

Published as a conference paper at ICLR 2024
Table Overview of the different Ablation Studies
Section
Test Case
Hyperparameter
Purpose
C.I.1
Number of Basis Terms
Symmetric learning task performance Anitsymmetric learning task performance
Number of Basis Terms
Basis Function Layout
Comparison against baseline methods
C.3.1
II
Basis Function
Comparison against baseline methods Influence of including certain sine/cosine terms Repeatability and statistical significance
II
Fourier Series Terms
Random Seed
Coordinate Mapping
Evaluation of prior inductive biases
II
Window Function
Evaluation of prior inductive biases
C.4.1
Number of Basis Terms
Validation of one-dimensional results
C.4.2
Network Layout
Search for optimal network
C.4.3
Basis Function
Comparison against baseline
C.5.1
IV
Basis Function
Comparison against baseline
C.5.2
IV
Network Architecture
Influence of Overparametrization
C.6
Computational Performance
Comparison against baseline
APPENDIX
the appendix, we will be providing additional information to the main paper. First we will be discussing additional details regarding the basis convolution model, e.g., defining all choices of basis functions we considered, see Appendix A Next, we will be discussing details regarding the simulation setup and data generation, e.g., how the random initial conditions were chosen, Ap- pendix B Finally, we will be providing additional results for all of our ablation studies, see Table for overview of the ablation studies, in Appendix C
A SUPPLEMENTARY MODEL DETAILS
Here, we will discuss the mathematical foundations of our convolutional approach and mathematical definitions of all used basis functions, window functions, and coordinate mappings.
A.1 MATHEMATICAL MODEL
In general, a convolution in ID is mathematical operation two functions / R + IR, the input function, and g R + R. the filter function, which are convolved using a convolutional operator o, ie., h(z) (f o (z). This convolutional operator in continuous form is defined as
(f (fag)(z) = [ f(x 7))((dd)
(11)
In Machine Learning problem, the goal is now that given an input function f. e.g., the input feature vector at a given set of positions, and target output function h, i.c., the ground truth, to find a filter function g such that is the convolution of f and g. To achieve this, g needs to be parametrized into learnable form ge. where we consider three primary approaches:
Using Multilayer Perceptron (MLP)
Using Radial Basis Functions, e.g., piece-wise linear functions
Using Traditional approximation techniques, e.g., Fourier series
Note that for the discussions here, will only focus the latter two approaches they very similar, and most of our evaluations concentrate on these approaches.
For Lagrangian simulation, e.g., using SPH, several inductive biases can be included to help network. As particle simulations only have finite set of particles, we can reformulate the convo- lutional operation where interactions based the edges of graph, with particles being
14

---

Published as a conference paper at ICLR 2024
the vertices and the edge features being pair- wise particle distances. Accordingly, we can rewrite the convolutional operation for n particles and apply it at the location of a single particle 2y to be
(J (o))
> ((r))go(<; xz):
(12)
which is a direct discretization. However, this formulation is not very practical as every convolution would require an interaction for every particle-pattille pairing, i.e., computational complexity of (n))). To resolve this issue, we apply an inductive bias motivated by SPH and similar methods in limiting the interactions being the support radius of the particles. Note that for convenience, we will assume h = for further discussions (which can be ensured through appropriate scaling of the positions z. This results in change in the convolutional operation as
compact domain, i.e., ge is zero outside of the interval hh, with
(f o (z)) = E f(zz)go(a; - xz), JEN,
(13)
where Mr; all the neighboring particles of i with la; xzg S h. Based on this convolution, some approaches introduce a normalization function ((x) as
(f a ge))
Σ f(x)) -x,): JEN;
(14)
<(x)
However, an inductive bias we apply is that the filter function g should behave similarly to an SPH interpolation, i.e., fewer neighbors lead smaller interpolation values. While this, on a numerical level, leads to a kernel deficiency problem in SPH for free surfaces, it is nevertheless widely used SPH technique. This can be achieved simply by setting (xx) 1, which has been done in prior work, e.g., in the LinCConv approach (Ummenhofer et al. 2019).
further inductive bias we apply to exclude the particle from interactions with itself and instead use fully connected layer for self-interactio 2022) and is primarily motivated by gradient interpolations SPH, which include no contribution from particle on itself, similar to central difference schemes. If learning task were to learn such a gradient interpolation, either explicitly or implicitly, as part of more involved learning problem, any weight describing the self inffrinteraction would have be zero. This, however, overly restricts the weights, and is more straightforward to exclude such an interaction, e.g., as
This bias has been utilized before (Prantl et al.)
(f o (z;
Σ f(z,)9o(1; - 2,) wf(xs). JEN
(15)
We will exclude this bias from further discussions for convenience and readability.
To parameterize ye, we define ga as a combination of ri basis terms (T()), with an associated weight 0, where the overall filter function is summation of these terms. For example, for the LinCCom approach, the basis terms b, describe a linear interpolation over the domain [-1, 1] with each basis function being piece-wise linear. For computational efficiency, as a dot product of a basis term vector b with the weight vector θ as
reformulate this summation
ge (x) (b(x), O).
(16)
We impose further inductive bias for two-dimensional convolution in that the filter functions are separable w.r.t. the coordinate axes, e.g., the LinCCom approach uses a bi-linear interpolation with piece-wise linear separable basis functions. Treating as a weight matrix with basis functions b,, with u components, and by, with v components, and x = Pyl, we can re-write the convolution:
(x) = ΣΣ B(r),, Ougi B(x) = by 8 (y
(17)
Which works analogously for three dimensions. At this point an important observation is that if b and by are orthogonal bases then the resulting two dimensional basis will also be orthogonal. To prove this first consider the definition of orthogonality for basis in dimension, which is
15

---

Published as a conference paper at ICLR 2024
defined based the existence of inner product such that (f.g) = 0. which is evaluated through an integral of the form
(.9) J (())(()))(()))
(18)
where w(x) is weight function. Given some functions, e.g., an orthogonal polynomial basis, the corresponding series can be written (x) - 42 azF3 (2)), where F is the basis polynomial and aj is a sequence of coefficients. For orthogonality to hold for this series any inner product of two basis terms needs to either non zero, if the same term multiplied with itself, or otherwise. Accordingly, we can write the orthogonality constraint as
S (()))(())))))) = O, Vi,jENAi
(19)
For example given the Chebyshev polynomials of the first kind (x) defined through the recursion
Tu(x)= = 1,
(xxx) I, (1) = (x) + (n))),),
with an according orthogonality constraint of
S. [(=)) =0, O, Wi.jNNiijj.
(20)
As an example, consider (xxx) and (2))), which yields
[ 18x4-8.2 8x² +
1
[Ax*
3r]
dr
(21)
x²
which can be readily integrated and yields 0 as its definite integral. The orthogonality constraint in two dimensional function is defined as an inner product over a square region and can be written as
Sf [(x
v)g(s, (a,,v0oo w) y)drdy O,
(22)
with a weight function w dependent on both parameters. For our basis convolution approach we now consider two, not necessarily identical, orthogonal bases a and b defined via outer product as
B(I, v) (. (y()), Vi
(23)
which means that for orthogonality hold any inner product of two basis term Bi, (I. (r.y) = a, (x())
(2) and (x, y) = a4 (z) (y() is either non zero, for i & k A & I, or 0. To show that this is true we first consider the orthogonality of the basis terms Brg along each cardinal direction, i.e., we want to verify that
[i fid (2, (x, (2)(x)dx=
[" (1)),,,))) (z)))) v)dz 0, Vy E I, (24)
and analogously along the other axes. Considering the integrand, b, and by are independent w.r.t. the variable of integration and can be moved out of the integration to yield
[ (z)bb(())) an (z) (b) br())] S" ax(x) (x(x))((, y)dz.
(25)
The right hand part of this equation is the same the orthogonality requirement for itself and, accordingly, the integral is zero if w(x, y) is identical to the weight function w, (z) for the orthogo- nality of a. The derivation along the other axes proceeds analogously. We consider the original
16

---

Published as a conference paper at ICLR 2024
orthogonality requirement again and refactor w(x, y) wa(x) mo()), which yields:
lf Ja w)... ())( yddddy,
/ʃ Jig ()n,,,)),n,, v) [wa (x) wa(v)) dady,
= L [ (,()),,))) ax (z) (ua(x) ws(y)) drdy
(26)
[ (a:(x) as (x) (x)) (v) bi())ws(y) drdy
S (x(x) a ())) [ bv) by (v)
Thus we have shown that given two orthogonal bases a and 6, with respective weight functions Wa and wp, the two dimensional basis resulting from an outer product is orthogonal with respect to w(x, y) = Wa (2(z) ((xx)//(())..)..))) Accordingly, using any basis along an axis and using this seperable construction results in orthogonal combined basis. Note that this can be analogously shown for three-dimensional bases that are the tensor product of three bases.
For many learning problems, is not a scalar function but describes a multi-dimennional input function f IR + Rimpur, with input features defined at all locations. Furthermore, the output of a convolution multi-dimensional i.e., f fog E Roupput where we impose the bias that each input feature should be connected with each output feature, i.c., the convolution should be fully connected, which is in line with prior work. Accordingly, we can redefine the convolutional operator as
inputFeetures
(f a (1) Joupout
Σ
Σ
g
input, (x..
J (Iz).
(27)
JENA
input=1
which can be reformulated by summarizing all relevant input feature associations for one output as jagufFoatuces.coupu... (x)
ourpat
(x)
[g
1output
(x),
88
(28)
(f o ge)
Σ
(gg
oumpur
(x: 22),f(r;)).
(29)
JENA
We now compute the messages M, which are an n X o tensor (with o being the number of output features and n the number of edges of the graph), based on the basis function tensors Bx and B of shape n x u and n x v, respectively, as well as the input feature tensor F, of shape n x (with input features), as well as the weight matrix W, of shape X v x x o, using Einsum notation as
M = B B W Wuvio FRS.
(30)
which can be efficiently implemented using either the built-in einsum function of PyTorch (Paszke et al. 2019) (and related frameworks) brary (Thakkar et al., 2023), as done by Symmenhoffer et al, (2019). Computing the gradients of such an operation is, mathematically, straightforward it is just sequence of matrix multiplications, and the shape of the actual base functions comprising B do not affect the shape of the gradients.
more direct implementation such as Nvidia's Curlass li-
However, relying on auto-diff gradients can be impractical for such operations as the intermediate matrices can be large and must be stored for every convolution operation in a Neural Network. We chose to implement this process through custom forward and backward operation that does not compute B for all edges at once. Instead it performs this operation in batches of size b, which limits the memory requirements significantly as there is no need for large intermediate matrix. Further- more, we recompute the basis terms 8 and B during backpropagation as this allows us to only store the particle distances, which are shared for all layers, and input features, per layer, for back- propagation, instead of large matrices that are potentially different per layer. Whilst this imposes some computational overhead, it can also significantly reduce memory requirements, allowing train- ing on GPUs with 4GByte of VRAM and less, for three dimensional networks, with the batch size parameter b trading off computational performance and memory requirements. For details on computational requirements Appendix C.6
17

---

Published as a conference paper at ICLR 2024
A.2 WINDOW FUNCTIONS
A further inductive bias applied by prior work is including so-called window functions in the convo- lution operator. These window functions are inspired by the shape of SPH kernel functions, which zero at the support radius and tend to be shaped like Gaussian functions. Imposing this bias can be achieved straightforwardly by including an additional term in the convolutional operation as
X
(f gx)
Σ JENS
go(x: xx))((x)) w
(31)
(xx))
where W is the window function, which is zero at and Gaussian-shaped. While many choices exist for these window functions, we limit evaluations to the following window functions:
Using no window function is the most straightforward choice as this can be implemented by simply not including the window function. This term could also be defined as
None:
11.
1 0, else.
Wone (r)
(32)
Linear: A naive choice for window function is a function that is 1 at the origin and linearly decays towards 0 at = 1. While this function not very Gaussian shaped, it does not significantly impact the shape of the learned convolutional operation besides tending towards and can be defined as
Linema (r) = [1 r)+
(33)
where F1 = max (-:0).
Parabolic: A slightly more complex variant of the prior Linear window function uses a parabolic decay instead of a linear one. This window function can thus be defined simply
W Paraholic (r) = [[1 r
(34)
Müller: This window function is based on Müller et al. (2003) and defined as a polynomial series of order 6. The advantage of this window function is that is Gaussian-shaped but does not require any square root operations as the distance is only used squared and is defined as
Wrllller (r) [1
(35)
Spiky:
This window function is based on Müller al. (2003) and is purposefully designed to be Gaussian shaped, i.e., the gradient of the kernel function does not tend towards 0 as r tends to This is imposed ensure that particles keep repulsing each other in SPH, which avoids particles clumping up unnaturally, also known as the pairing problem in SPH. This function is defined as
Spiky (r) = [1 r
(36)
Piece-wise Bezier splines are a popular choice in SPH for kernel functions and find wide usage, especially within Computer Graphics focused research (Thmsen et al. 2013). For these functions, three popular choices exist based on the degree of the spline, i.c., cubic, quartic, and quintic splines (Dehnen & Aly 2012). These are defined, respectively, as:
B-Splines:
1
3
Cubie (r) = [1 r)
4
(37)
2
3
WQQartie (r) = r|
10
(38)
5
5
2
5
1
5
M M-Quitice (r) = r
15
(39)
3
3
18

---

Published as a conference paper at ICLR 2024
A.3 COORDINATE MAPPINGS
So far, we only considered the coordinates to be given as Cartesian; however, it might be helpful to utilize other coordinate systems. Coordinate mappings A applied to the filter function ge
(f (x.)
Σ
(xxxxx))((xx))) w
[x - Xj
((x+)
h
JENM
The mapping function A thus is a function A : Rd - Rd, which could be defined in arbitrary ways, e.g., in one dimension used coordinate mappings. Accordingly, no functional mapping exists in one dimension besides an identity mapping, i.e., using A(x) x. For two-dimensions (which also can be expanded similarly to three-dimensions), we consider the following three mappings:
could utilize (xx) = =2? however, we only consider commonly
Identity: This mapping serves as the baseline approach of directly using the Cartesian
A Identity
(x) =
As the imputs to A are distances, limited by spherical support radius h, direct choice for coordinate mapping is to map the input to polar coordinates. Note that this implies some necessary changes to how the basis functions are evaluated; however, we will skip the details here for brevity as the results indicate no significant gain in using this mapping. This polar coordinate mapping can be defined straightforwardly using the atan2 function as
Polar:
= 2 |2|x||2 1, 1 atan2 (x Xx)
A Polar
(x)
where the scaling ensures that the domain remains unchanged, i.e., A (-1,111
Preserving: The preserving mapping, proposed by Ummenhofer et al. (2019) and based on Griepen- trog et al. (2008), is intended to remap the spherical support volume of the convolutions to a cubic domain ensure that each weight influences comparable amount of space. While this becomes more important in three dimensions, as the volume ratio between cube and sphere is much worse than that of square and circle, can still be applied two dimensions by setting the component when performing the mapping to zero. This preserving mapping works in two-stage process where a ball first mapped to a cylinder and then mapped to a cube. This mapping is defined as:
(0,0.0),
if[|a[12
alla V/ (T=.0) la 32)
Ioll
ir S x² x²
Aball bcll--cy((q))
(40)
FV [al
3/lll2
/y Talla+1a sgn(z))
3|alle
else
(0,0,)),
ifr = 0, =
Acst-rcubeq))
(sgn(x)]l<z, Wllk, # sgn(x)]((2,, v) arctan _. (gn())]l<r,))) v)1 aretan " sgn(v))l<z,))ll v)ll,.2)
iflyl S |x]
(41)
else.
A.4 BASE FUNCTIONS
So far, we only considered an arbitrary basis tensor (x) and now we would like to discuss all the utilized base functions within this paper with particular emphasis basis. For the versions based on Radial Basis Functions (RBFs) we assume a Cartesian coordinate system with an evenly spaced grid of central points c; among the æ and y axes, as x; -1+
our proposed Fourier
and = 1 + - respectively with a separation distance of Ar. Due to our assumption of scparability, the functions are defined equally regardless of which axis they are applied to, with the sole exception being the DMCF formulation, which requires some further processing. Furthermore, for each of these central points, we can compute the relative signed distance r; from the input point and the centroid, i.e.,
p C, Accordingly, each RBF is defined as b; (q 2
Finally, most of these basis functions are designed to act as partitions of unity, i.c., with a weight vector θ 1, the resulting filter function şe is everywhere. This is essential, want these
19

---

Published as a conference paper at ICLR 2024
Figure 10: The 5 x 5 basis function combinations for the DMCF basis (left) and SplineConv (right).
methods to act interpolation functions rather than approximations. Enforcing this property is possible by either normalizing the output of g. i.c., by introducing a corrective term > ba "
or by adjusting the definitions of the basis functions such that they do not require this corrective term. We chose the latter option the former introduces modifications the shapes of the basis functions, e.g., the basis functions for the corner terms might appear different than the central points, which is undesirable, and this approach has been used by prior work well (Fey et al. 2018).
Nearest Neighbor: Nearest Neighbor interpolation is simple baseline to compare all other meth- ods against and is constructed simply note here, however, that the function needs to be carefully designed such that there is no overlap of the basis terms on the edges of their influence radii, i.e., using [g] S 0.5 would lead multiple bases contributing to the same input p, which would violate the requirement of a partition of unity. This basis term can then be defined as
group of piece-wise constant functions. An important
[1, 0. else.
<
(g)
(42)
LinCConv: A higher-order interpolation scheme is linear interpolation, also used by (Ummen- hofer et al. 2019) and referred to by this name in this paper. Building a linear interpolation as a radial basis is straightforward using piece-wise linear definition, which be further simplified by using the [+], notation used before as
BlinCCom (g) = - 141l+
(43)
In prior work, modification of the linear basis was proposed that includes only antisym- metric terms, primarily focused on antisymmetries of the combined coordinate-axes Consequently, these terms have few symmetries along the individual axes, Figure TO but always antisym- metric overall. |Prantl et al. (2022) defined these terms as a bilinear basis and then modified the filter weights separately after each weight update be antisymmetric. However, this results in the network not seeing the correct gradients and losing half of the weights to enforce antisymmetry. Instead, we define the DMCF basis using the LinCConv basis with modification
DMCF:
(44) which utilizes all weights but increases the interpolation frequency in by factor of 2. For com- pleteness, an alternative formulation that also ignores half the weights but still yields correct gradi- ents without increasing the interpolation frequency could be formulated as
UMCF (q) ([sgn(q.) LinC'Conv (2/ 1) @ b LinCConn (sgn(qx) {gy)] 1),
DMCF2 q) gla
[[sgn(q. LinC'Con, ((q.)) @ LinCConn ,/ (sgn q.)] .)
(45)
B-Spline/SplineConv A natural extension of linear interpolation is using higher-order spline func- tions, where the cubic B-Spline function has been used before by Fey et al. (2018). However, ensur- ing that these methods are partitions of unity is more challenging requires adjusting the spacing
20

---

Published as a conference paper at ICLR 2024
of the centroids by scaling them by 1 i..., the new centroids are ci &plinee = (1- n) (-1+ n. Furthermore, the width of each basis function needs to be adjusted, compared to their definition as window functions used before. We evaluate these modified widths through an optimization pro- cess, i.e., we optimized this parameter such that the interpolation results in partition of unity with minimal width per basis. This results in three B-S B-Spline basis terms:
3
q
6 SplineCom (4)
=
1
4
(46)
2
C
QQuartic (q) = 1 l
3
1
9
5
+ 10
(47)
5
5
lg]
2
1
1
Quintic (4)
1
6
15
(48)
3
with normalization constants e of 1.732051, 1.936492 and 2.121321
(Dehnen & Aly 2012).
An alternative kernel often used within CFD CFD-oriented SPH applications is part of the Wendland series of kernel functions (Dehnen & Aly 2012 Sun et al. 2018). These functions are also polynomial but do not exhibit some of the numerical disadvantages as the B-Spline kernels and, as such, are an interesting alternative to evaluate for a machine learning context. These terms also require the modified spacing, identical to the B-Splines, and the Wendland-2 basis is defined
Wendland-2:
Wendland (a (()) (1)=111 c
-11-10 |u| ] (1 (1 + |0| c
(+) ; 1.620185.
(49)
Gaussian: A natural extension of the Spline bases is utilizing non non-compact Gaussian functions, i.e., using exponential functions. As these are locally defined, i.e., only as part of the filter function and not used to determine the graph edges, non non-compact functions do not impose any of the usual drawbacks. These terms also require the same spline centroid and are defined as
Gaussian (q) exp
(50)
Spiky: The Spiky kernel, discussed before, is primarily included as additional case for ablation
studies as, due to the shape of the function, this basis term cannot be normalized by adjusting the spacing and width of the basis. Nevertheless, this basis term is defined as (Müller et al., 2003)
b 8Spkk (q) [ - 14
(51)
RBF Bump: There is rich history of Radial Basis Functions within RBF Interpolation theory, and we chose one of these terms that is different from the classic bases. Note that for an actual RBF basis, and thus RBF network, the centroids would need to be learned as well, which is non-linear optimization task, to achieve proper interpolation qualities and not doing so, as done here, is unlikely to work but serves as a valuable baseline for ablation studies. Nevertheless, this basis is defined
exp 1 (1-00.1873961119444666666
-1
ifr else.
J Bump (q)
(52)
0.
In addition to this set of radial hase terms, we also consider two traditional approximation tech- niques that operate more globally, i.c., they are directly evaluated on p instead of g. While many approximation bases exist, we primarily focus on Fourier series terms and Chebyshev polynomials as they find wide application for interpolating periodic and non-periodic functions.
Fourier: This method is the primary focus of our paper and results from using Fourier-series as the basis. The first term of a Fourier series is always constant function, i... bp(p) = 1; however, for higher order terms, we could either first use the sine or cosine term. We will later include an ablation study of some variations of this ordering, but the following is the standard definition:
1,
i= 0,
Fourier (p)
COS π
|
+1pp)
odd,
(53)
sin
[1
+1 p). i even.
21

---

Published as a conference paper at ICLR 2024
Figure 11: The 5 x 5 basis function combinations for Fourier basis (left) and Chebyshev (right).
For some of the variants we evaluate in our ablation studies, we modify the first term, i.e., by F as this term is inherently symmetric, but want to evaluate purely asymmetric Fourier terms as well. To achieve this, we either (a) drop the term, (b) modify the term to be based on sign of x, i.e. bp
Fourier, sgn (x) = (g) or (c) use directly, i.e., by Fourier. linear (x) = x. Using the standard basis definition for = 5 in two dimensions,
Figure [ w.r.. [r, w/. in a variety of configurations. the r and coordinates as well as the combined coordinates
find symmetries and antisymmetriess,
Another popular choice in graph convolutions are Chebyshev basis functions, which smooth, inherently symmetric, and antisymmetric. For these terms, we primarily consider Cheby- shev polynomials of the first kind; see Figure II for visualization of this basis, defined as:
Chebyshev:
to
Cbebysbev
(p)
b1 Chebyshov (p)
(54)
by
Chebyyshv (p) 2pb,-1
2hcchchysher (p)
D
Chehyssee (p),
where the magnitude of the basis terms is bound by the domain [-1, Furthermore, for some ablation studies, we also considered Chebyshev polynomials of the second kind, defined as
by Chebyshev2
(p) =
b Chebyshev2 (p) = 2p
(55)
by Chebyshev2
(p) = 2pb Chebyshev2 (p)
Chebyshov2 (p),
6-2
which are not bound in magnitude by the domain [-1,
B EXPERIMENTAL DETAILS
We focus on several datasets with focus on quantifiable behavior, and describe in following how these datasets were generated, Existing datasets oftentimes involve behavior that is difficult to numerically quantify or behavior where loss metric difficult relate to the visual perception of the simulation.
We chose create our datasets to include a wide variety of SPH problems across one, two, and three dimensions to make them as versatile as possible. In this section, we will focus on the setup of the solvers for the generation of our datasets, as well as basic network parameters and data augmentation techniques. Accordingly,
will be discussing each test case in separate sub-section.
Overall, there are few similarities between our different test cases and datasets; however, some fa- miliarities exist. Notably, all of our datasets, well as the implementations of our used classical solvers, and the implementation Furthermore, https///gthhbb.ccm/uumppbboo
our network architecture is available as open source code at
utilized the Adam (Kingma & Ba 2015) op- timizer in all cases. We built all of our code, including the SPH simulations, using Py Torch and PyTorch Geometric for graph processing, e.g., neighborhood searches.
22

---

Published as a conference paper at ICLR 2024
590
500
1000
1000
1500
-11
1500
1.00
590
500
14
1000
1000
10
1500
1500
0.25
1.00
Visualization of the test dataset for test case consisting of simulations of 2048 timesteps each. Each plot represents a simulation from the dataset, with simulation time on the (from 0 to 2.048s) axis and partiele index the y axis (from to 2048), with color mapping indicating particle velocity (blue for up and red for down). The initial conditions were randomly generated with identical simulation parameters, c.g., viscosity, across all simulations and not spe- cially selected for good or bad testing performance for some methods.
Figure 12:
B. TEST-CASE ONE-DIMENSIONAL COMPRESSIBLE SPH
Lagrangian fluid simulations pose several problems for machine learning techniques compared to Eulerian simulations. One primary consideration is the particle spacing and distribution, which, for Lagrangian simulations, changes the flow evolves. This means that if the predictions of network lead to errors on one inference step, then the input to the next inference step is out of distribution, compared to the ground truth. Accordingly, changing and varied particle distributions and that the network solution does not learn behavior that only useful for very narrow set of distributions. Consequently, data augmentation techniques are essential during training: however, in this test case, want to specifically investigate how different particle distributions affect the network performance.
vital that machine learning techniques can deal with
For incompressible simulation, the spacing between particles is somewhat constrained due incompressibility limitations; however, for a compressible simulation, we can generate particle dis- tributions that cover much broader range of inter-particle spacings. Furthermore, by reducing the dimensionality of the simulation, we can perform much more focused investigation of this relation- ship. Accordingly, our case is a compressible one-dimensional SPH simulation that primarily investigates differences in methods' capabilities to handle broad range of particle distributions.
Our underlying simulation model, in this case, is simple Equation of State (EoS) based explicit time integrator using Runge-Ke Runge-Kutta Kutta integrator of fourth order (Antuono et al. 2012) combined with an explicit diffusion term for the velocity field. As the EoS, we chose an ideal gas equation with no influence of energy temperature, as the general NNTI approach only predicts changes in position and not in energy. For the diffusion term, we chose the approach of PPice (2012), although the exact choice of diffusion term does not matter for our evaluations. Finally, we chose periodic domain and simulation region of Ω = [-1, 1] as boundary conditions.
To set up the initial conditions, we use a two-step process where we first create random density profile for a fixed domain of [-1, and then place a set number particles (2048) such that their summation density P; = Σ, m, iy is identical to the random density profile. To produce the initial density profile, we utilize Perlin noise (Perlin 1985), perceptually isotropic gradient noise commonly used for procedural generation. Perlin noise, in general, is based on pseudo-randon process based seed that generates n-dimensional noise texture of frequency along each dimension, defined as Py I + [-1, 1]. A common technique is to combine multiple Perlin noise textures of increasing frequency generate so-called Octave noise using several octaves Moocavve
23

---

Published as a conference paper at ICLR 2024
32
1000
90
110
130
2000
31
2000
-22
1000
1000
400
1.1
2000
32
B-22
-22
2000
-22
-01
1:10
12
-0.00
-12
0.10
23
1.00
2.
-01
-11
2000
-3
35
0.
30
-22
-01
2000
0.0
0.5
1.0
13
2.0
0.0
0.5
1.0
1.5
2.0
0.0
0.5
1.0
2.0 0.0
0.5
1.0
1.5
2.0
Figure 13: Visualization of the training dataset for test case consisting of 32 simulations of 2048 timesteps each. Each plot represents simulation from the dataset, with simulation time the (from to 2.048s) axis and particle index on the y axis (from to 2048), with color mapping indicating particle velocity (blue for up and red for down). The initial conditions were randomly generated with identical simulation parameters, e.g., viscosity, across all simulations.
a lacunarity / denoting the increase in frequency per Octave and a mixing factor a, as
Postae
Octave M a' (x).
(56)
To generate our dataset, set Flceaae 4, and As the noise is the range [-1, 1], we need to modify the noise amplitude by scaling the noise by 1 and adding an offset of 2, meaning the density field is pscudo-random periodic density field in the domain [1 3. 2 1]-
While it is possible use this octave noise the density field directly, this would lead to disconti- nuities at the boundaries. One technique to sample periodic noise embedding a regular structure into higher dimensional noise function. By embedding circle, which has constant curvature, into a two-dimensional Perlin-noiss function and then using parametric description of circle, i.c., [r cos sin 0) for 0 [- R), we can generate periodic one-dimensional noise.
24

---

Published as a conference paper at ICLR 2024
This density field is then sampled on discrete grid, with 2048 grid points, yielding a discrete Prob- ability Density Function (PDF). We then compute a discretized Cumulative Density Function (CDF) by integrating the discrete PDF. This discrete CDF can then be used to construct approximate in- verse CDF by building a linear interpolation using the discrete CDF as z-coordinates and a regularly spaced range from [-,, as the y-coordinates. We then regularly sample the domain [0, 1] using pernilles points and sample the inverse CDF using these locations to find the initial particle locations.
After sampling the particles, we initialize the simulation with a constant initial velocity of 0, particle size a "punidss a support radius / = 4 * a, diffusion coefficients a = and / = 2. a stiffness coefficient = 10, a numerical speed of sound c₂ = 10, particle rest density po 1000 and with a fixed timestep of At = 10-3 We utilize 36 random initial conditions to generate our dataset and evaluate 2048 timesteps each. Out of these, we use 32 (chosen randomly) the training set, see Fig. 13 and the remaining 4 as the testing set, see Fig. 12
For the training task, we compute multi-step update, i.e., for each timepoint in the simulation, we compute the average velocity oftior over the prior s 16 timesteps and the next 16 timesteps eyea The network's ' goal is to predict Dyea based on u/ prior and the particle area, which are the only quantities that influence the partiele behavior in the underlying simulation. During training, compute the L2 difference between ground truth and network prediction for a batch size of 4 without temporal unrolling, where each batch is the result of picking 4 random samples across the entire training dataset where each training timestep is used at-most-once. Our training consists of 5 epochs, each consisting of 1000 weight updates, with an initial learning rate of 10-3 that is halved after every epoch. To evaluate the test error, we utilize all four test simulations and evaluate the Ly difference between ground truth and prediction for the time points t = [0- s, 128 + s, 256 + s, 1024 + s]
For the simulation and training on this dataset, we utilized a naïve neighbor search that computes the distance of all particles to all particles, with a computational complexity of O(n²) and calculated the distance of particles using floating point modulo operations to implement the distance checks. Due to this implementation, the positions of particles may lie outside of the domain [-1, 1]. and a modulo operation needs to be applied to find the actual position of particle within the simulation domain. Furthermore, this means that the number particles in the simulation stays constant, and no ghost particles utilized to implement the periodic boundary condition. Accordingly, the position is consistent, i.c., particles are never mapped into or out of existence based on where they are in the simulation, and their position is the direct result of the integrated velocity.
B.2 TOY PROBLEMS
SPH simulations are generally built on interpolation operations using either the kernel function or the gradient of the kernel function. Accordingly, if neural network is supposed learn entire simulation computed using SPH, then the network should be capable of learning interpolation and gradient interpolation tasks. Conversely, if network cannot learn either operation, the likelihood of the network learning an entire simulation step is low. Consequently, we devise simple set of two toy problems to evaluate the network's ability to learn these tasks.
As the basis of these evaluations, we utilize the first test case, i.e., the -ieedimmnsional compressible simulation data set, with modified learning task. In general, an SPH interpolation is defined as
(A.) = Σ
A my
Waj Pj
(57)
JEN;
with A being some quantity, mi being the mass of a particle, p being the density of a particle computed using SPH interpolation, and Waj = W(x, 24, h) is Gaussian-like kernel func- tion (Monaghan| 2005). The most straightforward SPH interpolation, then, is a density estimate that can be computed by setting A = p. i.e., (Koschier al. 2019
Pi Σ my Waj,
(58)
JEN
which can be modified to computing the number density d, Po
P
Solenthaler & Pajarola 2008
6,
Σ a,
a,Waj.
(59)
JEM
25

---

Published as a conference paper at ICLR 2024
Pe being the area of particle. This number density can be interpreted a quantity that depends single pessaag-passingg step using cutoff radius of h, vertex feature a, and symmetric convolutional filter function W. Note that in SPH, the magnitude of the kernel function scales inversely with the particle scale. Accordingly, is convenient to use a vertex feature of 1 and treat a aj as a single term. The first learning task is to learn a convolutional filter function fe as
with
Ea,Wos Σ 1 ge(zz Zz, h). JEN;
(60)
JEN
In addition to the SPH interpolation, we also consider an SPH gradient interpolation, which, in its most naive formulation, can be defined simply as (Koschier et al. 2019)
(v.A.))
Σ
A m
V.W.
(61)
Pj
JEN.
where Vi denotes the spatial derivative with respect the position of particle i. Note that many more advanced formulations for gradient interpolations are commonly used in SPH, which, e.g., are accurate for constant fields (Price 2012), but for our toy problem, a naive formulation suffices. We then derive a learning target analogous to the previous one as
v Σ V. JEN.
(62)
This task is, again, a single message-passing step; however, this time, the convolutional filter func- tion be learned has to have odd symmetry instead of even symmetry for the kernel interpolation.
For both toy problems, we evaluate the performance of a single message-passing step continuous convolutional network with no activation functions being used, as we want to investigate how the basis functions can work as interpolators for the SPH interpolations. In addition to the CConv approach, we consider more generic GNNs where we still utilize only a single message-passing step. Otherwise, the training setup is identical to the one in test case I.
B.3 TEST-CASE II: IWO-DIMENSIONAL WEAKLY-COMPRESTE
SPH
The focus of our second test case is quantifiability within a CFD-oriented context. Existing test cases often focus the visual appearance, e.g., collisions of random blobs with challenging to-quantify behavior. In contrast, our test case setup can be adjusted and has clearly quantifiable behavior.
We use closed wwo-dimensional domain Ω [_1 random initial velocity field to set up our test case. However, using a random per- per-particle velocity would not be very useful as this would lead to arbitrary amount of divergence in the first timestep that can readily lead to problems for numerical solvers. Accordingly, the valloitvelocity field needs to be divergence-free at initialization; however, this is easier said than done.
1J2 with no-slip boundary condition and
The initial particle configuration for all of our simulations is a regular spaced grid of 64 x 64 particles such that the summation density of all particles at
is equal the rest density. This sampling means that there exists a regular grid with cell centers identical to the particle positions, which makes the initialization more straightforward The initial velocity field is then computed using curl noise (Bridson et al. 2007) sampled on the grid and then resampled back to the particle set.
The idea behind curl noise is that the gradient of the curl of random potential field can be used to create a velocity field that, by definition, is divergence- free. The initial potential field is sampled based on a two-dimensional Octave noise field, using Perlin-noise, as this smooth, relatively low- frequency noise leads to large structures that can be easily simulated using direct numerical solvers. However, this potential field would still create issues in the actual simulation as it does not respect the boundary conditions. As we utilize no-slip boundary condition, the most straightforward way is to ensure that the potential field is constant at the boundary, e.g., O, such that the curl is also 0. Based on the suggestions by Bridson et al. (2007), this can be ensured by linearly blending the potential field to be 0 at the boundary and using a gradient along the boundary normals to blend the potential field from its actual value with 0 based on the boundary distance.
The next step is to compute the curl of the potential field, which is done on the grid representation of the simulation domain using first-order central difference scheme. This curl is then directly
26

---

Published as a conference paper at ICLR 2024
Visualization of the initial conditions of the 36 initial conditions used for the training dataset for test case II. The simulations here are run using weakly compressible SPH model with an initially divergence-free velocity field with solid boundary conditions outside the shown areas. Color mapping indicates velocity magnitude, and streamlines indicate the flowfield mapping the particle velocities to regularly sampled grid using SPH interpolants.
Figure 14:
0 computed by
mapped to the particles, where we utilize an SPH gradient interpolation to evaluate the initial veloc- ities of all particles. Note that it is crucial to compute the gradient using SPH gradient operator as, otherwise. the divergence of the velocity field is not quite zero.
We can then generate our dataset. To make the cases more interesting, quency potential field, i.e., the test simulations have larger structures that are not seen as such during training. Overall, generate 32 training samples, Fig. 14 and testing samples, Fig. 13
lower fre-
The simulations are then performed using a 8-SPH method (Marrone et al. 2011) with a Runge- Kutta time integration scheme of fourth order (Antuono al.) (d.2)) Each simulation contains 2048 timesteps where the learning task is to compute the position of particle for once. Note that the underlying 8-SPH simulation uses a so-called continuum density formulation, i.e. the density of any given particle is not the result of an SPH interpolation, as in test case but is integrated over time using the continuity equation. However, we only provide the velocity of particles features, making the learning task ambiguous and more complex, thus making it more interesting. Furthermore, due to a combination of the multi-step prediction task, the RK-4 time integration, and the particle scaling, the receptive field of the particles in the simulation covers the entire domain, i.c., each particle influences other particles; however, we utilize the architecture of Ummenhofer al. (2019), which only uses 4 rounds of message-passe passing. As such, the problem is
timesteps
27

---

Published as a conference paper at ICLR 2024
Figure 15: Visualization of the initial conditions of the initial conditions used for the test dataset for test case II. The simulations use weakly compressible SPH model with an initially divergence- free velocity field with solid boundary conditions outside the shown areas. Color mapping indicates velocity magnitude, and streamlines indicate the flowfield at t = 0 computed by mapping the particle velocities to a regularly sampled grid using SPH interpolants. Note that the frequencies are lower than the training samples to better evaluate generalization.
made even more difficult as the networks are not provided with all seemingly necessary information. All of these limitations make the learning task here exceedingly difficult the network has (a) all necessary features, (b) receptive field that is too small, and (c) the testing samples are out out-of-bar ba-d.
A further complication arises from the nature of the simulation. Vortices SPH can be challenging to simulate as they often lead to the formation of holes, i.e., regions of the simulation domain with no particles. While we avoid this during the simulation through particle shifting (Rastelli et al.) (022). the network must also learn this correction on top of the physics update. However, this is made more challenging as the dataset only considers already corrected particle positions, i.e., the network never sees correction explicitly and, thus, has to implicitly learn this behavior.
We compute various metries to evaluate this dataset and focus on metrics calculated on regular grid that spans the closed simulation domain. This regular grid data can be computed by performing an SPH interpolation for all grid centers, e.g., for velocity and density. The regularized data has the advantage of not overly punishing particle swaps, i.c., two particles may be located at each other's location in relation to the ground truth, leading to a significant per-particle error even though the underlying velocity field is accurately predicted. We compute the L2 of the density, velocity, and divergence field on the regularized data, where the latter calculated using an SPH gradient interpolation on the cell centers. Finally, we also compute a Power Spectrum Density (PSD) on the grid velocity field to compare the overall structures, irrespective of their spatial location.
For training, utilize 20 epochs, each consisting of 1000 weight updates, with an initial learning rate of 10-3 that is halved after every five epochs. We start the training with an initial rollout length of I, which increased every second epoch by 1, up of 10. Evaluations are computed on all test samples from the dataset for frames (s. 512 + 8, 1024 + s, 21715 + s] for an inference length of 64, where is the prediction distance chosen as 16. During training, we use batch size of 2 and augment all samples by including random jitter for the particle positions, set to be normally distributed with a standard deviation of 0.01h, and a random rotation (Brandstetter et al. (2022a). Note that we do not apply noise during the unroll while training and only use the augmentation on the first step, similar to Prantl et al. 2022).
maximum unroll length during training
B.4 TEST-CASE III: TWO-DIMENSIONAL INCOMPRESSIBLE SPH
Our other test cases far consisted only of particle distributions where particle neighborhoods were always fully occupied. However, one of the significant advantages of SPH-based simulations is the ability to handle deforming free surfaces where particle neighborhoods are not fully occupied. Ac- cordingly, it is essential that methods also applicable in simulations involving free surfaces. A common task for this problem is the collision of two or more blobs liquid free space, and while these simulations can be visually impressive, they are often not very challenging in their dy- namics. Our dataset aims to be a more challenging variation of such test cases by increasing the difficulty of the physical setup and using more accurate solver.
Generating the initial conditions in this test case is much more straightforward than the prior two test cases only generate spheres of liquid, with larger one located at the origin and smaller
28

---

Published as a conference paper at ICLR 2024
Figure 16: Visualization of the training dataset for test III showing all 64 constellations. Each initial condition chosen by randomly placing a smaller sphere around a larger sphere with random offsets, causing generally similar but different behavior. Color mapping indicates particle velocity from low (blue) to high (yellow) after 32 timesteps.
Figure 17: Visualization of the testing dataset for test case III showing all 4 constellations. Each initial condition chosen by randomly placing a smaller sphere around larger sphere with random offsets, causing generally similar but different behavior. Color mapping indicates particle velocity from low (blue) high (yellow) after 32 timesteps. Training samples are chosen in-band, relative to the training samples, but with unseen rotations and offsets.
randomly placed orbiting the larger We then give the orbiting sphere initial velocity towards the origin and apply an offset in an orthogonal direction the velocity to offset the impact location generate more variety. We then move the smaller sphere to collide with the larger sphere
29

---

Published as a conference paper at ICLR 2024
after a small number of timesteps, where this number of timesteps is consistent across all samples. Based on this setup, we now generate 68 simulations, with 4 random simulations, see Fig. 17 used as in-band testing samples and the remaining 64 used for training, see Fig. 16
As the underlying solver we utilize a divergence-free SPH solver set to 0.001% maximum density error and 0.01% maximum divergence error, resulting in, average, 64 substeps per simulation timestep. We base this solver on a summation density approach and use a simple kinematic viscosity formulation with low viscosity coefficient stabilize the simulation. Time integration is done using an RK-4 scheme with diffusion freezing, Note that no boundaries exist in this simulation, i.e., no boundary treatment necessary, and no such provisions are included in the network.
We then add a potential gravity field that is based on the distance d and direction of a particle relative to the origin, i.e., [x] and x, respectively as
dv
2 x
(63)
dt
for some gravitational constant g. This gravitation is applied the current velocity as an advection velocity, see Sec. [2.]] where the network task is to predict a particle's position change for a single timestep. Note that in contrast the prior scenarios, we train a single timestep prediction here as the scenario is already challenging as is. We then train all networks for 20 epochs, each consisting of 1000 weight updates with a batch size of 4. We start the training with an initial rollout length of 1 and increase the rollout every second epoch by 1. Finally, we start with an initial learning rate of 10-3 and decrease it every 100 iterations to reach 10-5 5 at the end of training.
B.5 TEST-CASE IV: THREE-DIMENSENE DIMENSIONAL Toy PROBLEMS
To expand our evaluations to three dimensions, we chose to create a simple setup to verify the findings from one and two dimensions. Accordingly, this dataset is designed to serve as a straight forward test case for SPH kernel and gradient interpolations, similar to the problems in test case I, but in three-dimensional space. Generating the data was performed by creating a regular grid of 16 4096 particles on a regular grid spanning the unit cube [-1, 1 with an initial volume per particle of 8/ 8/4096 and a support radius such that each particle has 32 neighbors in this initial constellation under periodic boundary conditions. For the dataset we apply a random three dimensional periodic Perlin noise to the volume of each particle by multiplying the particle volume with random scaling in [-1,1]. While this does result in negative volumes, this is mathematically not an issue as the SPH interpolant and gradient work regardless of particle volume, albeit they only make physical sense with positive masses. We then add an additional normal distributed offset ((0 0.05/) to each particle position add more variation to the dataset. This results in a configuration as shown in Fig. [8] where we generate 1024 such configurations for training and 4 for testing. For each of these configurations we then compute the SPH density and the gradient of the density using a naive gradient formulation, analogous to the toy problems for test case I.
In the training setup, we utilize 4000 weight updates, i.e., approximately epochs, with a learning rate set initially as 10-2 and ending at 10-4 with a reduction in learning rate every 25 weight updates. We perform no data augmentation for this toy problem and utilize as features either the particle density Pi evaluated using a standard SPH interpolant, or normalized volume V/2 = 1/2 chosen such that if a window function identical to the ground truth kernel function was used, the ideal network weights would be 1 for a linear basis.
C ABLATION STUDIES / EVALUATION
This section will discuss the various ablation studies we performed as part of our research and summarize the respective results. We will also provide more in-depth data for all the experiments discussed in the paper, both numerically and visually. See Table for an overview of all experiments.
C.1 TOY- PROBLEMS
Our first set of evaluations focused on test case I, see Appendix B.I where we wanted to investigate how different basis functions could learn simple symmetric and antisymmetric task. In addition
30

---

Published as a conference paper at ICLR 2024
Particle
Particle
Sliced Density
1.00
1.40
0.75
0.75
67
85
825
67
000
0.00
0.25
0.25
-.50
0.5
03
07
-..7
-1.00
-1.0
0.00
1. 20 20 o 100
-1.8 18:3
28 1
Figure 18: Visualization of the setup for test case IV showing initial particle volume (left), resulting particle density (middle) and slice through the density field within the simulation domain. Color mapping indicates, respectively, volume, density and density, with purple indicating low values and yellow high values.
to the evaluations shown in the main paper, we provide additional results for the basis function approaches in the form of tables and an additional ablation study for MLP--hasdd approaches. The latter primarily serves to verify that our implementations of the MLP-based approaches can learn something and how they relate in performance to the basis function approaches. It is important to note, however, that a direct comparison in this case of basis function approaches and MLP on parameter count basis would not be a fair comparison the inductive biases built into basis function test case
approaches would result in parameter counts that are orders of magnitude different specifically in-tune with the inductive biases. Section C.I.1 will discuss the results for learning the kernel function, and Section C.I.2 will discuss the results for learning the kernel gradient. Numerical results are provided in Table and
C.1.1 SPH KERNEL INTERPOLATION
Within SPH, kernel interpolations are important and central aspect, and hence, we investigate learning these interpolations as fundamental and central task. To investigate the ability different basis functions and network architectures to learn SPH kernel interpolation, toy problem, as discussed in Appendix [ We now performed two separate ablation studies, one focused on basis convolutions and one focused on MLP-based approaches.
up simple
Basis Convolutions: Continuous Convolutions using basis functions are inherently built around inductive biases that SPH exhibits, and, accordingly, they should perform relatively well in this setup. To perform our evaluations, we evaluate 10 different basis terms for a single message-passing step architecture with activation functions normalization. The goal here is to determine which basis functions best work to approximate an SPH kernel, and accordingly, the expectation would be that symmetric and smooth methods should perform ideally. Note that we used no window function as using window function, especially if it is identical to the SPH kernel function, would make the learning problem trivial.
Baseline Methods: We evaluated a few baseline methods in this problem, i.e., LinCConv, SplineConv, DMCF and also a naive network build using nearest neighbor interpolation. We ob- served that the best-performing basis in this case was LinCConv as this approach yielded an L error 2.5 10-5, which lower than either SplineConv 13.13- 10-1) and Nearest Neighbor (4.3- 10-1), not able to learn anything in this task due to its inherent antisymmetric nature.
while DMCF Furthermore, we observed that the number basis terms significantly impacted the results, i.c., us- ing only 4 basis terms resulted in noticeably worse result than using 32 terms due to the inherent lack of smoothness for these basis functions.
Using Fourier-basic resulted in much better result than the baseline methods; however, it becomes apparent that increasing the parameter count yielded a worse network perfor- mance. This effect is primarily due to the nature of Fourier Series, higher term counts do not
Fourier Methods:
31

---

Published as a conference paper at ICLR 2024
PODee
2311
3344
1344
2345
46
61
Figure 19: This figure evaluates the correlation between the hidden layer architecture of MLP-based graph networks and the ability of a single-layer pessage-passingg network to learn an SPH kernel function for various basis functions. The x-axes correspond to the number of features per hidden layer, and the y-axis corresponds to the L2 error evaluated using the test dataset from test case I for 4 different frames = [0, 0.512, 1.024, 2.048). with color indicating the number of hidden layers. The plots show the behavior of the MLPCConv (top left), Message Passing PDE (top right), PointNet (bottom left), and GNS (bottom right) network architectures.
improve the smoothness of the convolution but include higher and higher frequency terms. As the underlying kernel function here is relatively low in frequency, including high-frequency terms re- quires the network to learn that they have zero contribution, which difficult to achieve. This is supported by the fact that the Fourier series itself, with two terms by harmonic, shows an apparent behavior of decreasing in performance for every two terms added, i.e., every time a new harmonic is added, the performance decreases, while the Fourier series only consisting of even symmetry terms decreases in performance for every term added. Finally, the Fourier series consisting only of odd symmetry terms was not able to learn only very coarse result in this task due to its inherent anti- symmetry; however, the result is much better than the DMCF approach as the first term is a constant term that is still symmetric.
The Chebyshev basis variant performed reasonably well in this case but did not show much improvement as more terms were included; however, still managed outperform LinCCom for most parameter counts due to the inherent smoothness.
Chebyshev Basis:
MLP-based methods: Considering the results for the MLP-based GNNs, see Figure we can make several observations. Firstly, the PointNer approach fails to learn anything meaningful for any architecture, which is somewhat expected as this architecture does not consider the interconnectivity of the particles and is only provided the position and area of each particle, which is not enough in- formation to reconstruct an SPH kernel interpolation. Secondly, both GNS and MP-PDE performed somewhat comparably with a clear trend that highlights that increasing the number of hidden lay- ers beyond two does not provide much benefit while increasing the number of neurons per hidden layer shows clear and noticeable impact on the results. This is in line with the recommendations of Bunncee--Gozzalee et al.] (2020) to use a hidden architecture that is shallow, i.e., with 2 hidden layers, and wide, i.e., with 128 neurons per layer. Finally, MLPCCom performed better than all other MLP-based approaches primarily due to its inductive bias of representing convolutions.
32

---

Published as conference paper at ICLR 2024
Table 2: This table shows the results of ablation study to find the optimal number of basis terms for learning a purely symmetric SPH kernel interpolation. The results were obtained as the L2 loss for a simple SPH kernel interpolation evaluated on the test dataset of test case I at 4 different timepoints for 4 different network initialization seeds. Each column represents different number of basis terms, with the best behavior for each number of basis terms highlighted in bold.
2
3
5
6
7
8
16
32
Basis
Nearest Neighbor
7.07-10 3
7.06-10 33
3.38-10
1.81-10
1.79-10
1.45-10
1.20-10
1.73-10 3
1.16-10
4.31-10
LinCConv
7.08-10 3
7.08-10 3
2.64-10 -3
1.65-10
8.47-10
8.40-10
3.33-10 4
2.73-10-
2.50-10
8,40-10
DMCF
4.01
4.01
4.01
4.01
4.01
4.01
4.01
4.01
4.01
4.01
SplineConv
7.08-10-3
7.07-10-3
6.71-10-3
5.94-10-3
5.02-10
4.16-10-3
3.69-10-3
2.66-10-3
5.22-10-4
5.13-10
SFBC
7.08-10-3
1.61-10
1.61-10
8.56-10
8.56-10
4.00-10-7
4.00-10
7.15-10-
1.30-10
8.06-10
Fourier (even)
7.08-10-3
1.61-10
8.67-10-*
3.16-10-1
9.40-10-7
2.81-10-5
3.93-10
9.55-10-6
1.80-10-3
6.20-10-5
Fourier (odd)
7.07-10-3
7.07-10-3
7.07-10-3
7.07-10
7.07-10
7.07-10-3
7.07-10-3
7.08-10-3
7.07-10
7.07-10
Chebyshev
7.08-10-3
7.07-10-3
3.65-10-9
3.64-10-5
2.46-10-5
2.43.10-5
1.09-10-4
1.00-10-4
2.26-10-4
6.90-10-5
33
Table 3: This table shows the results of an ablation study to find the optimal number of basis terms for learning purely antisymmetric SPH gradient interpolation. The results were obtained as the L2 loss for a simple SPH gradient interpolation evaluated on the test dataset of test case at 4 different timepoints for different network initialization seeds. Each column represents a different number of basis terms, with the best behavior for each number of basis terms highlighted in bold.
1
2
4
5
6
16
32
Basis
LinCConv
1.00-10-1
1.47-10-1
1.51-10
1.24-10
1.30-10
1.19-10-1
1.32-10-
1.23-10
1.14-10-1
9.30-10
DMCF
9.82-10-2
2.07-10-3
8.36-10-4
9.02-10-4
5.28-10-4
2.34-10-4
2.42-10
7.51-10-4
1.86-10-3
3.50-10-4
SplineConv
1.00-10-1
1.00-10-1
1.06-10-1
1.20-10-1
1.16-10-1
1.07-10-1
1.04-10-1
1.07-10-1
9.73-10-2
9.17-10-2
SFBC
1.00-10-1
4.13-10-3
4.13-10-3
1.42-10-3
1.48-10-3
1.50-10-4
4.37-10-4
4.63-10-4
4.99-10-4
6.63-10-3
Fourier (even)
1.00-10-1
1.00-10-1
1.00-10-1
1.00-10-1
1.01-10-1
1.01-10-1
1.01-10-1
1.01-10-1
1.01-10-1
1.14-10-1
Fourier (odd)
1.00-10-1
4.13-10-3
1.44-10-3
1.09-10-4
1.38-10-4
2.38-10-1
2.88-10-4
3.15-10-4
2.57-10-4
2.06-10-4
Chebyshev
1.00-10-1
9.89-10-2
9.89-10-2
9.53-10-3
2.02-10-1
1.93-10-1
3.11-10-2
2.81-10-2
5.40-10-2
8.66-10-2

---

Published as a conference paper at ICLR 2024
314
114
PAOCOAR
10
This figure evaluates the correlation between the hidden layer architecture MLP-based graph networks and the ability of a single-layer message-passing network to learn an SPH kernel gradient for various basis functions. The x-axes correspond to the number of features per hidden layer, and the y-axis corresponds to the L2 error evaluated using the test dataset from test case for 4 different frames = [0, 0.512, 1.024, 2.048). with color indicating the number of hidden layers. The plots show the behavior of the Message Passing PDE (top left), PointNet (top right), MLPCConv (bottom left), and GNS (bottom right) network architectures.
Figure 20:
C.1.2 SPH KERNEL GRADIENT INTERPOLATION
In addition to kernel interpolations, gradient interpolations are another vital aspect of SPH simu- lation and, accordingly, if method is not able to learn these terms on their own, then the method might not be well suited to learn an overall SPH simulation as, e.g., physical forces generally rely on such gradient terms. Overall, the setup is equivalent to the first toy problem.
Compared to the first toy problem, the results here, see Fig. 2 and Table 3 are notably different. On the one hand, for the kernel interpolation, all baseline methods were able to learn something; in this case, both LinCConv and SplineConv reasonable approximation of the gradient interpolation. On the other hand, DMCF, which was unable learn anything in the prior task, can perform well as the Fourier-based approaches due to the inherent antisymmetry of this basis. These results indicate that learning to represent an antisymmetric interaction is very challenging without an inductive bias help the network. Note that the performance of DMCF did not improve with very high numbers of base terms, i.e., the best performance was achieved with 6 terms, indicating that smoothness in this problem not as necessary as for the first toy problem.
Baseline Methods:
not able to learn any
Chebyshev Basis: While this basis performed well for the first toy problem, it cannot learn any- thing meaningful in this task, even though this method exhibits inherent symmetries. However, it essential note that the fourth term is antisymmetric and clearly performs better than using one more or fewer term, which still indicates that including antisymmetry is important.
In this toy problem, the based methods still performed better than all other compared results; however, the difference between the Fourier methods and other methods is not nearly as significant. A crucial observation that some behaviors observed for the previous toy problem, especially the stepped change in performance, are still exhibited here but in an inverse relationship relative to before. As such, we can clearly observe that the Fourier basis improves notably higher harmonics included and that the performance peaks after including the third harmonic terms. Another important observation that the oddly symmetric Fourier basis does not degrade in performance as notably as the complete Fourier basis.
Fourier Methods:
MLP-based Methods: While for the previous toy problem GNS, MP-PDE and MLPCCom were able to learn the task, in this problem only MP-PDE shows promising results. This clearly highlights
34

---

important part of Published as a conference paper at ICLR 2024
SFBC
10
10
10
Festures Per
Features
Nessage- Passing Steps
Figure 21: This figure shows the relationship of features per layer (x-axis) and message passing steps (color coded) to the ability to learn single-step physics update. The Lg error is evaluated using the test dataset from test case for 4 different frames [0, 0.512, 1.024, 2.048] and across 4 random network initializations, errors bars indicating the lower 5-th to upper 95-th percentile. The DMCF approach is shown on the left, and our proposed Fourier basis approach is shown on the right.
that this problem is much more challenging to learn, especially without the inclusion of several in- ductive biases. While MLPCConv did not perform well, exhibited some learning behavior for larger architectures. Overall, only GNS was able to perform well the inclusion of several in- ductive biases, i.e., providing physical quantities to all MLP operations, did not aid the MP-PDE approach but rather added further additional terms that made the task more challenigng This highlights an important avenue of future research for basis convolutions as the usage of ver- tex MLPs with physical quantities as additional imputs might also be helpful in those architectures; however, this is beyond the scope of our current research.
learn.
C.2 TEST-CASEI
Learning the physics update, i.e., training NNTI, is significantly more challenging than learning a simple SPH kernel kernel gradient. The underlying simulation also requires multiple steps, e.g., to compute particle density and pressure forces, and, accordingly, we utilize multiple message- passing architecture here. We are now interested in how different methods perform, especially with regard to network size, and if there is clear and significant difference between methods that incorporate symmetry and those that do not.
Regarding network scaling, varied the number of message-passing steps and the number of features per layer; see Figure On the one hand, the results indicate that the number of pesssag--passinge steps, once the number of steps is least 3, does not significantly impact perfor- mance. However, for some feature per layer sizes, there is still a trend, e.g., for 4 features per layer. On the other hand, the number of features per layer has much more pronounced impact the network's 's performance and appears to be the driving force of the network performance. Overall, we still see the expected behavior insofar as larger networks perform better than smaller ones. Accord- ingly, we focus the further discussions on a network with 6 passsage-passing steps and 32 features per layer. We also investigated the number of basis terms and found that, like the antisymmetric case, 6 terms perform optimally.
Network Scaling:
Baseline Methods: Compared to the toy problems, the only baseline method that performs reason- ably well is DMCF, which is also the only baseline method that performed well for the antisymmetric case and the only baseline method that includes antisymmetric terms. This result supports our prior observations, i.c., that antisymmetry is essential property for convolutions and that most bascline methods struggle with learning this aspect.
The Chebyshev-basee based convolution, in this case, performs worse than all other methods and, notably, performs worse with larger networks. This generally poor performance is pri- marily due the difficulty in training these networks, i.e., the initialization significantly influences the final network performance, and most initialization does not lead to a well-converging network. Accordingly, further investigations of the initialization of these methods
Chebyshev Basis:
35

---

Published as a conference paper at ICLR 2024
Table 4: This figure shows the quantitative results for an ablation study of basis terms in test case II. All entries are computed by initializing 4 networks for each case at all 4 testing samples from the dataset different time points and then performing density, velocity, divergence, and PSD errors on grid resampling of the particle quantities and finally computing the mean across inference period. Bold indicates the lowest value per row.
inference steps while computing the
Configuration
Window: Müller Div. +
Window: None
Dens.
Vel. +
PSD +
Dens.
Vel. +
Div. 4
PSD +
Basis
Nearest Neighbor
0.040
0.126
2.31
0.170
0.261
5.13
148
LinCConv
0.041
0.130
3.44
45
0.089
0.147
2.32
80
DMCF
0.052
0.124
2.75
52
0.422
0.496
8.62
285
SplineConv
0.043
0.129
2.47
55
0,258
0.285
2.41
152
SFBC
0.031
0.106
2.46
34
0.035
0.098
1.29
41
Chebyshev
0.051
0.135
3.23
73
0.468
0.332
4.20
123
Chebyshev (2nd)
0.058
0.142
3.48
75
0.418
0.296
3.55
106
Gaussian Kernel
0.312
0.377
5.81
203
0.491
0.326
3.60
104
Müller Kernel
0.683
0.637
5.81
354
0.493
0.344
3.96
126
Wendland-2
0.057
0.133
2.53
0.489
0.329
3.71
106
Quartic Spline
0.045
0.151
2.20
82
0.397
0.417
3.99
259
Spiky
0.113
0.165
1.49
0.170
0.199
2.42
111
Bump RBF
0.480
0.433
2.49
125
0.530
0.432
3.48
170
future research, especially for physical simulations. Note that these observations align with prior observations pattern recognition, e.g., by et al. (2022)
Fourier Methods: In this scenario. the complete Fourier Series approach performed better than all other methods and outperforms DMCF across a broad range network sizes, with a more signif- icant gap for smaller network sizes. Interestingly, the oddly symmetrical Fourier Series performs better than all methods except for DMCF, but still significantly worse than the complete Fourier Series. This indicates that simply using antisymmetric basis may work for some basis functions, including symmetric terms, which can significantly improve the overall network performance.
Finally, for the MLP-based GNNs, lines, e.g., LinCConv, but notably worse performance than the antisymmetric basis convolutions. This, again, indicates significant boost the learning capabilities of network by including use- ful inductive biases at the core of the method. Furthermore, these networks required significantly more parameters for these networks, as discussed in Section 4.2 An important final note is the performance of PointNer. While PointNer did not perform well for either toy problem. it performs surprisingly well for this problem and even outperforms MLPCComv, even though it does not con- sider graph connectivity explicitly. This indicates that while this problem is challenging, basing guess for the new particle velocity on the current position and velocity, due the smooth nature of the simulation, can already give reasonable performance.
MLP-based:
performance comparable to most base-
C.3 TEST-CASE
While the main paper already discussed this test case in-depth, the results shown here provide more quantifiable foundation for our conclusions. In this section, we will discuss all the ablation studies and provide numerical results in tables and bar charts, well as visualizations of the inference behavior. The latter serves as a helpful visualization of the difference of the basis functions on top of the quantifiable metrics. Finally, we will expand some evaluations, e.g., we include separate section regarding network initialization stability and how different basis functions perform when evaluated using different time points the test simulation verify that our results are representative.
C.3.1 BASIS FUNCTION EVALUATION
Basis functions are an essential part of our proposed basis convolution methodology, and accord- ingly, evaluating broad set of potential basis functions is crucial part of the evaluation. As
36

---

Published as a conference paper at ICLR 2024
030
Figure 22: This figure shows the inference stability for an ablation study of basis terms in test case II. The lines are evaluated by initializing the network at all 4 testing samples from the dataset at different time points and then performing 64 inference steps while computing the (f..tr.)) density, velocity, divergence, and PSD errors grid resampling of the particle quantities, with error bars indicating lower 5-th to upper 95-th percentile. Note that network initialization seeds are used here well. The color here indicates different basis terms, with the top row showing behavior for using no window function and the bottom row showing the results when using the Spiky window function.
Window: Cartesian
107
Window
Map: Cartesian
307
Figure 23: This figure shows the quantitative results for ablation study of basis terms in test II. All entries are computed by initializing 4 networks for each case at all 4 testing samples from the dataset different time points and then performing (f.l.t.r.) density, velocity, divergence, and PSD errors on a grid resampling of the particle quantities and finally computing the behavior across the inference period. The color indicates the basis term, and error bars indicate the lower 5-th to upper 95-th percentile. The top row indicates the results of not using a window function, and the bottom row indicates the results of using the Spiky window.
inference steps while computing the
37

---

Published as a conference paper at ICLR 2024
formulation is generic and can be used to implement broad range of methods, we focus the most intuitively interesting methods that represent several classes of methods, see Appendix [A.4] The first group of methods is piece-wise bases, which include LinCConv, important prior base- line to compare against, and a naïve Nearest Neighbor based approach as a useful baseline. The second group of methods is Gaussian hases consisting of SplineConv, using cubic B-Spline basis, a quartic B-spline basis, the polynomial Müller Kernel, the SPH-based Wendland-2 Kernel, and a Gaussian, where we chose the methods to represent broad set of choices common in SPH litera- ture. The third group of methods are generic Radial Basis Functions, which include the Bump RBF and the Spiky Kernel, neither of which is Gaussian Gaussian-shaped. The fourth and final group are bases with inherent symmetries and include our Fourier basis, Chebyshev bases of first and second kind, and the DMCF basis. For our evaluations here we consider the inference behavior, see Fig. 22 as well as the average performance over an inference length of 64 steps both numerically, see Table and regarding their distributions, see Fig. 23
Regarding piece- bases, we see notable difference between nearest neigh- bor and linear interpolation; however, this difference is heavily reliant on the choice of window function. When not using a window function, linear interpolation exhibits lower error in all metrics compared to nearest neighbor, which is strongly influenced by the nearest neighbor basis becoming unstable for some experiments. However, when using the Miller window function, the nearest neighbor basis exhibits lower error than the linear basis for all metrics besides the PSD er- ror. Notably, all errors for the nearest neighbor interpolation with a window function are below the values for the linear basis without a window function. Overall, this indicates that nearest neighbor interpolation is useful due to its simplicity but relies heavily on external biases to perform stably. In contrast, linear interpolation performs reasonably well even without a window function.
Piece-Wiie Wise Bases:
Gaussian Bases: Regarding these basis functions, the most crucial observation is that none of these methods work without a window function. i.e., they all exhibit density errors that are an order of magnitude worse than the best method. However, some of these methods perform comparably well without a window function, while others do not. Interestingly, both B-Spline bases and the very similar Wendland-2 basis perform well, while the polynomial Müller basis and Gaussian kernels do not. This difference indicates that similar to SPH literature (Dehnen & Aly 2012), there are significant differences between kernel functions that might not be apparent from their shapes, i.e., the Müllerr basis is very similar in shape the cubic B Spliine but works notably worse. Furthermore, these results indicate that the compactness of the basis function is still a valuable property, as the non-compact Gaussian basis does not work well.
Generic RBFs: Regarding generic basis functions, we did not see strong behavior in any case, which indicates that these functions not ideal bases for our formulation. However, this is hardly surprising as RBF interpolation does not solely adjust the weights for each basis function but instead shifts the centroids and modifies the shapes of the basis terms, neither of which we include in our network. These additions pose an interesting direction for future work; however, optimizing them as part of a neural network may be challenging as these terms are non-linear.
For symmetric without a window function, we find that only our proposed bosdder-based approach remains stable while all other options quickly become unstable. This is not very surprising as, on the one hand, the Chebyshev basis is not compact, and, as seen for the Gaussian basis, compactness appears to be a useful property. On the other hand, the Fourier basis contains terms that decay towards zero at the outside edges. Furthermore, with a window function, all of these bases become stable and show behavior comparable to all other methods, with our pro- posed Fourier-based network performing better than all other bases. Finally, while the divergence error only increases for the Fourier-basse method, all other symmetric bases exhibit a lower error; however, do note that none of the other bases were stable without a window function.
Symmetric Bases:
Overall, we could clearly observe behavior that is consistent for different classes of functions where our proposed network architecture outperforms all other bases in all but single metric. Moreover, these evaluations indicate that other base terms, such as nearest neighbor inter- polation, can still be useful even though they appear less capable than more complex bases, such as cubic B-Splines. Finally, regarding the inference behavior, we observed that methods were either stable, i.e., they remained at low error the entire inference period or steadily got worse over time with no method randomly becoming unstable late into the inference period.
Conclusions:
38

---

Published as a conference paper at ICLR 2024
Table 5: This figure shows the quantitative results for an ablation study of window functions in test case II. All entries are computed by initializing 4 networks for each case at all 4 testing samples from the dataset different time points and then performing 64 inference steps while computing the density, velocity, divergence, and PSD errors grid resampling of the particle quantities and finally computing the mean across the inference period. Bold indicates the lowest value per row.
Basis
LinCConv
Ours
Dens. +
Vel. +
Div. +
PSD +
Dens. +
Vel. 4
Div. 4
PSD +
window
None
0.080
0.137
2.11
67
0.038
0.103
1.24
44
Müller
0.041
0.130
3.44
45
0.031
0.106
2.46
34
Spiky Kernel
0.056
0.180
4.81
87
0.031
0.113
3.06
29
Cubic Spline
0.038
0.135
3.11
54
0.032
0.113
3.11
31
Quartic Spline
0.044
0.150
3.53
63
0.034
0.117
3.36
35
0.037
0.119
2.31
43
0.026
0.097
2.04
32
0.063
0.147
3.13
67
0.025
0.095
1.87
31
C.3.2 ABLATION STUDY: WINDOW FUNCTIONS
Window functions are strong inductive bias in CConv approaches built the ideas of SPH. These window functions aim to (a) ensure compactness, (b) ensure smoothness, and (c) as bias towards the underlying simulation. Compactness is, generally, ensured by defining window functions such that they are equal to for 1g] = 1, e.g., f(g) = [g]. smoothness is, generally, enforced by using higher order polynomials, e.g..., cubic B-splines, and the bias is achieved by using SPH kernels as window functions. Accordingly, it makes sense to evaluate the strength of these individual motivations by evaluating set of basis functions that exhibit varying levels of smoothness and similarity to the underlying SPH simulation, which, for this problem, is the cubic B-Spline kernel.
LinCConv:
For this approach, observe that the best-performing window function, overall, is simple linear window, which only ensures compactness but does not impose any further restrictions on smoothness. The second best kernel function is the cubic B-Spline kernel, indicating that the third motivation, e.g., similarity to the underlying SPH kernel, is also an important effect. Still, this additional bias might introduce other effects that are not desirable. Finally, the other kernel functions, especially the smoother quartic B-Spline kernel, perform worse than the cubic B B-Spline. erucial observation here, however, that the best method regarding the divergence error is not using any window function, which clearly indicates that using window function is not universally beneficial choice and, instead, opens up the ability to fine-tune the network based on the task.
Contrary to the LinCConv approach, the results for our proposed basis are much closer to each other, e.g., the deviation for the density is #0.06 for our approach while for LinCCom the spread is +0.21, indicating that the method is inherently more capable as it outperforms the LinCConv base regarding all metrics with all basis functions. Furthermore, barring the divergence error, the worst choice of window function for each metric for our proposed basis still performs as well, or better, than LinCConv with the respective best choice. For our proposed basis, we also observe that the linear window performs very close the best window function regarding all metrics except for divergence, while not using a window function still results in the lowest divergence.
Fourier:
Based on the results observed here, we conclude that not using window function is an important choice to consider, as not using window function leads to the lowest divergence еттоr in all cases; see also Appendix C.3.I Furthermore, simple choices such as a linear window can outperform more complex and intuitive window functions when choosing an actual basis function. In contrast, some SPH kernels, such as the Spiky window, do not perform as well as others.
C.3.3 ABLATION STUDY: COORDINATE MAPPINGS
Coordinate systems play an essential role in many systems and can encode useful information through their definitions. For example, SPH is built rotationally symmetric kemell functions, ie., functions that solely depend on the distance and not the angle of two particles. Accordingly, using polar coordinates is useful in defining SPH interactions, and using the same coordinate sys-
39

---

Published as a conference paper at ICLR 2024
Figure 24: This figure shows the quantitative results for an ablation study of window functions in test case II. All entries are computed by initializing networks for each case at all testing samples from the dataset different time points and then performing 64 inference steps while computing the (f.It.r.) density, velocity, divergence, and PSD errors on a grid resampling and finally computing the behavior across the inference period. The color indicates window function, and error bars indicate lower 5-th to upper 95-th percentile. The top row indicates results for the LinCConv basis, and the bottom row indicates results of using our Fourier basis.
the particle quantities
Table 6: This figure shows the quantitative results for an ablation study of coordinate mappings in test case II. All entries are computed by initializing networks for each case all 4 testing samples from the dataset at 4 different time points and then performing 64 inference steps while computing the density, velocity, divergence, and PSD errors on grid resampling of the particle quantities and finally computing the mean across the inference period. Bold indicates the lowest value per row.
Window: None
Window:Spiks
Dens. +
Vel. 4
Div.
PSD
Dens.
Vel.
Div.
PSD 4
Map
Basis
Cart.
LinCConv
0.099
0.159
2.88
85
0.035
0.129
3.16
49
NN
0.166
0.291
6.78
150
0.026
0.107
1.49
31
SFBC
0.033
0.097
1.17
35
0.033
0.114
3.19
29
Pol.
LinCConv
0.028
0.095
1.91
28
0.026
0.105
2.33
30
NN
0.125
0.206
3.81
115
0.025
0.108
1.50
33
SFBC
0.059
0.162
0.99
62
0.051
0.172
1.63
49
Pres.
LinCConv
0.281
0.312
6.36
168
0.043
0.138
3.13
64
NN
0.172
0.257
5.14
120
0.027
0.110
1.64
32
SFBC
0.036
0.099
1.35
44
0.034
0.119
3.31
36
tem for neural network that is learn an SPH simulator would also make sense. Furthermore, coordinate systems can be used to make the use of coefficients more efficient, as proposed by Um- menhoter al. (2019), as using regular grid in higher dimensions is not ideal for spherically limited interactions as some weights will have reduced, or even no, influence, However, as we focus our evaluations on vro-dimensional systems, there merit in re-evaluating these biases on our spe- cific problems. Accordingly, we evaluate the influence of choosing coordinate mapping, polar coordinates or a volume-presse peeeeeeeesservnng mapping for the LinCConv approach, our proposed approach and nearest neighbor interpolation approach baseline for window function and the Spiky window function, see Fig. 225 and Table 6
LinCConv:
For the LinCCom approach, we saw the best results without window function for the polar coordinate mapping and the worst outcomes for the volume preserving mapping, while when using window function, the results only marginally different. These results indicate clear benefit to changing the coordinate mapping in some cases, ie.. without a window function. Still, when using window function, the differences negligible. While this is somewhat contrary
40

---

Published as a conference paper at ICLR 2024
Figure 25: This figure shows the quantitative results for an ablation study of coordinate mappings in test II. All entries are computed by initializing networks for each from the dataset at 4 different time points and then performing 64 inference steps while computing the (f.I.t..)) density, velocity, divergence, and PSD errors on grid resampling of the particle quanti- ties and finally computing the behavior across the inference period, Color here indicates coordinate mapping and window function combinations and error indicating lower 5-th upper 95-th per- centile. From top to bottom, the subplots represent using a Nearest Neighbor Basis, LinCCom and SFBC.
all 4 testing samples
Table 7: This figure shows the quantitative results for an ablation study of Fourier terms in test case II. All entries are computed by initializing 4 networks for each case at all 4 testing samples from the dataset at 4 different time points and then performing 64 inference steps while computing the density, velocity, divergence, and PSD errors grid resampling of the particle quantities and finally computing the mean across the inference period. Bold indicates the lowest value per row.
Configuration
Window: None Vel.
Window: Spiky Vel.
Density
Div.
PSD
Density
Div.
PSD
Basis
SFBC
0.046
0.118
1.65
57
0.029
0.109
2.90
28
Fourier (4-Terms)
0.072
0.126
2.33
67
0.026
0.105
2.46
27
Fourier (erme))))
0.041
0.125
2.56
55
0.028
0.106
2.64
30
Fourier (even)
0.373
0.369
5.50
212
0.308
0.336
4.97
177
Fourier (odd)
0.095
0.188
3.81
87
0.032
0.113
3.10
37
Fourier (odd)
0.829
0.523
3.87
NaN
0.328
0.278
2.96
101
Fourier (odd) sgn(x)
0.896
0.561
5.46
NaN
0.431
0.306
3.39
104
prior work, the volume-preserving mapping is mostly intended for three-dimensional applications and not two-dimensional tasks.
Fourier Basis: Similar to observations regarding the window function, Appendix C.3.2 there only a very small difference between different choices of coordinate mappings. This indicates that (a) the influence of different hyperparameters is very low for our proposed basis, and (b) ideal performance can often be achieved while removing inductive biases from the network.
Overall, as
very little difference with window function and for method in general, propose to remove this inductive bias in most cases from network architectures but still include it as can be helpful in corner cases, e.g., without window functions for LinCConv.
41

---

Published as a conference paper at ICLR 2024
Fourier basis
Our Basis
(c) Only even terms
(d) Only odd terms
Figure 26: Visualization of some Fourier basis setups we evaluated as part of our ablation studies.
C.3.4 ABLATION STUDY: FOURIER TERMS
Fourier Series have a long and rich history of over 200 years of research, but at the core, a Fourier series is sequence of alternative sine and cosine terms of increasing frequencies with additional constant function added on top. As such, Fourier series can be defined as
1. cOS
= O, iodd,
6 Fourier (p)
2 + 1] p) 1
(64)
VT sin (TT
ieven.
However, while for an infinite Fourier series, equally many cosine and sine terms exist, this is not necessarily true for finite sequence terms of Fourier Series. For odd number of terms the constant term and (n 1)/2 harmonics of a sine and cosine term each can be utilized. However, common choice in baseline methods, e.g., LinCCom (Ummenhofer al. 2019). using only four terms, which is even and, accordingly. some term needs to be excluded. As representative choices, we investigate either excluding the second harmonic sine term (an odd symmetry term) or the first harmonic cosine term (an symmetry term). Furthermore, we investigate four more series that consists of (a) only even symmetry terms, (b) only odd symmetry terms with an additional constant term, (c) only odd symmetry terms with the constant term replaced by p, and (d) only odd symmetry terms with the constant term replaced by sgn(p). Note that many more choices exist, e.g., an odd symmetry series with offset term at all, but chose limit ourselves this set of series.
Complete Series: Considering the definition from Eqn. 64 we evaluate this series with 4 and terms, Fig. 26a where the 4 term variant drops the second harmonic sine term. As such,
42

---

Published as a conference paper at ICLR 2024
W woo Nowe Map Cortie
30
Wade Spay
044
307
Figure 27: This figure shows the quantitative results for an ablation study Fourier terms test case II. All entries are computed by initializing 4 networks for each case at all 4 testing samples from the dataset at different time points and then performing 64 inference steps while computing the (f.It.r.) density, velocity, divergence, and PSD errors grid resampling of the particle quantities and finally computing the behavior across the inference period. The color indicates Fourier term combinations and error bars indicate lower 5-th upper 95-th percentile. The top row indicates the results of not using a window function, and the bottom row indicates the results of using the Spiky window function.
would expect the term variant to perform notably better regarding symmetric quantities, such as density error, but comparable for antisymmetric quantities. Considering the results in Table and Fig. 27 we observed the expected behavior when not using window function, whereas without window function, the results are much closer to another.
Ours: As antisymmetry is an essential aspect of physical systems, including antisymmetrie terms should improve the performance, and their behavior should dominate the overall learned behav- ior. Accordingly, removing the first harmonic even symmetry term and trading this off with sec- ond harmonic symmetric term is an interesting choice; see Fig. 26a For four-term basis, the two-dimensional series consists symmetric and 8 antisymmetric terms, whereas the alternative choice consists of 10 symmetric and 6 antisymmetric terms. Using this series, we observe signif- icantly improved behavior regarding all quantities without a window function, whereas without window function, this choice of terms performs worse in all quantities. A crucial observation here that the divergence error of this formulation increases by 75% when using a window function, whereas the divergence error only increases by 6% for the alternate formulation. However, this formulation leads to the lowest overall divergence and better behavior without a window function, used this function as the basis for our approach.
Symmetric Series: Using only even terms, i.e., the constant term and cosine terms of increas- ing harmonies, see Fig. 26c leads to 16 symmetric and antisymmetric terms, whereas using the constant term and sine terms of increasing harmonics, see Fig. 26d leads to 10 symmetric and antisymmetric terms. Note that with our separable basis formulation. is not easily possible only construct antisymmetric terms. Regarding the even variant and the modified odd variants, we saw no stable behavior regardless of whether window function was used. This indicates that neither basis is useful choice overall. In contrast, the results regarding the even basis indicate that using a purely symmetric basis not a useful choice for learning physical behavior. However, for the odd series with the constant term, we saw reasonable behavior with and without window function; however, the results were worse than using a modified Fourier Series. This primarily due to the antisymmetric terms, Fig. being only for the sole and terms, i.e.. none of the mixed terms are antisymmetric. This indicates that having antisymmetry alone is insufficient and that having a reasonable ratio of antisymmetric to symmetric terms is essential.
43

---

Published as a conference paper at ICLR 2024
Wices Now Map Catiion
300
Figure 28: This figure shows the quantitative results for an ablation study network initialization stability in test case II. All entries are computed by initializing 32 networks for each testing samples from the dataset at 4 different time points and then performing 64 inference steps while computing the (f.l.trr.) density, velocity, divergence, and PSD errors grid resampling of the particle quantities and finally computing the behavior across the inference period. The color indicates the basis term, and error bars indicate the lower 5-th upper 95-th percentile. The top row indicates the results of not using a window function, and the bottom row indicates the results of using the Spiky window function.
all
Table 8: This figure shows the quantitative results for an ablation study regarding network initial- ization stability test case II. All entries computed by initializing 32 networks for each case all testing samples from the dataset at 4 different time points and then performing 64 inference steps while computing the density, velocity, divergence, and PSD errors grid resampling of the particle quantities and finally computing the mean across the inference period.
Configuration
Window: None Vel.
Window: Spiky Vel.
Density
Div.
PSD
Density
Div.
PSD
Basis
LinCConv
0.114
0.193
3.86
103
0.043
0.147
3.57
74
SFBC
0.037
0.103
1.30
41
0.031
0.112
3.04
30
Nearest Neighbor
0.160
0.266
5.42
150
0.030
0.115
1.67
41
DMCF
0.435
0.479
7.58
268
0.370
0.426
6.52
213
SplineConv
0.260
0.327
3.27
194
0.364
0.732
8.14
332
Chebyshev
0.409
0.351
5.80
157
0.038
0.123
2.86
51
C.3.5 ABLATION STUDY: NETWORK INITIALIZATION
Network initialization and training stability are important aspects machine learning tasks and especially important for repeatability. To investigate the influence of network initialization and our training setup, performed two separate ablation studies where first evaluated 32 different initialization seeds for a set of core methods. Secondly, we investigate if our choice of test frames influences the result, i.e., if the performance of NNTIs is significantly different when applied starting at timestep 0 of simulation or much later in the simulation and if by choice of this starting point, our results were biased towards specific methods.
Considering random initialization, see Fig. 28 and Table g we did not observe any significant differences between this broad set of evaluations and the narrower choice of seeds in the respective ablation studies. While this does indicate that our training setup and initialization process results in repeatable and fair evaluation, we can still observe that the spread of some methods, especially the Chebyshev basis, is much greater than the other methods, indicating that some base functions are more difficult to initialization. We already observed a similar behavior before, e.g., in test case I, where the Chebyshev bases did not perform well many other methods
Random Initialization:
44

---

Published as a conference paper at ICLR 2024
Net
Figure 29: This figure shows the quantitative results for an ablation study of network initialization stability in test case II. All entries are computed by initializing networks for each case at all 4 testing samples from the dataset at a single time point and then performing 64 inference steps while computing the (f.l.trr.) density, velocity, divergence, and PSD errors on grid resampling of the particle quantities and finally computing the behavior across the inference period. The color indicates the basis term, and error bars indicate the lower 5-th to upper 95-th percentile. Each row indicates different initial timestamp in the simulation where the inference started.
most cases, partially due to the initialization. In the future, we hope address the initialization of Chebyshev-based convolutional networks in more depth, but this is beyond current scope.
Testing Frame Dependence: When considering the initial frame choice, see Fig. 29 we did not observe a strong dependence of the initial frame choice to the performance of the network, i.c.. if network performed poorly, it performed poorly regardless of the initial frame. However, a crucial observation be made here is how the errors increase differently based initial frame choice, e.g., when considering the density error, the increase of error per inference step is much lower the further into the simulation the inference is started. This behavior is likely due the underlying SPH data, especially during initialization. As mentioned before, see Appendix B.3] particles in this test case are initialized on rectangular grid of particles with even spacing, whereas particle distributions become disordered and chaotic during the simulation. Accordingly, few training samples include the regular particle sampling from the beginning of simulations, while most of the training samples are during the more random stages of the simulation. Note that this particle disorder is general SPH problem often addressed by particle shifting (Rastelli et al., 2022), a process that aims to regularize particle distributions during simulation to improve numerical accuracy. Regarding the inference behavior, we observed that several methods very quickly diverge away from the ground truth solution but that the error metrics stabilize after some time, where this stabilization happens significantly after the training rollout period. This kind of behavior could indicate that the methods find new stable equilibriums that change depending on the basis function and may be due the inherent numerical
45

---

Published as a conference paper at ICLR 2024
Table This table shows an ablation study performed in test case III find the optimal number of basis terms. These numbers were computed using differently initialized networks, training them to overfit to single training sample and evaluating the point-wise wise distance based solely on this single training sample. The columns indicate the number of basis terms, with lower values being better.
4
8
Configuration
LinCConv Spiky
0.054923
0.050374
0.049219
SFBC @ None
0,052832
0.050328
0.050279
SFBC @ Spiky
0.045858
0.042443
0.040794
Table 10: This table shows the result of an ablation study performed on test case 3 in an overfitting setup find an optimal network layout and show the point-wise distance for proposed Fourier- based architecture with no window function for different network layouts. Note that all values are premultiplied by 102 for legibility, with the lowest value per column being boldened.
Steps
1
2
3
4
5
6
16
32
Features
6.07
5.88
5.83
5.77
5.99
6.02
6.07
6.15
6.37
6.61
5.99
5.76
5.57
5.45
5.37
5.31
5.24
5.21
5.07
4.97
4
5.77
5.35
5.22
5.09
5.06
4.92
4.90
4.83
4.70
4.92
5.63
5.18
4.93
4.85
4.71
4.66
4.61
4.48
4.35
4.54
16
5.42
4.90
4.65
4.46
4.37
4.32
4.31
4.30
4.27
4.42
32
5.34
4.75
4.53
4.40
4.22
4.24
4.20
4.20
4.31
4.36
64
5.17
4.55
4.33
4.24
4.25
4.27
4.94
4.38
4.30
5.27
properties of each basis function. However, investigating this complex relationship is beyond the scope of our work but an interesting direction for future, more theoretical, research.
C.4 TEST-CASE III
This section will provide additional details regarding the architccural ablations for the final test case. The main goal of these results is to show how the network be scaled for better performance and quantify the differences. Finally, we provide results for methods not discussed in the main paper, e.g., DMCF without a window function provide additional context.
C.4.1 ABLATION STUDY: BASIS TERM COUNT
Analogously, test case I, we investigated the influence of the number of basis mance of the networks. To simplify this ablation study, we focused solely on evaluating LinCCom and Fourier-based approach, with and without window. Furthermore, we limited the training to simple overfitting case as we were only concerned with the learning abilities in this problem, the generalization capabilities. To set up this overfitting test, we chose a random frame from the training set and trained neural network with 4 message-passing steps and 32 features per layer with an otherwise identical training setup to the other evaluations.
the perfor-
In this evaluation, saw a minor improvement when increasing the base terms for the LinCCom approach improved by 9% from 4 to 6 base terms and by 2% from 6 to 8 terms. For our Fourier- based approach without window function, performance only increased by 4% from to 6 base terms but only by 0.1% from 6 to terms, see Table I1 Accordingly, we used 6 base terms for all further evaluations. A crucial observation here is that base terms also the choice for optimal behavior the oie-dimensional compressible test case, which indicates that this choice is an overall beneficial choice for improving the learning behavior of basis convolutions.
46

---

Published as a conference paper at ICLR 2024
LiKCCrra
None
SFBC Salay
Forrer None
Fourter Salay
Figure 30: This figure shows the result of an ablation study performed on test case 3 in an over- fitting setup to find optimal network layout. From left to right, the subplots represent using the LinCCon approach with a window function, our modified Fourier basis without and with window function, and using an proper Fourier basis without and with window function. The x-axes repre- sent the number of message-passinnn steps, and the y-axes the number of features per layer. The color indicates the spint-wise distance error from high (white) to low (black), with lower being better.
C.4.2 ABLATION STUDY: NETWORK LAYOUT
The following ablation study we performed was regarding the network layout, similar the ablation study performed in test case I; however, we used the same overfitting setup as for the prior ablation study due to the complexity of this training task. For this ablation study, we evaluated three basis functions, i... LinCConv, complete Fourier Series, and our modified Fourier series from test case with a removed first harmonic cosine term, as well as with no window and with a window for the Fourier based methods, Table I0 and Fig. 30 Note that these results were obtained using single initialization seed, and, accordingly, there might be some notable outliers in the data that are not representative of the overall behavior; however, the general trends still hold.
LinCConv: For the LinCCon approach, we saw a clear and expected trend, i.e., improving the number of ressagg-passing steps and increasing the features per layer improves performance. Con- trary to the results in the one-dimensional case, however, increasing the number of message-passing steps to more than two still significantly improves performance. This difference primarily due to the inherently increased complexity of the problem, i.e., the underlying physical system has a much larger receptive field and requires significantly more SPH interpolants per simulation step. Consequently, increasing both has a notable benefit for the network's performance.
Regarding the Fourier-based methods, we observed similar behavior as for LinCConv; however, the overall performance of the Fourier-based networks is notably better at all sizes of network, especially when not using a window function. Furthermore, for this problem, using complete Fourier Series for small network architectures significantly improves over the modified Fourier Series, indicating that including the symmetric first harmonic term beneficial. However, as the performance is only marginally different for larger networks, continued using this Series.
Fourier Methods:
These ablations show that network with 6 message- paattgge-paiingggggg steps and 32 features per layer, provides reasonable performance for all tested methods without requiring excessive parameter counts. The similarities with test case I indicate that there are general trends in learning behaviors that are worth considering when setting up neural network for given learning task.
C.4.3 BASIS FUNCTION EVALUATION
As the final ablation study, we considered the influence of different basis functions for test case III regarding inference performance. To do this, we trained the respective networks on the entire test case III dataset with an incremental rollout of up to ten timesteps during training, a batch size of four, and decreasing learning rate. As basis functions, we chose the five most important baselines and proposed configurations, i.e., LinCConv DMCF and Nearest Neighbor as well as Fourier and Chebyshev Series based networks. For the comparisons, mean distance of each particle in the prediction relative to the closest ground truth particle, regarding temporal inference behavior, see Fig. 31 and as average inference length of 96 timesteps.
chose the point-wise distance, i.e., the
47

---

Published as a conference paper at ICLR 2024
Table 11: Quantiative results for the different basis functions evaluated in case when learning the physics update. Each value is computed for 4 different network initializations evaluated on all testing samples for 4 different timesteps using an inference length of 96 steps and then computing the mean point-wise distance across the inference period. Lower values are better, and the lowest value per column boldened.
window Basis
Müller
None
Spiky
LinCConv
2.739.10
3.047-10
9.014-10
Chebyshev
2.871-10-3
3.508-10-1
2.941-10
DMCF
2.899-10-3
3.618
5.534-10-1
Nearest Neighbor
3.221-10-3
3.298-10-3
6.197-10-3
SFBC
2.850-10-3
2.876-10")
3.224-10-3
10a
10a
Pointwie Disiance
10
10-3
20
40
80
Inference Step
Basis
LinCConv
Chebyshev
Nearest Neighbor
SplineConv
SPBC
Window
Müller
Figure 31: This figure shows the inference stability for test case 3 when learning the physics update. Each line represents different network initializations evaluated on 4 testing samples for 4 dif- ferent timesteps using an inference length of 96 steps and computing the point-wise distance to the ground truth data after every inference step. The color indicates the basis function, style indicates the window function and error bars indicate the lower 5-th upper 95-th percentile.
LinCConv:
Contrary the prior test cases, the basic approach by Ummenhofer al. 2019) per- formed better than all other approaches, on average; however, during the first 40 inference steps, the performance slightly worse than our proposed Fourier-based approach. An important observation here is the apparent strong influence of the window function on the performance of the basis func- tion, i.e., changing from the Spiky window to the Müller window improved performance by factor of 3.3, mostly due to instabilities during longer inference periods, see Fig. 31
The approach by Pranti et al. (2022) performed very close to other baselines when Müller window was used; however, when using the Spiky window as suggested, the performance decreased by two orders of magnitude. Furthermore, not using a window function in this case did not lead to stable prediction behavior for any inductive bias of antisymmetry can helpful in many cases, it does not guarantee stable behavior and clearly demonstrates how this approach is susceptible to hyperparameter changes.
DMCF:
These observations highlight that while the
Nearest Neighbor: This basis function performs somewhat similar to LinCComv as the choice of window function notably impacts the performance, but overall, performance is pretty comparable.
Fourier-based: While our proposed Fourier-based approach outperformed all other bases in most other comparisons, even regarding the overfitting behavior in this test case, there is significant
48

---

Published as a conference paper at ICLR 2024
Protuate (Ml Tagge
350
18
Esss Tancice Court
Edax Fandiaan
Figure 32: Results of the ablation study for the density (left) and density gradient (right) toy prob- lems in test case IV. Base function color coded.
benefit in this test case. However, our proposed technique still performs better than all other ap- proaches without a window function and, regarding average behavior, only marginally worse than LinCConv. Furthermore, during the first 40 inference steps, our proposed technique outperforms LinCCon, highlighting that there still is a benefit to choosing our proposed technique.
Chebyshev: While the Chebyshev basis struggled in many prior test cases, in this case, it demon- strates performance in line with all other approaches. Furthermore, it shows the best behavior with the Spiky window function, which also exhibits better performance than not using window function for this approach. This clearly highlights that while this method is useful, its performance is still dependent on hyperparameter choices such as window functions.
C.5 TEST CASE IV
This section provides the results of our evaluations in the three dimensional problems in test case IV regarding SPH kernel and gradient interpolation. Section C.5. will discuss results for single-layer networks. setup analogous to Sections C.2 and Section C.5.2 will discuss results regarding three and four-layer networks to evaluate overparametrization. i.c., how the networks learn for a task that requires fewer message passing steps in the ground truth than the network uses.
C.5.1 SINGLE-LAYER NETWORKS
In this toy problem we want to evaluate the ability of different base functions to learn simple SPH kernel and gradient interpolation using single-layer network setup with no bias, analogous to the one-dimensional toy problem, see Fig. 32
Kernel Interpolation: Comparing the results with the results from test case I, see Fig. 2 we observe significant overlap in behavior. For small numbers of terms, i.e., 1 term per axis, all basis functions perform equally poorly with an improvement as the number of terms increases, except for the odd dimensional test case, the three dimensional
Fourier basis. While this was expected in the basis of an odd Fourier basis contains many terms that are symmetric, however, none of the terms along a cardinal direction are symmetric, which prevents the network from learning any reasonable behavior. Similarly, the even Fourier basis performs better than the full Fourier basis for 3 terms, due the inclusion of more symmetric terms, and, analogous the one dimensional problem, the full Fourier basis outperforms the even basis for larger numbers of terms, e.g., for = 8. Furthermore, LinCCom performs better than the nearest neighbor basis for most configurations, whilst both are significantly outperformed by SFBC. However, it is important to note that the overall achieved loss terms here are much higher than in the one-dimensional case, i.e., the lowest loss here is on the order of magnitude of 10 whereas in the one-dimensional case the lowest loss was on the order of magnitude of 10-8
Comparing the results with the results from test case I, see Fig. can also observe a significant overlap in behavior. However, whilst for the one-dimensional case the odd Fourier basis performed better than all other methods, in this case it performs much worse than the full Fourier basis. This, similar to the kernel interpolation task, is due to the three dimen- sional odd Fourier basis containing significant number of symmetric product terms and is not
Gradient Interpolation:
49

---

Published as a conference paper at ICLR 2024
Nevertk varrynee
Figure 33: Results the ablation study for the gradient (left) and density (right) toy problems in test case IV for network with 2 hidden layers using 32 features each. Base function color coded.
Nevertk verrinee
Figure 34: Results of the ablation study for the gradient (left) and density (right) toy problems in test case IV for network with 3 hidden layers using 32 features each. Base function color coded.
purely antisymmetric basis. Furthermore, while SFBC outperforms LinCConv for some parame- ters, the difference is much smaller than for the one-dimensional case and the overall loss terms significantly worse than for the one-dimensional case for a comparable network setup.
Overall, can observe significant overlap between the one-dimensional and three-dimensional toy problems, indicating that the insights can be transferred readily across dimensionalities. Further- more, similar the two-dimensional evaluations, separable bases cannot be fully antisymmetric, leading to notably worse performance for purely antisymmetric tasks, whereas symmetric tasks are still straight forward learn.
C.5.2 MULTI-LAYER NETWORKS
Thus far we only considered the performance of networks in cases in which the perceptual field, and number of message passing steps, between the ground truth and the graph network were identical, for the toy problems, or cases which the network is significantly reduced compared the simulation, for all other problems. However, it is useful to evaluate how the networks perform they are overparametrized in this regard, i.e., what kind of behavior can be observed the graph network performs more message-passee
steps than necessary. To evaluate this behavior we utilize the same setup as for the single layer toy problem, see Appendix [ but utilize a network set up with 2, see Fig. 33 and 3, see Fig. 34 additional layers using 32 features each. As network inputs we consider two variants using either just the normalized particle volume or the normalized particle volume and particle density, where in the latter
the particle density input is not required to compute the ground truth and serves to highlight behavior in case an unneeded feature is provided, i.e., how well network reject false signal.
50

---

Published as a conference paper at ICLR 2024
Two-Layer Network: For the Kernel interpolation we significant drop in performance com- pared to the single Layer case where for two and four base function terms only SFBC was able to learn the correct behavior, and reached similar loss to the single-layer setup. For eight base func- tion terms all networks perform comparably well, but with a loss that is notably worse than for the single-layer setup. For the gradient interpolation terms for SFBC (by approximately an order of magnitude), showing that this problem might be more straight forward learn with a deeper network with non-linearities, through activation functions, than with a single-layer linear setup. However, this only holds for SFBC for four basis terms as no other method was capable of learning any useful behavior. For eight base function terms LinCConv, SFBC and DMCF all out performed the single-layer performance, whereas the Nearest Neighbor base function did not learn any meaningful output.
observed notably lower error for four base
The overall behavior observed here is very similar to the two-layer network regarding SFBC, i.e., the performance in all cases is as good, or better, than for the single-layer network. However, while for the two-layer setup LinCConv and was still able to learn plausible be- havior for four and eight base terms, for the three-layer setup LinCConv did not learn any meaningful behavior for the kernel interpolation.
Three-Layer Network:
When we add the particle density as an unneeded input to the network can clearly observe that for the kernel interpolation only SFBC is able to learn the correct behavior (for two and four basis terms), whereas for eight basis terms only the Nearer Neighhor method learns any meaningful output. For the gradient interpolation we furthermore observed that no configuration was able to learn the gradient for a two-layer setup and eight terms, whilst for a three-layer setup and eight terms only DMCF changed from learning to not reliably learning (as indicated by the large error bar) the problem task. For four terms only SFBC was able to learn any meaningful behavior, consistent with not having the additional input.
Additional Inputs:
Based on these evaluations we can observe that SFBC is more resilient to bes--thanideea learning setups that utilizes either suboptimal, i.c., too deep, network architectures or utilize unneeded net- work inputs. While some methods, e.g., LinCCOnv show some resilience to these changes, they are not as resilient. Finally, while the three-layer setup with density as an additional feature is not a, mathematically, optimal setup for learning the density, SFBC showed lower loss for both two and four basis terms by up to an order of magnitude.
C.6 COMPUTATIONAL PERFORMANCE
So far we only considered the performance of all models regarding their learning efficacy without discussing computational requirements, which are of vital importance for applications in scientific machine learning. To evaluate the computational performance of our method, and baselines, we first constructed set of dummy data points, i.e., 4096 particles in ID, 64 4096 particles in 2D sampled on the unit square and 16 = 4096 particles in 3D sampled on the unit cube. We then computed support radius for all particles such that each particle, under periodic boundary conditions, has 32 neighbors and added a random offset to the particle positions. Based on this set of particles then sampled random feature per particle using a unit normal distribution, (NO, 1), and a respective random ground truth per particle using a unit normal distribution as well. For each set of hyperparameters we then constructed respective neural network and fed the dummy features in, performed a forward pass, compute a loss against the ground truth and performed a backward pass and updated the weights of the network. For each such weight update we computed the computational cost of the forward and backwards pass, as well as the entire weight update process. We also evaluated these updates for a single-layer architecture and an architecture with one and two hidden layers of 32 features each and used 64 measurements per hyperparameter set (totaling 401, 600 timed weight updates total across all constellations). All measurements were done a system with an Nvidia RTX A5000 GPU with 24 GiB of VRAM and an Intel Xeon 6242R CPU with 754 GiB of RAM. We then evaluated (a) the computational scaling for different basis functions, numbers of basis terms and dimensionality. (b) the computational cost of coordinate mapping and (c) the computational cost of using a window function.
To investigate the computational scaling we evaluated set of basis functions in one, two and three dimensions using identity, polar and preserving coordinate mappings and using the Muller window function, well as using no window function, see Fig. 35 and 36 for the results for the forward
51

---

Published as a conference paper at ICLR 2024
(13.3
21.00
Figure 35: Performance evaluation of the forward pass of various basis functions and approaches using an identity coordinate mapping and no window function in one, two, and three spatial dimen- sions. Basis function color mapped and row indicating dimensionality, shown for network using single layer (left), single hidden layer with 32 features (middle), and two hidden layers with 32 features each (right). All measurements milliseconds.
and backward passes. Considering the forward and backward passes, overall, we observed no sig- nificant difference in overall behavior but did measure a slightly faster backward pass than forward pass, which is the expected behavior as the backwards pass also involves the same convolutions, but transposed and in reverse, but can reuse some results from the forward pass. Consequently, we will not consider the passes separately for all remaining discussions. Regarding the aggregate results, see Fig. 38ltt we see very small difference in performance in all dimensions between the basis, with the only significant outlier being the DMCF method in two dimensions. Otherwise, all aggre- gate values are within +6%, indicating only a small difference in performance between methods. Looking the influence of basis terms on computational performance, see Fig, 36] can observe that most configurations for eight and fewer base terms perform very similarly, i.c., within measure- ment accuracy for most measurements, with the only outlier being the DMCF method. For larger networks, i.e., in three dimensions, and for the network with two hidden layers, we see increase based on the number of basis terms that is approximately cubic, i.e., ( which expected as the number of weights increases with (n2)). However, it is important to note that for smaller num- bers of base terms, e.g., and fewer the scaling sublinear due to the GPU being underutilized for smaller term counts. Similarly, observe that in dimension only SFBC and Chebyshev show a significant increase in computational requirements for large numbers of base terms as they are the only bases that require measurable increase all other methods are bottle-necked in performance based on other parts the network.
Regarding computational requirements of the coordinate mapping, Fig. 37 marginal increase in computational requirements in two dimensions and three dimensions for all numbers of basis terms. Considering the aggregate results, see Fig. 38 middle, the overall differ-
clear but
52

---

Published as a conference paper at ICLR 2024
14.07
Figure 36: Performance evaluation of the backward pass of various basis functions and approaches using an identity coordinate mapping and no window function in one, two, and three spatial dimen- sions. Basis function color mapped and row indicating dimensionality, shown for network using single layer (left), single hidden layer with 32 features (middle), and two hidden layers with 32 features each (right). All measurements milliseconds.
Dimensionally
Dimenslorally
102
3x
Idintit
Coordinaeee
Uimmenhater
Figure 37: Performance evaluation of weight update for different coordinate mappings for SFBC and no window function two and three spatial dimensions. Coordinate Mapping color mapped, results aggregate over all three network architectures, and both window functions tested. All mea- surements milliseconds.
ence in computational cost between different coordinate mappings is relatively low at 7% in two dimensions and 0..4% in three dimensions.
53

---

Published as a conference paper at ICLR 2024
96-7
Figure 38: Performance evaluation of weight update aggregate results across all relevant hyper- parameters for one, two, and three dimensions (indicated on the x-axis). Color mapping indicates the basis function (left), coordinate mapping (middle), and window function (right). Note that for one dimension, we only evaluated an identity coordinate mapping. All measurements in millisec- onds.
12
Figure 39: Performance evaluation of weight update for different window functions for SFBC and an identity coordinate mapping in two and three spatial dimensions. Window function color mapped, values aggregated over all coordinate mappings and architectures. All measurements in milliseconds.
Regarding computational requirements of window functions, see Fig. 39 we observed that for one dimensional networks there is slight, but consistent, overhead imposed by the window function. However, this increase relatively small and disappears for two and three dimensions, resulting in mean difference of performance of at most +0.7%, with and three-dimensional networks showing a difference less than +..3%.
Based these results, there is notable overhead of using SFBC approach, compared to LinCConv, for any network architecture used in our other evaluations. Furthermore, there is no significant overhead imposed by using either coordinate mapping or window function.
D CONCLUDING THOUGHTS
In our appendix, we provided a broad range of data regarding many different ablation studies high- light the strengths and weaknesses of several methods. Overall, our proposed Symmetric Fourier Basis Convolution (SFBC) approach performs either ideally or close ideally across broad range of problems and scenarios regardless of the hyperparameters chosen, indicating that our method is versatile and useful method overall. We hope that our evaluations and datasets inspire future research and help with developing novel and exciting solutions to these challenging physical problems.
54


Reviews:
[{"review_id": "rxOJJsOEql", "reviewer": "Reviewer_Ac8w", "summary": "This paper introduces Symmetric Fourier Basis Convolutions (SFBC), a new method for learning Lagrangian fluid simulations using graph neural networks. The key idea is to construct convolutional filters using Fourier series as separable basis functions, incorporating both even and odd symmetries. Through extensive experiments on three fluid simulation datasets, the authors demonstrate that SFBC outperforms prior methods like LinCConv and DMCF in terms of accuracy and stability. A major contribution is showing that previously used inductive biases like window functions are no longer needed with the Fourier basis. The paper provides a generalized framework for continuous convolutions, rigorous evaluation of design choices, and highlights the benefits of symmetry and smoothness in the convolutional bases.", "questions": "N/A", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The paper aims to tackle an important and challenging problem in ML for physics - modeling Lagrangian fluid mechanics.\n\n2. Leveraging ideas like symmetry, smoothness, and Fourier bases to inject useful inductive biases into graph networks is logically sound and extends prior work nicely. In addition, it is nice to unify symmetry and antisymmetry under the same framework.\n\n3. The extensive experimental methodology covering diverse design choices and evaluations on three distinct test cases is a major strength. The baselines are strong enough and the evaluation is rigorous. While the theoretical novelty is incremental, the engineering rigor and state-of-the-art results are valuable and significant contributions.", "weaknesses": "1. The theoretical novelty is somewhat limited as the core concepts like symmetric bases are adapted from prior work in other domains. Additional theoretical analysis could further strengthen the approach.\n\n2. The evaluations are extensive but restricted to a specific type of fluid simulation problems. It is hard to extract insights that can generalize to a broader scope. Testing generalization on more diverse physical systems could better establish applicability.\n\n[1] Li, Z., Huang, D. Z., Liu, B., & Anandkumar, A. (2022). Fourier neural operator with learned deformations for pdes on general geometries. arXiv preprint arXiv:2207.05209.\n[2] Pfaff, T., Fortunato, M., Sanchez-Gonzalez, A., & Battaglia, P. W. (2020). Learning mesh-based simulation with graph networks. arXiv preprint arXiv:2010.03409.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "O9KVwinZRD", "reviewer": "Reviewer_akhW", "summary": "The study introduces a novel approach employing continuous convolution with symmetric Fourier basis functions, effectively leveraging the problem's inherent biases. The research demonstrates the superiority of using symmetric/anti-symmetric basis functions over the previously suggested explicit weight-tying mechanism. Additionally, the paper offers comprehensive ablation studies on various hyper-parameters, providing in-depth insights into their impact on the outcomes.", "questions": "1. Regarding the proposed technique for obtaining symmetric/anti-symmetric basis functions (Eqn 9 and 10), it is essential to clarify whether these equations pertain to symmetry along the x-axis, y-axis, or a combined symmetry along both axes [x,y].\n Additionally, what is the assumed domain of the basis functions $b_x$ and $b_y$? Furthermore, it is crucial to determine whether a given set of orthonormal basis functions retains their orthonormality after applying this technique.\n\n2. How does having more parameters than DMCF help here? Because even if you are learning more coefficients, the basis functions are constrained. So, in terms of expressivity, wouldn't the two methods be the same?\n\n3. Did you directly use the even and odd Fourier Basis for symmetric and anti-symmetric basis functions (after proper modification of the 0 frequency)? If not, what is the explicit equation of the SFB? And if so, Do you think, with proper regularization, filters with Fourier basis will be able to learn proper symmetric and anti-symmetric kernel?\n\n4. In Figure 2, we can see that with a high number of basis functions, Fourier even (or odd) is outperforming DMCF. What might be the reason for it?", "limitations": "", "rating": 5, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "The research presents a method for acquiring symmetric/anti-symmetric basis functions successfully applied to selected problems. The study conducts a detailed ablation study to compare the proposed technique with related methods and the choice of hyperparameters. The authors also additionally present a novel dataset that can be used for the aforementioned problem.", "weaknesses": "1. The contribution of the work requires a clearer elucidation. Although the work introduces a symmetric Fourier Basis, it does not explicitly define the analytical form of the basis functions being utilized or distinguish them from Fourier sine and cosine series. Additionally, essential questions pertaining to the rationale behind the new proposed technique remain unaddressed (please refer to the listed questions).\n\n2. I find it puzzling that the study did not incorporate \"WaterRamps\" and \"Liquid3d\" from prior research, considering their complexity. Utilizing these challenging scenarios would have been more suitable for demonstrating the model's effectiveness.\n\n3. Improvements in some cases seem marginal.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "0nsJxkuEwn", "reviewer": "Reviewer_T2aE", "summary": "The authors propose a general formulation for continuous convolutions using separable basis functions and evaluated on 3 different datasets. They demonstrate that even and odd symmetry are the critical components for success. With the proposed Fourier-based network, the inductive bias, like the window function, is no longer necessary.", "questions": "1. \bThe result of the proposed separable basis is demonstrated on the 1D and 2D tests. Since the method uses the outer product, will it cause issues with the scalability to higher dimensions?\n\n2. What is the computational overhead compared with other baseline models?\n\n3. What is the limitation of the current methods?\n\n4. The author cites Fey 2018 works several times. In the cited work, the high dimensional spline basis is also the product of the 1d spline basis, sharing similar properties of the proposed separable basis. Can you elaborate on the difference between the proposed one and Fey's method?", "limitations": "", "rating": 6, "confidence": 5, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "Quality: the paper is rich in detail, more than 40 pages in total and give sufficient background and relevant math for understanding the problem.\n\nOriginality: The author proposed an SFBC approach that works better than other CConv-based methods. The author also shows that with the proposed structure, the window function is not necessary.\n\nSignificance: The proposed Fouries-basis network is part of a larger group of symmetric methods that opens up future research.\n\nClarity: The content is self-contained and easy to follow.", "weaknesses": "Some details need to be clarified. See questions. Moreover, there are no movies for the learned dynamics to check for temporal coherence. The unique contribution compared with previous research needs to be elaborated.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "AhEoBAVMYG", "reviewer": "Reviewer_Acuo", "summary": "- The paper proposes a generalized formulation of continuous convolution approaches using separable basis functions. The coefficients of the basis functions are learned from the training data.\n- By using the Fourier series as the basis functions with both even and odd symmetry, the symmetries are built in the in the convolutions which are beneficial for accurately learning Lagrangian flow physics. \n- The method is evaluated on multiple cases:\n    - The Symmetric Fourier Basis Convolution approach is compared against against other CConv approaches in a toy problem.\n    - It’s compared against multiple MLP-based GNNs methods, using a compressible one-dimensional problem \n    - It’s evaluated on a two-dimensional closed domain simulation, which shows the proposed method can predict over a long time frame with inference stability. In this 2D case, the paper also explores the influence of window functions, Fourier terms, and coordinate mappings.\n    - The free-surfaces is also evaluated, using a fluid blob collision scenario, to investigate how different basis terms perform with partially occupied support domains.", "questions": "- When applying this method to 3D scenarios, is there a concern that the weight matrix (W) might contain too many parameters, and intermediate matrix B (in Eq. 20) might become too large, causing challenges for training?\n- When comparing the proposed method with related works, it might be beneficial to include runtime performance metrics as well, which can include information such as the number of iterations required for training, the time taken for each training iteration, and the time take for generating one prediction.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The proposed method of using Fourier basis functions provides an inherently symmetric and smooth continuous convolution approach, which shows good results on several testing cases.\n- The paper provides additional details, analysis and ablation studies in the appendix to better show case the method and the comparison with other continuous convolutional approaches.", "weaknesses": "- Most real-world applications of fluid simulation are in 3D cases. It would be good if the paper has some 3D examples or discussions regarding how to extend the method to 3D cases.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (poster)