Boundary Graph Neural Networks for 3D Simulations
Authors: Andreas Mayr
Sebastian Lehner
Arno Mayrhofer
Christoph Kloss
Sepp Hochreiter
Johannes Brandstetter
Venue: ICLR 2022 Submitted (2022)
Keywords: Simulation
Graph Neural Network
Boundary Conditions
Granular Flow
Physics Application

Abstract: The abundance of data has given machine learning considerable momentum in natural sciences and engineering. However, the modeling of simulated physical processes remains difficult. A key problem is the correct handling of geometric boundaries. While triangularized geometric boundaries are very common in engineering applications, they are notoriously difficult to model by machine learning approaches due to their heterogeneity with respect to size and orientation. In this work, we introduce Boundary Graph Neural Networks (BGNNs), which dynamically modify graph structures to address boundary conditions. Boundary graph structures are constructed via modifying edges, augmenting node features, and dynamically inserting virtual nodes. The new BGNNs are tested on complex 3D granular flow processes of hoppers and rotating drums which are standard components of industrial machinery. Using precise simulations that are obtained by an expensive and complex discrete element method, BGNNs are evaluated in terms of computational efficiency as well as prediction accuracy of particle flows and mixing entropies. Even if complex boundaries are present, BGNNs are able to accurately reproduce 3D granular flows within simulation uncertainties over hundreds of thousands of simulation timesteps, and most notably particles completely stay within the geometric objects without using handcrafted conditions or restrictions.

Under review as a conference paper at ICLR 2022
BOUNDARY GRAPH NEURAL NETWORKS FOR 3D SIM- ULATIONS
Anonymous authors Paper under double-blind review
ABSTRACT
The abundance of data has given machine learning considerable momentum in natural sciences and engineering. However, the modeling of simulated physical processes remains difficult. A key problem is the correct handling of geometric boundaries. While triangularized geometric boundaries are very common engi- neering applications, they are notoriously difficult to model by machine learning approaches due to their heterogeneity with respect to size and orientation. In this work, we introduce Boundary Graph Neural Networks (BGNNs), which dynam- ically modify graph structures to address boundary conditions. Boundary graph structures are constructed via modifying edges, augmenting node features, and dynamically inserting virtual nodes. The new BGNNs are tested on complex 3D granular flow processes of hoppers and rotating drums which are standard compo- nents of industrial machinery. Using precise simulations that are obtained by an expensive and complex discrete element method, BGNNS are evaluated in terms of computational efficiency as well as prediction accuracy of particle flows and mixing entropies. Even if complex boundaries are present, BGNNs are able to ac- curately reproduce 3D granular flows within simulation uncertainties over hundreds thousands of simulation timesteps, and most notably particles completely stay within the geometric objects without using handcrafted conditions restrictions.
INTRODUCTION
Deep learning (Krizhevsky et al. 2012) dramatically changed scientific fields such computer vision, natural language processing, or the medical sciences. More recently, deep learning research has been expanded towards physical simulations such as fluid dynamics, deformable materials, or
acrodynamics (Li et al. 2018 Ummenhofer et al. 2019 Sanchez--Gnzzler et al. 2020 Pfaff et al. 2020). The progress of deep learning in physical simulations was often driven by Graph Neural Networks (GNNs) (Scarselli et al. 2009: Defferrard et al.] 2016, Kipf & Welling 2017). which
proved effective when modeling interactions between many entities via forward dynamics (Battaglia et al. 2018).
Here want to focus learning practically relevant granular flow simulations. Granular flows ubiquitous in nature and industrial processes. Pharmaceutical powders, plastic granulates, or rocks obtained by mining are just some examples of granular media that are used in industries and which are processed in multitude of different flow states. Consequently, simulations of granular flow processes are required for the design of many industrial processes. Such simulations allow the optimization of devices and machinery in which particle flow is essential. These simulations are needed for wide range of materials and material mixtures, which can have different cohesion or friction properties. The utilized machines constitute the boundary conditions for the considered granular flow processes. Many conventional simulation approaches, that are based on solid mathematical theory and which take complex boundaries into account, model these boundaries by triangularizations. Therefore, industrial machinery is often represented by triangular meshes, which can be considered a standard geometric description in the engineering field. Figure I exemplarily visualizes triangularized surfaces for two standard components in industrial setups, namely hopper (left) and a rotating drum (right).
In this work, try to tackle the problem of accurately modeling granular processes when trian- gularized boundary surfaces are present. We therefore suggest a simple extension conventional GNNs: Boundary Graph Neural Networks (BGNNs). BGNNs dynamically modify graph structures

---

Under review as a conference paper at ICLR 2022
to model interactions of particles and triangulated boundaries. This is achieved by inserting virtual nodes, adding edges and modifying node and edge attributes.
We test the effeetiienness of BGNNs complex 3D granular flow simulations of hoppers and rotating drums. The data for BGNN training is obtained by precise but potentially time-consuming simulations. The accuracy of BGNN models is measured via various aggregate quantities, namely averaged particle positions, flows, and mixing entropies. We ensure that our BGNN models have high simulation quality by requiring high accuracy of the models in terms of these aggregate quantities. BGNN predictions have the potential to be considerably faster than traditional granular flow simulation methods while keeping the same precision. Code, training and test data will be made public upon publication.
The main contributions of this work are:
We introduce Boundary Graph Neural Networks (BGNNs), which dynamically modify graph structures via boundary graph structures to allow an accurate modeling of particle interactions with triangularized boundaries.
We implement BGNNs for learning complex 3D granular flow simulations of hoppers and rotating drums.
We assess the BGNN simulation quality via the difference of relevant physical quantities between model predictions and simulations. We show that BGNNs able generalize granular flow dynamics over hundreds of thousands of timesteps while having the potential be considerably faster than state-of of-the-rrt simulation methods.
Figure 1: Triangularized boundary surfaces the hopper (left) and the rotating drum (right). The accurate description of rather simple, curved geometries requires relatively large number of triangles of different shapes and sizes.
BACKGROUND
Dataset. Due the fact that there are no governing equations for granular flow like the Navier- Stokes equations for fluid flow (Faccanoni & Mangeney 2013). method similar to molecular dynamics, the Discrete Flement Method (DEM) (Cundall Strack 1979) is used for simulating granular flow. The key idea of DEM is to represent granular media by discrete particles (e.g. spheres or polyhedra), which interact by exchanging momentum using established contact models. LIGGGHTS (Kloss al. 2012 to simulate particle flow for a wide range of materials and complex mesh- based wall geometries and, which therefore enables the simulation of many industrial processes. Consequently. used LIGGGHTS to generate training, validation and test set simulation trajectories for granular flow within different hopper and rotating drum environments.
App. A is an open-source DEM implementation, which is able
Time transition model. We build our method upon Sanchez-Gonzaler et al. (2020) and use the semi-implicit Euler method to numerically integrate the equations of motion using model-predicted acceleration. The time-transition from time to time + 1 is given by gr+1 o Ar &f and xftl = + At d where x is the particle location, and & the particle velocity. To calculate the time-transition f+1 the particle acceleration & is predicted.
2

---

Under review as a conference paper at ICLR 2022
Graph Neural Networks. We consider graphs (V. &), with nodes V, V and edges eij E &, where N-dimensional node features pv; R attached to each of the nodes. Whether the graph G contains an edge between a pair of nodes (v;,V,) depends on the distance between the nodes:
Eij E = djvi.vy) & cut-off,
(I)
where the cut-off' radius is usually a hyperparamter of the model. Edges might have M-dimensional edge features duj E RM attached to each edge etz, Graph networks are designed to learn from graph- data (Scarselli et al. 2209 Kipf & Welling 2017; Defferrard et al., 2016 Battaglia ct al. (018). Message passing networks (Gilmer al. 2017) are specific type of graph neural networks and usually consist of three different types of layers: node and edge feature embedding layers, the core message passing layers, and read-out layers. Message passing iteratively updates the
embeddings of edges (m:j) and nodes (hz),, i.e., the embeddings of aij and pv,, at edge eij and node
V/ via:
mij = (br; fej. mij),
h/ = ds (h;, Dei,ee mit).
(2)
where the aggregation De node Vy in Eq. is across all nodes that connected node Vi via an edge éij. Typically O represents mean or max operation. The learnable functions & and di commonly presented by Multilayer Perceptrons (MLPs). Equation ( describes in compact form the computation and aggregation of messages, and the subsequent update of node embeddings. The final node embeddings are used for predictions via read-out layers.
BOUNDARY GRAPH NEURAL NETWORKS
Boundary Graph Neural Networks. We introduce Boundary Graph Neural Networks (BGNNs) for modeling the time transition dynamics in simulations within complex geometries. In BGNNs, each graph node v, associated to a particle with location Vy, velocity कv; and acceleration w. which is similar to SSancaez-Gonzzzz to include boundaries (see Fig. 2p. BGNNS dynamically add virtual nodes v V for boundary regions, iff the corresponding boundary region is within cut-off radius to any other particle. We augment the set of edges & by boundary edges ēرj giving an enhanced edge set & with etj E & and ēij e &. Analogously to Eq. (1), the existence of patticleppaticle edges eij and particle-boundary edges ěij is determined via:
et al. (2020). BGNNs modify and enhance the graph structure
e1 E c = d(v,,)) vi) < cut-offfe
(3)
ει E E c = d(vi, v,) & ccttofffe.
(4)
Note, that the cut-off radii cut-off, and cut-off, not necessarily the same, and, : V R. while d: V x Φ - R, i.e. bidirectional edges are used between real nodes and unidirectional edges used between real and virtual nodes.
In order to include more information about boundary surfaces into particle-boundee interactions, V-diiensional node features that encode information about the inclination of triangles in space are concatenated with the existing node features pv; E R. Additionally, coordinate information is used both for existing nodes (X = {Z x),-)) () well as for virtual nodes (_²- = {8 sp. ... "..,,,,,,.. & )) For virtual nodes, the additional coordinates žp, are chosen such that they minimize the distance between points from boundaries and real particles. The resulting set of node features P and node coordinates X are:
P_ (P P Pay. )
x == (Frg. (.P Bay 3.p-1 (5)
where P:: E I N+N and i;; RR3 denote the elements of P and X, respectively. Similarly to above, message passing updates the embeddings of edges (mizj) and the embeddings of nodes (h₂) via
m, El hr fuj. mi)) f,.
h = (h., Daye mij)
(6)
where the aggregation paye & at node v; in Eq. (6) is across all real or virtual nodes that are connected to vi via an edge ēij. Similar to Gilmer et al. (0077) and Satorras al. (202)), make of
3

---

Under review as a conference paper at ICLR 2022
Per Xa1
1/
Pige Nep
Prp Xep
Ag
d( 0u)
Piys Rio
Piy Xep
Figure 2: Dynamic modification of the graph edges (red lines) and nodes (red points). Calculation of the distance d((22, Vo) between a real particle node v2 and the triangle corresponding to virtual particle node Vo- Right: Insertion of an additional edge between Vo and v2 and representation of the nodes in terms of the corresponding node features Puy. X1, and Piy. Ki, for real and virtual node features, respectively.
pairwise distances (ja;; - , and d; 2ty and deterministic functions thereof). These are for BGNNs between real and between real and virtual particles and we pass this information to the graph network edge attributes å.j. for which an initial edge embedding my is determined via an edge embedding layer. The final node embeddings are used for the predictions via the read-out layers. For aggregation D. we the mean. The main challenges of the implementation of BGNNS the modification of the graph structure towards a dynamic boundary graph structure to include boundary information via the introduction of modified edges ēij. modified node features pi, and coordinates for virtual nodes Z... In the following section, we discuss how ēij. Pi, and i; are calculated.
DYNAMIC BOUNDARY GRAPH STRUCTURE
boundary graph structure, which a dynamic modification of the graph structure needed for two reasons: (i) A static graph structure, which inserts many particles (number proportional to the surface area) for every boundary surface, may result in large computational costs 3D scenes. (ii) For certain time frames, only some parts of the mesh might be relevant and computations can be saved. For example, as long as particles are in free fall in a container far away from the bottom, the mesh part describing the bottom is irrelevant for the next time step.
One virtual particle is enough to describe particle boundary interactions. Every particle "sees" most one virtual particle representing the boundary surface area, namely that particle which has the shortest distance. Since interaction strength decreases continuously with the distance, it is ensured that for every barudclebboundry interaction that boundary point with the largest contribution is considered. Table shows average numbers of nodes /V/. as well as average numbers of boundary edges | and the relative increase in edges (ratio of the number of added wall edges the total number of particle edges). The scalability of BGNNs would suffer if more than one particle per particle-boue boundary interaction surface was considered.
Table I: Growth of the number of edges due to boundaries in the graph. The table shows statistics across the training trajectories of non-cohesive particles in standard setting in hopper and drum experiments. For each trajectory the frame with maximum relative increase in the number of edges due to virtual particles [ has been selected representative frame. This is done since are interested in the maximum effect additional virtual particles have on the memory requirements. Number of particles |v]. number of additional virtual edges (8]. and % increase are listed.
Experiment
[v]
% increase
Hopper
1113 738
5475 3547
72.2
Drum
3283 + 282
1678 = 188
54.8
4

---

Under review as a conference paper at ICLR 2022
A dynamic boundary graph structure is obtained in three stages:
Calculating distances between real particles and triangular boundary surface areas order to decide if virtual particle node needs to be inserted into the graph, and subsequently obtain additional edges ěιj
Obtaining positional coordinates for virtual nodes Be, as representatives of the relevant fraction of the triangular surface area B,,
Modifying node features pi by including normal directions to facilitate the learning of geometric relationships between particles and boundary surfaces
Modified edges ěι via calculation of distances to boundaries. To obtain modified edges ēij. boundary particles are dynamically inserted, but only if real particle is close to the corresponding boundary. The insertion of unnecessary edges into the graph is avoided. Such edges would connect nodes of real and virtual particles, although they far apart. To decide whether virtual particle has to be inserted requires the calculation of distances between pairs of real particles and mesh triangles. Specifically, the squared distance between the particle center and the closest point on the mesh triangles is calculated (adopted from (Eberly 1999)). For this purpose, a location on a triangle t is parametrized by two scalar values и and v
t(u v)
u e + vej.
where u 2 0, v 2 O, and + & 1, b represents one of the nodes of the triangle, and, ep and e are vectors from b towards the other two nodes (see Fig. 3) The minimal Euclidean squared distance of the point p to the triangle given by the optimization problem:
d min qu,, v) = [[ v) -pll
S.L. u 0 v 2 0, u+v <
The minimizing arguments u' and v' parametrize the closest point t(u', r') of the triangle to the point p. As indicated in Fig. 3 seven cases have to be to distinguished: one case (cO) in which t(u', v') located within the (closed) triangle, three cases (c1, c3, c5) which t(u', v') is located on one edge of the triangle (including the edge corner points as special cases), and three cases (c2, c4, c6) in which t(u,, v') located on one of two edges (including the triangle corner points as special cases). Whether a virtual particle is inserted is determined by Eq. 44) and the particle-triangle distance d.
(u, 1)
case
ffu, n) = const
case 1
case 3 5
x))
case
b
Co
case 6
case 4
case
Figure 3: Visualization of point- triangle distance calculations in 3D. The triangle is represented by a parameterized function & (u, v) + u eo + v ei with u 2 0, v 2 O, u + v (indicated by the orange area). Level sets of (u, v) are indicated by ellipses and describe the squared Euclidean distance of a triangle point (u, v) to the point p. for which we compute the minimum distance.
Positional coordinates T., for virtual nodes. After having modified the edges ēij via the calcu- lation of distances between real particles and triangular boundary surface areas, we now have to decide where to place the virtual particles on these boundary surface areas. For virtual particles the coordinates é,, correspond to (u', v"), i.c., the closest points on the triangular boundary surface areas with respect the interacting real particle.
5

---

Under review as a conference paper at ICLR 2022
Modifying node features P. by including boundary normal directions. Typical granular flow simulations comprise substantially more particle particle interactions than particle boundary interactions, which may impede the learning of particle- boundary interactions. In Kipt et al. (2018) the problem of qualitatively different interactions is addressed by introducing a dedicated message generating network for each interaction type. We avoid such extensions of our model by means of the following two approaches. First, we introduce additional node features, such that the neural network able to distinguish the different types of nodes. Second, we adapt the weight initialization of the node feature embedding V such that the embedding network can be trained with larger values for the additional features. Consequently, the network can learn different dynamics for particle particle and partiele- boundary interactions. The additional node features are: (i) type feature, i.e., a binary indicator of whether a node represents a particle that is real or virtual, and, in the latter case, (ii) the components of the normal vector (see App. B for more information on an orientation-indeeendent representation of the normal vectors) of the triangular surface areas (null vectors for real particles).
5 RELATED WORK
There a rich body of literature on applications of Deep Learning in the context of physics simulations. Most notably related to BGNNs are the works of Saanhezz-Gozzaler et al. (2020). |mmennofer et al. (2019). and, [i et al. (2018), all of which propose methods of learning particle simulations without enforcing constraints. These approaches can be contrasted to works like Ladický et al. (2015) or Schenck & Fox (018) that utilize strong inductive biases. Ladicky et al. (2015) construct features for Random Forest Regression that are influenced by Smooth Particle Hydrodynamies (Gingold & Monaghan| 1977 Lucy 1977). [Schenck & Fox (2018) construct a differentiable fluid dynamics
network that is closely related the Position Based Fluids method (Macklin & Müller (013). Importantly, both methods are built on the assumption that the governing equations of the system are known, which typically not the case for granular flow dynamics.
Complex mesh-based wall geometries have been employed to compute updates for nodes of the mesh itself (Plaff al. 2020). In contrast, our simulator utilizes the mesh to represent static boundaries in a highly efficient way. We share the opinion of Saanaann- Gonzalez et al./ (2020) that
the network architecture with continuous convolutions as suggested by Ummenhofer et al. (2019) can be interpreted as GNNs. In doing so, a difference to Sannche--Gonzaler et al. (2020) and our work is that Ummennhofrr et al.) (2019) use static particles as special nodes the first message passing step only. Consequently, the framework of oSSannhez-Gonzalee et al. (2020), which is based on |Battaglia et al. (018), ). appears to be the most general to us. performing well even without explicit suggested in DPI-Net (Li et al. 2018). Experiments of [Sanchez-Gonzale jet al. (2200) further suggest that their simulation of sand particles are superior to the implementation al. (2020) only consider simple cuboid boundaries for their 3D simulations, leaving more realistic complex geometries as an open and yet untouched challenge. Furthermore, they use sampled, static particles represent boundaries for 2D simulations, which in general does not scale well for 3D simulations due to the quadratic increase of boundary particles (square areas instead of lines).
hierarchical clustering
of Ummenhoter et al. (2019). However, Sanchez-Gonzalez
EXPERIMENTS
We test the effectiveness of BGNNs on complex 3D granular flow simulations. The development, design, and construction of many mechanical devices is based on granular flow simulations. These devices can have very different geometries and must be designed for a wide range of materials with highly varying properties. For example, cohesion properties can range from dry, wet, to oily. In the simulations, we consider very common device geometries and different cohesion properties cover a wide range of situations with our available computational resources. The two common geometries are hoppers and rotating drums (see Fig. [] Fig. and Fig. 5) The two different cohesion properties are non-cohesive describing liquid-like, oily materials and cohesive describing dry, sand-like materials. We compare the BGNN predictions to the simulations in two aspects: speed and accuracy.
Simulation Details. For all experiments, gravitation acts along the z-direction. The upper part of the hopper is delimited along the y-axis by two planes, which are parallel to the x-z plane (see Fig. 4.. The x-axis is delimited by two planes, that are inclined at certain angles a, 180 the

---

Under review as a conference paper at ICLR 2022
X-y plane and at corresponding angles a - 9000000 9000 to the y-2 plane. The hopper has an initially closed hole at the bottom, which has an adjustable radius. The rotation axis of the drum is the y-axis (see Fig.55) The initial filling of the hopper and drum is done by randomly inserting particles into a predefined region. More information can be found in App. C We use around 1000 and around 3000 partieles for hopper and rotating drum simulations, respectively. In order to have trajectories with non-cohesive and cohesive particles, we use the simplified JKR model (Roessler & Katterfeld 2019) with a cohesion energy density of 0 J/m and 10 J/m² for non-cohesive and cohesive particles. The training data consists of 30 simulation trajectories, where each trajectory consists of 100.000 (250.000) simulation timesteps for hopper (rotating drum). For BGNN training every 40 (100)-th timestep is used. Trajectories have different angles and different hole radii (hopper) and different initial particle placement (drum). Moreover, the number of particles is varied by +25%.
Implementation Details. We use 5 message passing layers, with 128 and 512 nodes for interme- diate node and edge representation. The cut-off radii strongly depend on the particle size. We use cut-off radii of 0.02 and 0.008 for rotating drum and hopper, respectively. Cut-off radii have been treated as hyperparameter of our model. More details can be found in App.
Assessment Of Physical Quantities. Granular flow simulations should correctly describe systems on macroscopic scales in terms of particle-averaged positions R(1) and particle flows v()) for particles as a function of time: x(1) = , Σ, (1) and v(r) = J Σ, va( Hoppers are devices that aim at adjusting the flow of particles along the direction of gravity, which coincides with the z-axis in our experiments. Rotating drums are commonly utilized as mixing devices for various applications in c.g. industry, research, and agriculture. They are essentially rotating cylinders that are partially filled with a granular material. The mixing property of these devices is result of numerous particle interactions under time-varying boundary conditions. For rotating drum experiments, we quantify the extend of particle mixing via the mixing entropy (Lai & Fan 1975). If the z-coordinate of particle's ' initial position (0) is above (below) the median z-coordinate of all particles in the initial state, we assign it class c +1 (-1). Based on this assignment local entropies (ggimm 1) at grid cells 9kim are calculated, where the indices klm identify an individual grid cell. The local entropies s(gkim. f) are computed from particle counts ne (gkin. r). of the respective classes c 11. The total number of particles in a grid cell is obtained by (gghmm 1) = 141(gh.m- r) + nー (g r). Calculating the particle-nume number weighted average of the local mixing entropies yields the mixing entropy S(1) of the entire system:
S(1) = Σ n(ghtm.
ΣΣ
n(gkims 1) (fe(gklms 1) log fc (ghtm::)))) 1)),
(8)
1)
klm
kim c=+1
where fo (gkim, 1) denotes the relative fraction of class e particles in cell gklm at time r.
Results. In Fig. and Fig- resuuts for the hopper and the rotating drum simulations are presented. The left parts visualize granular flow snapshots at different time steps. both for cohesive and non- cohesive matcrials. The right parts of the figures include average position and particle flow plots for hopper, well as particle flow and mixing entropy plots for rotating drum simulations. The simulation uncertainties arise due to the different distributions of the initial filling and due to a 25% variation the number of particles across simulations. Short video clips of the trajectories are added as supplement. The difference between cohesive and non-cohesive particles is evident.
BGNNs have learned model granular flow simulations over thousands of time steps. Most notably, hardly any particle leaves the geometric boundaries. This is achieved without using handcrafted conditions or restrictions on the positions of the particles. Furthermore, BGNNs have learned to model barticle-boundary interactions and in doing so correctly represent the dynamics within the system. The predicted quantities are within uncertainties of the simulations. Therefore, we consider the BGNN predictions as sufficiently precise to substitute the simulations. Figure shows of distribution (OOD) scenarios, where the devices are changed with respect to the training data. The hole size of the hopper is decreased in mean by 50%, while side wall inclination angles have been increased by ~ 15a. For the drum the length of the corresponding cylinder was increased in mean by 50%. Our experiments show that our model generalizes well across variations in the geometry. This finding demonstrates that trained BGNNs could be used in the design process of a device to study variations of geometries without retraining the model.
7

---

Under review as a conference paper at ICLR 2022
Table 2 gives run-time comparison of the LIGGGHTS simulation versus forward pass of BGNNs, which only predict every 100 time steps. The highly optimized CPU algorithm (LIGGGHTSS and a non-optimized GPU compatible algorithm (BGNNs) compared via their wall-clock times since the hardware settings are quite different. Nevertheless, Tab. 2 shows that the wall-clock time of BGNNs is shorter than the wall-clock time of the simulation. The usage of more particles, would further increase the lead of BGNNs over the simulation in terms of wall-clock time. For the time comparison, we use typical simulation trajectory from our datasets with 3,408 particles, which needs approximately 2 GB GPU memory for one forward pass. There is potentially even more space for improvement of the BGNN predictions over simulations due to the so called Young's modulus. For simulations, it is often assumed that energy is purely transmitted through Rayleigh waves. Thus the time step DEM simulations is targeted be a fraction of the propagation time through a single, solid particle. As such the propagation time depends on material parameters, most notably the Young's 's modulus. However, for several materials the Young's moduli that reflect the true material properties, would lead to extremely small propagation times, which in turn means much more simulation steps. Consequently, much smaller Young's ' moduli are considered as an approximation, which is valid for gravity driven flows (Coetzee 2017). However, for many cases, e.g. the penetration of a particle bed by an object, this approximation breaks down (Lommen et al.) 2014). BGNNs have the potential be trained on very small time steps reflecting the true Young's moduli and consequently generalize over much more than "just" 40 or 100 time steps.
Table 2: Runtime comparison for one granular flow process consisting of 250.000 simulation timesteps, which are 2500 BGNN predictions.
method
device
specification
time steps
wall-clock time [s]
LIGGGHTS
CPU
AMD EPYC 7H12
356
BGNNs
GPU
NVIDIA A100
100
158
hupper cohesive
hopper non-cohesive
Prodiction
prediction
ground truth
prediction
ground truth
Position chylc
Ground Truth
0.0
02
114
16
6.
Time
Time
Flow zhyit
Prediction
Ground Truth
1A
Time
Figure 4: Hopper dynamics. Left: Distributions for cohesive and non-cohesive particles. Simulation data and BGNN predictions are compared. Particles are indicated by green spheres, triangular wall areas are yellow, the edges of these triangles are indicated by grey lines. In contrast to liquid-like non-cohesive particles, cohesive particles lead to congestion the hopper. Right: Position (upper right) and flow profile (lower right) for non-cohesive particles. Corresponding plots for cohesive particles can be found in App. compared. Simulation uncertainties different initial conditions. To support the reviewing process, we provide simulation predictions for a hopper with more timesteps in animations at attttps: tattppp://bbgndddee / m. bi titbucket...
C Simulation data (solid lines) and BGNN predictions (dashed lines) due to change of the particle numbers (+:25% and to
8

---

Under review as a conference paper at ICLR 2022
rotating drum cohesive prediction
rotating drum non-cohesive prediction
ground truth
ground truth
Flow six
Prediction
Ground Trath
at
0z
Pasition x
Time
Mixing Entropy
Prediction
Ground Truth
Time
Figure 5: Rotating drum dynamics. Left: Particle distributions for cohesive and non-cohesive particles. Simulation data and BGNN predictions are compared. Particles are indicated by green spheres, triangular wall areas are yellow, the edges of these triangles are indicated by grey lines. The circular arrow indicates the rotation direction of the drum. contrast to liquid-like non-cohesive particles, cohesive particles stick together much stronger. Right: position (upper right) and entropy plot (lower right) for non-cohesive particles. The entropy is shown for particle class assignment according to the (blue) and (red) position. Corresponding plots for cohesive particles can be found in App. C Simulation data (solid lines) and BGNN predictions (dashed lines) are compared. Simulation uncertainties conditions. To support the reviewing process, we provide simulation predictions for a rotating drum with more timesteps animations hulhttps: htttps://bgnndsimm... www.ssimmmiibbuu...e 1c/
due change of the particle numbers (+25% and to different initial
7 CONCLUSION AND FUTURE DIRECTIONS
We have introduced Boundary Graph Neural Networks (BGNNs) in order achieve an ac- curate neural network modeling of simulated physical processes with complex geometries. BGNNs dynamically modify graph structures modifying edges, augmenting node features, and dynamically inserting virtual nodes. We have tested BGNNs on complex 3D granular flow processes of hoppers and rotating drums, where BGNNs are able to accurately reproduce these flows within simulation uncertainties over hundreds of thousands of timesteps. Most no- tably particles stay within the geometric objects without using handcrafted conditions or restric- tions. However, should be mentioned that successful simulations often require precise hy- perparameter tuning. far we have not inves- tigated materials with high Young's moduli, and thus we have not yet tested the full generaliza- tion properties of BGNNs for such scenarios described in Sec. 6 Another interesting exten-
0OD hopper non-cohesive prediction
non-cohesive ground truth
ground
prediction
Time
Figure OOD generalization behavior for the hop- per (left) and the rotating drum (right). In contrast to the training and validation data the outlet size of the hopper was decreased, the inclination angles of the hopper side walls are enlarged, and, the length of the rotating drum is increased.
sion of our work might be introduce velocity dependent cut-off radius, and in doing to also consider those interactions which are going to happen within the next timesteps although the spatial distance for particles and respective boundaries is still large.

---

Under review as a conference paper at ICLR 2022
8 REPRODUCIBILITY STATEMENT
Nearly all data used in this work is generated by ourselves. The generalization properties strongly depend on the variability of the training samples since we used a limited number of training samples due to computational reasons. We made sure to keep track of the used parameters to be able to reproduce these datasets. To support the reviewing process. we provide simulation predictions for hopper with more timesteps in animations at
bi tbucket....... .ioeeettto//.
We have included error bars and uncertainty estimates wherever found it necessary and appropriate. For example, for the hopper and drum example, we considered simulation uncertainties. We selected hyperparameters based on a separate validation set. Besides noise terms for the sack regularization, one might fairly similar models if applying the same training and hyperparameter selection procedure. We have described our architecture and implementation details in Appendix [C] We have further provided supporting concepts and experiments in the appendix. For reproducibility, we plan to provide our code upon acceptance.
For reproducibility, we plan to provide our code upon acceptance.
ETHICAL STATEMENT
BGNNs might serve as a valuable tool to obtain speedups for the simulation of industrial processes. However, models obtained, may strongly depend on the variability of the training data. Although we performed some OOD experiments, that showed that our models might still perform well under certain geometric aberrations, it quite clear that a single model will be able to properly predict the flow dynamics of any arbitrary conceivable granular flow setups without ensuring that this behaviour might have been successfully derived from the limited training data. Therefore, the employment of BGNNs needs careful validation wrt. the specific application in mind.
REFERENCES
J.L. Ba, J.R. Kiros, and G.E. Hinton. Layer normalization. arXiv preprint
2016.
P.W. Battaglia, J.B. Hamrick, V. Bapst, A. Sanchez-Gonzalee V.F. Zambaldi, M. Malinowski, A. Tacchetti, D. Raposo, A. Santoro, R. Faulkner, Ç. Gülçehre, H.F. Song, A.J. Ballard, J. Gilmer, G.E. Dahl, A. Vaswani, K.R. Allen, C. Nash, V. Langston, C. Dyer, N. Heess, D. Wierstra, P. Kohli, M. Botvinick, O. Vinyals, Y. Li, and R. Pascanu. Relational inductive biases, deep learning, and graph networks. arXiv preprint 1806v/8806.0/220
2018.
N. Bonneel, M. van de Panne, S. Paris, and W. Heidrich. Displacement interpolation using lagrangian mass transport. 30(6):1-12 December 2011. ISSN 0730-0301. doi: 10.1 10.1145/2070781.2000 2024112.
C.J. Coetzee. Calibration of the discrete element method. Powder Technology. 310:104-142, 2017.
P.A. Cundall and O.D.L. Strack. A discrete numerical model for granular assemblies. Geotechnique, 29(1):47-65, 1979.
M. Defferrard, X. Bresson, and P. Vandergheynst. Convolutional neural networks on graphs with fast localized spectral filtering. In D. Lee, M. Sugiyama, U. Luxburg. I. Guyon, and R. Garnett (eds.), Advances in Neural Information Processing Systems, volume 29, Curran Associates, Inc., 2016.
D. Eberly. Distance between point and triangle in 3d. Retrieved from http://www.magcc- hwt.ware.com/Documnnti oel/varece//occmmntuuiiiiniiisio 1999.
G. Faccanoni and A. Mangeney. Exact solution for granular flows. International Journal for Numerical and Analytical Methods in Geomechanics, 17810);140-1133,, 2013.
R. Flamary and N. Courty. Pot python optimal transport library, 2017. URL htttps://pythone github.ion 1c7
J. Gilmer, S.S. Schoenholz, P.F. Riley, O. Vinyals, and G.E. Dahl. Neural message passing for quantum chemistry. In International Conference on Machine Learning, pp. 1263-1272. PMLR, 2017.
10

---

Under review as a conference paper at ICLR 2022
R.A. Gingold and J.J. Monaghan. Smoothed particle hydrodynamics: theory and application to non-spherical stars. Monthly Notices of the Royal Astronomical Sociery, 181(3):375 389, 12 1977. ISSN 0035-8711. doi: 10. 10.10993/mnrs/188...333 htthttss:///aoioooo 1093/
mnras/181 1..11....77 3.375
K. He, X. Zhang, S. Ren, and J. Sun. Delving deep into rectifiers: Surpassing human-level perfor- mance on imagenet classification. In Proceedings of the 2015 IEEE International Conference on Computer Vision (ICCV), pp. 1026-1034, 2015.
K. He, X. Zhang, S. Ren, and J. Sun. Deep residual learning for image recognition. In Proceedings of the IEEE Conference Computer Vision and Pattern Recognition ( (CVPR), June 2016.
W. Kahan. Pracniques: further remarks on reducing truncation errors. Communications of the ACM, 8(1 8(1):40, 1965.
D.P. Kingma and J.L. Ba. Adam: method for stochastic gradient descent. pp. 1-15, 2015.
T. Kipf and M. Welling. Semi-supervised classification with graph convolutional networks. In International Conference Learning Representations (ICLR), 2017.
Kipf, E. Fetaya, K.C. Wang, M. Welling, and R. Zemel. Neural relational inference for interacting systems, In International Conference on Machine Learning, pp. 2688-2697. PMLR, 2018.
A. Klein. A generalized
Computing, 76(3-4):279-299 2006.
C. Kloss, C. Goniva, A. Hager, S. Amberger, and S. Pirker. Models, algorithms and validation for opensource dem and cfd-dem. Progress in Computational Fluid Dynamics, An International Journal, 12(2/3):140, 2012. ISSN 1468-4349.
C. Kloss, Aigner, Mayrhofer, and C. Goniva. fastdem: method for faster dem simulations of granular media. In Particles 2017,2017.
A. Krizhevsky, I. Sutskever, and G.E. Hinton. Imagenet classification with deep convolutional neural networks. Advances in Neural Information Processing Systems, 25:1097-1105, 2012.
I. Ladický, S. Jeong, B. Solenthaler, M. Pollefeys, and M. Gross. Data-driven fluid simulations using regression forests. ACM Transactions on Graphics (TOG), 34(6):1-9, 2015.
F.S. Lai and L.T. Fan. Application of discrete mixing model the study of mixing of multicompo- nent solid particles. Industrial & Engineering Chemistry Process Design and Development, 14(4): 403-411, 1975. doi: 10. 1...02///260056a000
Y. Li, J. Wu, R. Tedrake, J.B. Tenenbaum, and A. Torralba. Learning particle dynamics for manipu- lating rigid bodies, deformable objects, and fluids. arXiv preprint arXiv: 1.10000.00 2018.
Lommen, D. Schott, and G. Lodewijks. Dem speedup: Stiffness effects on behavior of bulk material. Particuology, 12:107-112, 2014.
L.B. Lucy. numerical approach to the testing of the fission hypothesis. Astronomical Journal, 82: 1013-1024, December 1977. doi: 10.1086/1 12164.
M. Macklin and M. Müller. Position based fluids. ACM Transactions on Graphics (TOG), 32(4): 1-12, 2013.
J. Mellmann, K.I. Iroba, T. Metzger, E. Tsotsas, C. Mészáros, and I. Farkas. Moisture content and residence time distributions in mixed-flow grain dryers. Biosystems Engineering. 109(4):297-307, 2011.
V. Nair and G.E. Hinton. Rectified linear units improve restricted boltzmann machines. In Interna- tional Conference on Machine Learning, pp. 807-814, 2010.
T. Pfaff, M. Fortunato, A. Sanchez-Gonzalez, and P.W. Battaglia. Learning mesh-brr based simulation with graph networks, 2020.
11

---

Under review as a conference paper at ICLR 2022
S. Plimpton. Fast parallel algorithms for short-range molecular dynamics. Journal of computational physics, 117(1):1-19, 1995.
T. Roessler and A. Katterfeld. Dem parameter calibration of cohesive bulk materials using a simple angle of repose test. Particuology. 45:105-1155 2019.
2019.
A. Sanchez-Gonzaler J. Godwin, Pfaff, Ying. J. Leskovec, and P. Battaglia. Learning to simulate complex physics with graph networks. In Proceedings of the 37th International Conference on Machine Learning, volume 119, pp. 8459-8468, 2020.
V.G. Satorras, E. Hoogeboom, and M. Welling. (n) equivariant graph neural networks. arXiv preprint wrXiv:2/02.099844 2021.
Scarselli, M. Gori, A.C. Tsoi, M. Hagenbuchner, and G. Monfardini. The graph neural network model. IEEE Transactions on Neural Networks, 20(1):61-80, 2009.
C. Schenck and D. Fox. Spnets: Differentiable fluid dynamics for deep neural networks. In Conference on Robot Learning, pp. 317-335. PMLR, 2018.
B. Ummenhofer, L. Prantl, N. Thuerey, and Koltun. Lagrangian fluid simulation with continuous convolutions. In International Conference Learning Representations, 2019.
12

---

Under review as a conference paper at ICLR 2022
A DISCRETE ELEMENT METHOD (DEM) SIMULATOR LIGGGHTS
The open source DEM software LIGGGHTS (Kloss et al. 2012) is based on the Molecular Dynamics code LAMMPS (Plimpton| 1995 developed by Sandia National Labs. Due to the similarity of the underlying algorithms for neighbor list construction, output and parallelism this provided a stable basis for the contact models required for DEM. LIGGGHTS added support for triangular mesh walls, particle insertion and new particle shapes (multispheres and superquadrics). Several of those changes resulted in upstream contributions in LAMMPS.
Over the years LIGGGHTS has become a widely used software in both academia and industry that supports both cutting edge research and industrial applications. Support for several physical phenomena as, e.g. liquid transfer on particles, was instrumental in its success. However, it also highlighted requirements for additional research. In industrial applications there often is the need to study physical phenomena which occur on different time scales, e.g. particle collisions (O (100ss)) vs. moisture content particles (1s),, Mellmann et al., 2011) ). which can lead to weeks of simulation
time. While advances have been made to overcome such issues (e.g. (Kloss et al. (0177)) they remain limited in their application, due to the fact that they rely on prior simulation of the exact setup and cannot be used for interpolation of quantities directly related to the flow behavior.
In 2019 LIGGGHTS was again forked and forms the basis of the commercial DEM software Aspherix@ which expands the capabilities of LIGGGHTS with polyhedral particles, a significantly simplified input language and graphical user interface.
B NORMAL VECTOR REPRESENTATIONS
There is ambiguity in the representation of planes via normal vectors due the two possible orientations of the normal vectors, which correspond to the same geometry. In general, there are two possibilities of including normal vector information into the model: (i) encoding always that triangle plane normal vector which always points towards or away from the corresponding particle, (ii) including positively and negatively oriented versions of the normal vector, and order them. We decided for option (ii) since pathological cases where the particle is in the same plane as the triangle are avoided and training is further stabilized. In doing so, we have deal with the fact that the network predictions should be invariant with respect to the orientation of the normal vectors. Therefore, we define a partial ordering which is able to sort the normal vectors with respect to their orientations. For given normal vector (n,,12,
n2, n13) R'', we use the following partial order function
fo(n)= = Σ 3 (sgn (n:) 1)
(B.1)
1=1
to retrieve the scalar values (n) and fo(nn) and sort the two vectors according to their correspond-
ing mapped values. Different sign combinations of the normal vectors are shown in App. Fig. B.] To test the performance of our approach, we conduct a toy experiment as well as a simulation experiment with different representations of normal vectors. We describe both experiments in App. B. and App.
B.2
B.1 REFLECTION TOY EXAMPLE
We conduct a toy experiment to showcase that a partial ordering of normal vectors is helpful for learning 3D simulations. In detail, we consider reflection (Re at plane n as given by
v n
Re (v) = --
n.
(B.2)
n
and try to learn the reflection formula by a simple ReI.U network, which takes the components of and a as input features and predicts the 3 components of Ref (v). The training data consists of reflections at four fixed walls: the top, the bottom, the left, and, the right side of simple cube.
We use normal vectors of these walls, that point towards the inner of the cube. When evaluating the performance of the trained models, we observe decent predictions, if the orientation of the normal vectors describing the inclination of the walls was equal to the training data (see R1-R4 in Fig. B.2). However, for inverted normal vectors in the test set, only networks which take a partial ordering of the normal vectors into account predict the reflection correctly (see R3, R4 in Fig. B.3
13

---

Under review as a conference paper at ICLR 2022
1y 1/11
1/-1
11/0//11
OF 1-11
11/0/11
0/11/11
1/00111
0/0/11
1444
1/-1
1/41/0
0/-1/0
1/0000
+1/1700
+1/11/0
41/0/0
0/4-1/0
11/11-1
HA 41/0-1
0/-1/11
-100-1
0/41711
0/00
Figure B.1: Partial ordering of normal vectors. The numbers indicate the signs of normal vector components, which are used in the partial order function. The figures from left right visualize different (ordered) sign combinations. Sets of sign combinations without zero values form volumes (left), sets with one zero value form planar areas (middle), and sets with sections (right). The 0/0/0 combination forms a point at the origin.
values form line
B.2 SIMULATION EXPERIMENT
We compare three different versions of how to include normal vector information for the hopper particle flow experiments:
not including normal vector information, and filling six node features up with zero entries instead (V1)
including single normal vector orientation, which is given by the triangle corner point order of the mesh (V2)
including both normal vector orientations (six features) (V3).
From an information perspective, it should be noted that (i) distance information (scalar distance and distance vectors) the walls is present in the edge features the graph and (ii) in most cases the used normal vectors are oriented towards the outside of relevant border walls.
The different particle distribution trajectories obtained by the three versions are compared by comput- ing the Earth Movers distances (Bonneel et al.] 2011] Flamary & Courty| (2177, EMD) of predicted and simulated trajectories. We use Euclidean distances for the cost matrix, which we compute at time steps 20, 2,, ... 2, for training trajectories and 5 test trajectories. Table B.1 shows the means (μ) and standard deviations (σr) of EMD values at different time steps and from different training and test trajectories. A paired Wilcoxon test the concatenated trajectories, shows that V3 significantly outperforms V1 (p-value 2.42e-04) and V2 (p-value 1.50e-03) on the test data.
Interestingly, there is less significance the training data, which might indicate that the usage of iitentationiindppende features to represent walls, helps to improve generalization performance, while it might not be that helpful for optimization purposes alone.
Table B.I: Usage of different normal vector information in hopper particle flow experiment. The table summarizes means (Jr) and standard deviations (r)) of the EMD for the different versions and shows the results of paired Wilcoxon test.
Train
Test
Version
p-value Row
p-value
Row V3
VI
No normal vector
5.06e-05
1.17e-04
2.36e-02 3.40e-03
6.80e-05
1.59e-04
2.42e-04 1.50e-03
V2
Single normal vector
1.15e-04
3.84c-04
1.21e-04
4.33e-04
V3
Both orientations
5.99e-05
1.77e-04
6.36e-05
2.06e-04
14

---

Under review as a conference paper at ICLR 2022
ground truth:
wall normal vectors
incident rays
reflected rays
predicted reflected rays
R1:
R2: n. -
R3:
if fo (n) < fo (-n)
R4:
n
if fo (n) & fo (-n)
otherwise
n2
otherwise
Figure B.2: Reflection of rays at four different walls (left, right, bottom, top). Wall normal vectors are visualized by bold orange arrows. The incident rays are visualized by blue arrows, reflected rays by green arrows, and neural network predictions by red arrows. Neural network predictions are based on wall representations that are oriented the same way as in the training phase. The caption above each plot indicates the wall input features used for training each of the networks.
15

---

Under review as a conference paper at ICLR 2022
ground truth:
wall normal vectors
incident rays
reflected rays
predicted reflected rays
R1: n
R2: -n, m
R3:
if fo (-n) & fo (n)
R4:
n.
if fo (-n) & fo (n)
otherwise
n. -2
otherwise
Figure B.3: Reflection of rays four different walls (left, right, bottom, top). Wall normal vectors are visualized by bold orange arrows. The incident rays are visualized by blue arrows, reflected rays by green arrows, neural network predictions by red arrows. Neural network predictions are based on wall representations that are inversely oriented compared to the training phase. The caption above each plot indicates the wall input features used for training each of the networks.
16

---

consider the Under review as a conference paper at ICLR 2022
C EXPERIMENTS
C.I SIMULATION DETAILS
Hopper The initialization consists of two phases. In first step, particles are randomly inserted into a small cuboid which is positioned at certain height above the closed hole of the hopper. This cuboid is continuously filled with particles during the initialization phase and afterwards particles freely move downwards (along the direction of gravity). In this way, the hopper is filled up to certain height with 20,000 particles.
second phase, we cut out particles from the filled mass of particles. We do this (i) by applying randomly selected functions and by (ii) randomly filtering out particles from the whole particle mass. The randomly selected functions are e.g. hyperplanes, where only keep particles if they are at the same side of the hyperplane. The inserted particles have radius 0.002 m.
Drum For initialization we assume that the direction of gravity different than the usual gravitation direction. We insert particles at two random fixed regions within the drum. After the particles are inserted, they can move according the gravitation direction during the initialization phase. this way, we obtain different initial particle distributions within the drum. The inserted particles have radius of 0.01 m.
C.2 IMPLEMENTATION DETAILS
Graph Neural Network Raw inputs to graph networks are initial particle positions and the particle positions from the 5 previous frames of the simulations. From these positions velocities are computed. Further inputs include the particle type and the coordinates of the triangle mesh of the respective time frame. We use residual connections (He et al. (0166) for both node and edge updates. For both updates, we use simple two-layer MLP networks, ReLUs (Nair Hinton 2010) after the first layer, and layer normalization (Ba et al. 2016) without an additional activation after the second layer. For layer normalization we consider the e-parameter as a hyperparameter and set The networks for input embedding and read-out are similar to the message passing layers without layer normalization. The network weights are initialized similar to He et al. (2015); for the input embeddings assume increased number of input neurons for fan in. where additional neurons as virtual copies of c.g. the wall indication feature order to be able to upweight the influence of these features. We use the mean-squared error objective and train with Adam optimization (Kingma & Ba 2013). In order to facilitate learning, we provide as hyperparameter options not only | & Z; ④ as features to the network, but also
1.0.
2, and
T ||
|
2 reflecting the inverse distance law and the inverse-square law, which are present
| ||
many physical laws. We normalize input and target vectors and use variant of Kahan summation
Kahan dataset.
1965
Klein 2006) in order to compute numerically stable statistics across particles of our
Hyperparameter Selection We keep 5 trajectories for each setting aside for validation. Criterions for hyperparameter selection are (i) that particles stay within the geometric object, and (ii) that the ground truth trajectory is reproduced.
C.3 EXPERIMENTAL RESULTS FOR COHESIVE MATERIAL.
In the following, experimental results for a cohesive material are shown. The results in the main paper obtained for non-cohesive granular material, i.e. material with a cohesion energy density of 0 J//3, which results in liquid-like behaviour. Increasing the cohesion energy density to 10 3/m² corresponds to cohesive granular material, i.e. the particles have strong tendency to clump together. Figure C. shows the corresponding comparison of physical quantities for the cohesive granular material. Like in the non-cohesive case, the predictions for cohesive granular material widely in agreement with the ground truth simulation.
17

---

Under review as a conference paper at ICLR 2022
Prediction
Position adyfz
Ground Truth
Flow
Prediction
Ground Truth
0.0
02
Position
Time
Roymi
Mixing Egtropy $
xtyle
Prediction
a.1
Ground Truth
Prediction
Ground Truth
Time
Figure C.1: Experimental results for cohesive granular meterial. Left, top: Time integrated particles flow in rotating drum in the and z direction function of the position along the axis. Left, bottom: Mixing entropies in the rotating drum as function of time for particle class assignments according to the x (blue) and (red) position. Right, top and bottom: Average particle position and particle flows for the hopper as function of time.
18

---

Under review as a conference paper at ICLR 2022
C.4 A WORD ON HOPPER OOD EXPERIMENTS
For hopper geometries (see e.g. Fig. OOD experiments are characterized by an increase of the side wall inclination angles and by an decrease of the radii of the outlet sizes. Especially due to the latter, expect fewer particles to hit the ground for OOD architectures if particle-partied particle and particle-bounday interactions are correctly modeled. In order to statistically test OOD trajectories against in-distribution trajectories, we consider the proportion of particles which have traversed through the outlet of the hopper. We therefore create 15 in-distribution and 15 OOD trajectories for both cohesive and non-cohesive materials. We then apply a Kann-Whitney U test which assesses the proportion of in-distribution against the proportion of OOD particles traversing the outlet. The null hypotheses that the same or higher proportion of particles traverses the outlet for the case of an OOD trajectory compared to an in-distribution trajectory.
Table C.I: Comparison of the proportion (mean μ and std σ)) of particles beyond the outlet of the hopper
cohesive μ
non-cohesive
domain
in in-distribution
0.34
0.09
0.89
0.03
OOD
0.14
0.11
0.73
0.15
The Mann- Woin-Whitney U test shows that the predicted proportion values are significantly lower for OOD than for in-distribution trajectories (p-value < 1.4 10-4 for the cohesive material and p-value 1.2 * 10-+ for the non-cohesive material). We remind the reader that this is expected due to the on average reduced outlet size in OOD geometries. Furthermore, values of the cohesive and the non-cohesive model under the null hypothesis that the same or a higher proportion of particles traverses the outlet for the case of cohesive trajectories compared to non-cohesive trajectories. The applied Mann- WhinnWy U test yields a p-value 1.70 10 06 for the alternate hypothesis that the proportion is lower for the cohesive model, which is also in agreement with rational arguments (cohesive particles tend to clump together) and observations.
compare the predicted proportion
C.5 MIXER
Additionally to hopper and rotating drum geometries, we apply BGNNs to the geometry of rotating mixers. Interactions within the mixer geometry especially challenging, since large number of particles affected indirectly by the blades at the border and especially those in the middle of the mixer. Figure C.2 shows mixer dynamics The left parts visualize granular flow snapshots different time steps, where BGNN predictions and ground truth data are contrasted. The right parts of the figure include flow profile and an entropy curve plot.
C.6 ABLATION STUDIES
For ablating BGNNs, three design choices verified:
1. Is our sampling procedure more effective than dense sampling of the triangularized boundary surface areas?
2. Is a unidirectional particle-wai wall interaction sufficient to learn the corresponding particle-wad wall dynamics?
3. Is the mean node aggregation the right choice for the message passing node updates as outlined in Eq. (2).
In Figs. C.3
C.6 these three points are assessed for cohesive and non-cohesive particles in hopper and rotating drum geometries, respectively. Consistent results across both geometries and particle types are: (i) dense sampling of the boundary surface areas did not learn particle wall interactions well, and, (ii) max node aggregation did not seem to work. For including bidirectional edges instead of unidirectional ones, the behaviour was completely consistent.
C.4
C.5
19

---

Under review as a conference paper at ICLR 2022
mixer
prediction
ground truth
Fore
Protiimn
Time
Miking Entogs
Mising Patray
Entray Toreetti
Figure C.2: Mixer dynamics. Left: Particle distributions for particles in mixer. Simulation data and BGNN predictions contrasted. Particles are indicated by green spheres, triangular wall areas yellow, the edges these triangles are indicated by grey lines. Right: Flow profile (upper right) and entropy plot (lower right) for the used particles. The flow profiles are average velocities of particles a given z coordinate and are shown for the (blue)and the z (red) coordinate. The mixing entropies obtained by splitting particles into two partitions according to threshold the respective coordinate. The prediction entropy curve is shown in yellow while the ground truth simulation entropy curve is shown in blue.
A reason, why we have not been able to successfully train models for dense sampling may be an increased number particle-wall interactions, which might make more complicated to successfully train such model, especially if the number of training trajectories is limited.
For unidirectional particle- wall edges, we could successfully learn a model for cohesive particles of the hopper and for both particle types of the rotating drum; the ablation experiment however not successful for later stages of non-cohesive particles in the hopper. From methodological point of view it should be noted that the reason why skipped bidirectional edges that unidirectional particle-wai wall interactions are simpler and easier to train. In any case, the ablation experiments did not show ad advantages over default model architecture, but seemed confirm choice.
Throughout all ablation experiments max-node aggregation did not work well compared to mean aggregation.
20

---

Under review as a conference paper at ICLR 2022
Cohesive Hopper ablations
boundary edge
sampling
max-age
ground
bidirectionality
Time
Figure C.3: Ablation experiments for cohesive particles within hopper geometries. Different ablations from our default architecture (first three columns, see text) are compared the ground truth (last column). Particles are indicated by green spheres, triangular wall areas are yellow, the edges of these triangles are indicated by grey lines.
21

---

Under review as a conference paper at ICLR 2022
Kon-Cohesivve Hopper Ablations
boundary edge
sampling
max-age
ground
bidirectionality
Time
Figure C.4: Ablation experiments for non non-cohesive particles within hopper geometries. Different ablations from our default architecture (first three columns, see text) are compared to the ground truth (last column). Particles are indicated by green spheres, triangular wall areas are yellow, the edges of these triangles are indicated by grey lines.
22

---

Under review as a conference paper at ICLR 2022
Cohesive Rotating Drum ublations
boundary edge bidirectionality
sampling
max- agg
ground truth
Time
Figure C.5: Ablation experiments for cohesive particles within rotating drum geometries. Different ablations from our default architecture (first three columns, see text) are compared to the ground truth (last column). Particles are indicated by green spheres, triangular wall areas are yellow, the edges of these triangles are indicated by grey lines. The circular arrow indicates the rotation direction of the drum.
23

---

Under review as a conference paper at ICLR 2022
Non-Cohesive Rotating Drum ablations
boundary edge bidirectionality
sampling
max- agg
ground truth
Time
Figure C.6: Ablation experiments for non-cohesive partieles withing rotating drum geometries. Different ablations from our default architecture (first three columns, see text) are compared to the ground truth (last column). Particles are indicated by green spheres, triangular wall areas are yellow, the edges of these triangles are indicated by grey lines. The circular arrow indicates the rotation direction of the drum.
24


Reviews:
[{"review_id": "bHFcveKMZ4", "reviewer": "Reviewer_93ce", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "Strengths:\n- The paper is well written and fairly easy to follow and their is sufficient detail in the supplementary material to suggest that reproducibility is not an issue.\n- The simulation results and their presentation (both the discrete element method used as ground truth and the simulation of the boundary graph GNN) are well done.\n- With regard to the simulation community this method could have practical value since it doesn't require analytical modeling. How easily it could be adapted to other types of machinery though is not entirely clear. The claim that the heuristic method matches ground truth leaves me feeling a bit uneasy though, given the lack of a theoretical analysis or model. Perhaps this is just intractable for granular flows.\n\nWeaknesses:\n- A significant weakness in my opinion is the somewhat incremental methodological contribution. Whereas I can accept the practical value of boundary modifications for the application in mind, the idea of adding virtual nodes to a simulation based on a heuristic ball radius, and the consideration of additional node features (as detailed in Section 3) is modest as a contribution to either a \"representation\" or the \"learning\" of it, which is the general requirement of ICLR.\n- The message passing to update node embeddings seems to be directly based on past work. The specifics of what these messages contain (the new node features, i.e., beyond triangle inclination) are difficult to tease out. Computational considerations, such as how additional coordinate locations are chosen to minimize distance between boundary points and \"real\" nodes are not addressed. There's an attempt to illustrate the idea in Fig. 2 but it falls short of an algorithm and lacks a complexity analysis. \n- Much of Section 4, which discusses further details of how the graph is made \"dynamic\" appears to be a coverage of implementation issues. Here the additional features features are listed: 1) real or virtual particle, 2) normal vectors but how these features are related to the physics of granular flow is not discussed at any length.\n- I found the title to be a bit misleading. While the paper agues for a dynamic modification of a GNN at its boundary for 3D simulations, the experiments and the particular instantiation appear to be almost entirely for granular flow.", "summary_of_the_paper": "This paper considers dynamic modifications of a GNN, where nodes are inserted, node features are augmented and edges are modified, near the boundary of a triangulated surface mesh, using heuristics. The main intended application area is that of simulating granular flows, which are a challenge to model analytically. Experimental results are shown, comparing the boundary graph GNN outputs with simulations obtained using a compute heavy discrete element method, for two cases of interest in industrial machine design: hoppers and rotating drums. The experimental results include both qualitative illustrations and animations and quantitative comparisons (mixing entropies and flow positions).", "summary_of_the_review": "Overall, whereas the ideas are easy to follow and the illustrations are well done and the experimental results demonstrate a match between boundary graph GNN outputs and compute heavy simulations (Figs 4 - 6) the contribution of this paper in terms of new methods or representations or representation learning are modest with regard to where the bar is, in my view, for ICLR. I do not doubt that the method has practical value, and clearly a lot of work went into organizing the experiments, conducting them, and putting together the results. Perhaps this paper would be a better fit to a simulation conference or one with a more engineering oriented focus. The idea to dynamically alter graph structure - and as far as I could tell there are not too many other methodological innovations - is a bit thin with regard to novel content. I'm also puzzled by the choice of title - the choice \"3D simulations\" appears to be quite general, and it is not adequately supported by the paper's content (the focus on granular flows).", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "T72vJSkfOfO", "reviewer": "Reviewer_kKo7", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "Strengths:\n\n1. The paper proposes a solution to a less studied but important issue. To this end, it is potentially of broad interest in the community. \n\n2. The BGNN seems efficient in that every time it just insert one virtual particle due to its specific mechanism while some other SOTA methods for GNN-based 3D simulation have to calculate/update the positions of several particles. \n\n3. Experimental results demonstrate that BGNN works well on the challenging 3D granular flow processes of hoppers and rotating drums.\n\nWeaknesses:\n\n1. The construction of the BGNN largely follows the framework proposed by Sanchez-Gonzalez et al. (ICML’20), which makes the technical novelty of the work a bit limited.\n\n2. The proposed method is merely tested on two the hoppers and the rotating drums. More diverse geometries should have been considered for experiments to demonstrate the generalisation capability of the BGNN. And there is no comparison with previous methods such as Li et al., 2018, Ummenhofer et al., 2019, Sanchez-Gonzalez et al., 2020 and Pfaff et al.,2020.\n\n3. For Equation 9, the detail of the aggregation is not provided. At least some analysis and ideally an ablation study for a mean or max operation should be presented for this vital step.\n\nIn addition, there also exist some minor issues:\n\n1. The title of the paper is inappropriate as the term '3D simulation' is exremely general. Perhaps replace it with '3D granular flow simulation'.\n\n2. Figure 2, as the illustration of one of the key ideas of the paper, is poorly designed and barely informative. It is expected to illustrate how the proposed BGNN dynamically modifies and enhances the graph structure to include boundaries, rather than just exhibits the static local configuration at one time step.\n", "summary_of_the_paper": "This paper presents BGNN to dynamically model geometric boundaries for high quality 3D simulation, with the particular focus on the granular flow simulation. For a geometry represented as a 3D surface mesh, the BGNN updates the edges and insert virtual nodes to allow an accurate modeling of particle interactions with triangularised boundaries. ", "summary_of_the_review": "Based on the strengths and weaknesses of the paper listed above, I tend to take a slightly negative stance.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "Not applicable"}, "extra_text": {}}, {"review_id": "FQwNm5mYlO1", "reviewer": "Reviewer_bSwm", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "The paper is very interesting and well written. Despite its difficulty, it is easy to understand. The background and related work are quite complete and enough to understand the paper. However, there are some parts of the paper that should be explained better to improve its understandability and reproducibility:\n\n1.- The input to the GNN are node coordinates and features but, which kind of features are we talking about?\n\n2.- As commented in point 2, the GNN is composed of 1 encoder, the core message passing layers and the read-out layers. Authors do not specify the configuration of the 5 core-message passing layers nor encoder and read-out layers.\n\n3.- mij appears in (4) and (9) but authors do not explain its meaning. It could be the message passing of Gilmer et al 2017. However, in this case, the mij is a linear combination of the neighbors embedings that have already introduced in both formulas. It would follow the Gilmer proposal if it was eij instead of mij. Please, this part must be explained better.\n\n4.- In page 7, the authors state that \"We use 5 message passing layers with 128 and 512 nodes for intermediate node and edge representation\". But, what structure does it has?. This is the first time that authors talk about using edges in the GNN model. Which information do the edges have? How are their embedings built?\n\n5.- In the results, the authors comment that they introduce a 25% of variability in the number of particles but do not justify this decision.\n\n6.- A comparison with Sanchez-Gonzalez 2020 would be interesting to show the limits of the improvement introduced. ", "summary_of_the_paper": "This paper presents a novel contribution to the use of Graph Neural Networks (GNN) in the simulation of 3D physical phenomena. This contribution consists of the optimization in the calculation of the interactions of the particles with the contour surfaces that allows to optimize the need to add nodes in these limits. The results show an acceleration of the simulation calculation without affecting the accuracy of the result. ", "summary_of_the_review": "Pros: \n1.- A very interesting contribution to GNN applied to 3D simulation.\n2.- It could be a first step towards using GNNs in complex mechanical and biomechanical simulation.\n3.- The increase in speed without affecting too much the precision is a very interesting contribution.\n\nCons:\n1.- The model must be explained better for easier understanding (check the main review).\n2.- It is not clear the differences between the base GNN model (without the dynamic boundary nodes/edges addition) and the one proposed by Sanchez-Gonzalez 2020. (Maybe because cons 1).\n3.- Lack the comparison with other similar approximations.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "4bAhr_OMXjM", "reviewer": "Reviewer_QRXr", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "Strengths:\n\nS1: The proposed approach sounds like a very plausible way to encode triangularized boundary surfaces as nodes and edges in a graph. \n\nS2: The idea is very simple, and can be combined in a very straightforward way with existing models, without any detrimental effects to them (e.g. models can still learn a variety of dynamics, generalize to large systems, etc.).\n\nWeaknesses:\n\nW1: The technical contribution of this paper is sort of a \"featurization\" contribution, rather than a \"modelling\" contribution, e.g. it focuses on the problem of how to build an input graph that you can give to a GNN for a particle simulation with triangular boundaries. While this is the first paper (that I know off) that focuses on that specifically, some of the previous papers were already using sort of similar strategies. For example, Li et al. were already using a single particle to represent each wall in the container (rather than multiple), also adding edges according to the normal/closest distance between the surface and the particles (which does not present the scalability problems indicated by the authors). Ummenhofer et al. also used normals as boundary node features. Sanchez-Gonzalez et. al, included clipped normal distances to walls as additional node features to account for the cubic container to avoid the need to use large numbers of particles for the container walls. The approach proposed in this paper unifies all of those approaches in a sensible and general way, however, the technical novelty feels perhaps a bit too incremental. For comparison, in those 3 previous baselines, the choice of how to build the input graph is sort of a minor implementation detail compared to other contributions, while in this case, the paper focuses solely on that encoding detail.\n\nW2: The paper does not compare against any of those previous baselines using alternatives ways to encode boundaries (e.g. approaches that model surfaces as dense layers of particles). This is probably fine since the authors are not trying to claim that the model accuracy is better at modeling the dynamics, just that it is more efficient as it can do the same with less edges and nodes for the boundaries (and you don't need to actually run a baseline to believe that less nodes and edges will yield a faster model). However, it does take a bit away from the experimental results from the paper, since it probably indicates many of the previous baselines can also model those domains to similar accuracy, and achieve similar flow, mixing entropy and generalization performance.\n\nW3: Authors claim that the model can be faster than state of the art simulation methods. I believe this is probably true, but I feel it is not very well supported as it is a very small part of the paper. E.g. I am not sure if it is fair to say that BGNNs are non-optimized for GPUs, since neural networks are sort of optimized for GPUs.\n\nOther comments:\n\nO1: The description of the contribution could be a bit more clear. If I understood correctly, the idea is quite simple: For every pair of real particle and boundary triangle that are close, add an extra virtual particle at the closest point of the boundary triangle, including the normal of the boundary as an extra feature. Once you have added those particles, just apply *almost* the same model as in the baselines on teh resulting graph. Currently the description of this takes two sections and almost 3 pages (middle of page 3 to near-top of page 6), and there is quite a bit of redundant info between section 3 and 4. I think the paper could benefit from a simpler and shorter description. Or perhaps I am missing something about the methods, in which case I apologize.\n\nO2: Some of the math seems missing some descriptions and parts. For example I believe m_ij in equation 4 and 9 are not defined anywhere. Also, I believe the paper does not currently indicate which of the node features p_i and x_i are used to build the initial set of node embedded features h_i. I assumed absolute positions x_i are *not* used to compute the embedding and just used to build the relative positional features via subtraction at the edges (so the model is translation equivariant), but I am not sure whether absolute positions are used as features or not is clarified anywhere.\n\nO3: Table 2, it is unclear what the time steps column indicates. At first glance it seems to read that BGNN takes 158s to run 100 timesteps, while LIGGGHTS takes 356s to run 1 timestep. Consider renaming the columns to \"timesteps/per step\", \"step size\" or something like that, and \"full trajectory wall clock time\".", "summary_of_the_paper": "The paper proposes BGNNs, which are an extension to existing GNNs for particle simulation (e.g. Li et al, Ummenhofer et al, Sanchez-Gonzalez et al), aimed at improving the efficiency of the model when particles are near boundaries with complex triangular geometries.\n\nSpecifically, the authors propose to parametrize triangular boundaries such that, for any \"real\" particle/node close to any triangle of the boundary surface, a new \"virtual\" particle/node is added at the intersection between the boundary surface and the normal between the real particle and the surface. Furthermore this \"virtual\" particle/node has additional features such as the normal vector to the surface. The authors then implement the model in two geometries (drum and hopper) for two types of simulation (cohesive and non cohesive particles).\n\nThe main claims given the model are:\n\nC1. This approach to describe boundaries is more efficient than those approaches used in prior work.\n\nC2. The model makes predictions for the chosen domains that match quantities of interest to the engineering problem, for each simulation type, and can generalize to larger systems at test time.\n\nC3. The model can be faster than state of the art hardcoded simulation methods.", "summary_of_the_review": "The paper proposal for how to create an input graph for a particle simulation interacting with complex geometries is plausible. However, considering how large overlap with prior work there is, except for the \"graph encoder\" part of the model, I am not sure if this contribution on its own is large enough to meet the ICLR bar. Specifically, I am not sure the paper will bring many new insights to the average ICLR attendee, beyond seeing it as an implementation detail for existing models. In this sense I am borderline between weak accept and weak reject.\n\nI would be much more inclined towards a stronger accept if the focus was a bit more applied:\n1. Here's an actual application, with realistic data, that engineers really care about, and also here's a specific set of generalization settings that engineers really care about.\n2. Here's how you adapt an existing method (e.g. one of Li et al, Ummenhofer et al, Sanchez-Gonzalez et al, Pfaff et al.) to this problem including changes to the encoder to avoid wasting compute on large regions consisting only of boundaries, which ideally really matters if you want to get good running times in that application.\n3. Here's a detailed study of why engineers in those fields should prefer the learned model over other state of the art non-learned model. Including plausibility of generating a training dataset of the required size, limits of generalization, running time, trade offs, etc.\n\nOf course, it may be that a different venue could be more appropriate for a paper like that. Although as an ML conference attendee, I would personally be interested on the specifics of any evidence of the promise of end-to-end learned models replacing traditional simulators in any domain.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper describes means of incorporating boundary conditions into graph neural networks used for simulation.  Assorted techniques dynamically adjust computations near a boundary.  Reviewers agreed this work is well written but disagreed regarding whether the scope of the contribution (and application, which focuses mainly on granular flow despite the title of the paper) merits publication at a top conference like ICLR.

The authors rebut the claims of limited scope strongly, but the experiments in the paper somewhat belie the claims of broad scope.  Some comparisons are also missing to Sanchez-Gonzalez et al. 2020, which appears to be closely related.

The AC agrees that the scope of the work and contribution here have not met the bar for publication.  Reviewer QRXr has some thoughtful suggestions for ways to improve this work in future submissions, or sharing with an audience that can better appreciate the application-oriented contributions might be a reasonable direction to take.
Decision: Reject