Year: 2021

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.

heading: 1 INTRODUCTION
text: Deep learning (Krizhevsky et al., 2012) dramatically changed scientific fields such as 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 aerodynamics (Li et al., 2018; Ummenhofer et al., 2019; Sanchez-Gonzalez 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 we want to focus on learning practically relevant granular flow simulations. Granular flows are 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 a 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 a 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 a 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 1 exemplarily visualizes triangularized surfaces for two standard components in industrial setups, namely a hopper (left) and a rotating drum (right).
In this work, we try to tackle the problem of accurately modeling granular processes when triangularized boundary surfaces are present. We therefore suggest a simple extension of conventional GNNs: Boundary Graph Neural Networks (BGNNs). BGNNs dynamically modify graph structures
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 effectiveness of BGNNs on 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 our 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 are able to generalize granular flow dynamics over hundreds of thousands of timesteps while having the potential to be considerably faster than state-of-the-art simulation methods.

heading: 2 BACKGROUND
text: Dataset. Due to the fact that there are no governing equations for granular flow like the NavierStokes equations for fluid flow (Faccanoni & Mangeney, 2013), a method similar to molecular dynamics, the Discrete Element 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 et al., 2012, see App. A) is an open-source DEM implementation, which is able 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, we used LIGGGHTS to generate training, validation and test set simulation trajectories for granular flow within different hopper and rotating drum environments.
Time transition model. We build our method upon Sanchez-Gonzalez 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 C to time C + 1 is given by ¤xC+1 = ¤xC + ΔC ¥xC and xC+1 = xC + ΔC ¤xC+1, where x is the particle location, and ¤x the particle velocity. To calculate the time-transition xC+1 the particle acceleration ¥xC is predicted.
Graph Neural Networks. We consider graphs G = (V, E), with nodes E8 ∈ V and edges 48 9 ∈ E, where #-dimensional node features pE8 ∈ R# are attached to each of the nodes. Whether the graph G contains an edge between a pair of nodes (E8 , E 9 ) depends on the distance between the nodes:
48 9 ∈ E ⇐⇒ 3 (E8 , E 9 ) 6 cut-off, (1) where the cut-off radius is usually a hyperparamter of the model. Edges might have "-dimensional edge features a8 9 ∈ R" attached to each edge 48 9 . Graph networks are designed to learn from graph-structured data (Scarselli et al., 2009; Kipf & Welling, 2017; Defferrard et al., 2016; Battaglia et al., 2018). Message passing networks (Gilmer et al., 2017) are a 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 (m8 9 ) and nodes (h8), i.e., the embeddings of a8 9 and pE8 , at edge 48 9 and node E8 via:
m′8 9 = q(h8 ,h 9 ,m8 9 ), h′8 = k ( h8 , 48 9 ∈E m ′ 8 9 ) , (2)
where the aggregation 48 9 ∈E at node E8 in Eq. (2) is across all nodes that are connected to node E8 via an edge 48 9 . Typically represents a mean or max operation. The learnable functions q and k are commonly presented by Multilayer Perceptrons (MLPs). Equation (2) describes in a 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.

heading: 3 BOUNDARY GRAPH NEURAL NETWORKS
text: 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 E8 is associated to a particle with location xE8 , velocity ¤xE8 and acceleration ¥xE8 , which is similar to Sanchez-Gonzalez et al. (2020). BGNNs modify and enhance the graph structure to include boundaries (see Fig. 2). BGNNs dynamically add =̃ virtual nodes Ẽ 9 ∈ Ṽ for boundary regions, iff the corresponding boundary region is within a cut-off radius to any other particle. We augment the set of edges E by boundary edges 4̃8 9 giving an enhanced edge set Ê with 48 9 ∈ E and 4̃8 9 ∈ Ẽ. Analogously to Eq. (1), the existence of particle-particle edges 48 9 and particle-boundary edges 4̃8 9 is determined via:
48 9 ∈ E ⊆ Ê ⇐⇒ 3 (E8 , E 9 ) 6 cut-off4, (3) 4̃8 9 ∈ Ẽ ⊆ Ê ⇐⇒ 3̃ (E8 , Ẽ 9 ) 6 cut-off4̃ . (4)
Note, that the cut-off radii cut-off4 and cut-off4̃ are not necessarily the same, and, 3 : V ×V → R, while 3̃ : V × Ṽ → R, i.e. bidirectional edges are used between real nodes and unidirectional edges are used between real and virtual nodes.
In order to include more information about boundary surfaces into particle-boundary interactions, #̃-dimensional node features that encode information about the inclination of triangles in space are concatenated with the existing node features pE8 ∈ R# . Additionally, coordinate information is used both for existing nodes (X = {xE0 , . . . ,xE=−1 }) as well as for virtual nodes (X̃ = {x̃Ẽ0 , . . . , x̃Ẽ=̃−1 }). For virtual nodes, the additional coordinates x̃Ẽ9 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̂ = {pE0 , . . . ,pE=−1 , p̃Ẽ0 , . . . , p̃Ẽ=̃−1 }, X̂ = {xE0 , . . . ,xE=−1 , x̃Ẽ0 , . . . , x̃Ẽ=̃−1 } , (5)
where p̂8 ∈ R#+#̃ and x̂8 ∈ R3 denote the elements of P̂ and X̂ , respectively. Similarly to above, message passing updates the embeddings of edges (m̂8 9 ) and the embeddings of nodes (ĥ8) via
m̂′8 9 = q̂ ( ĥ8 , ĥ 9 , m̂8 9 ) , ĥ′8 = k̂ ( ĥ8 , 4̂8 9 ∈ Ê m̂ ′ 8 9 ) , (6)
where the aggregation 4̂8 9 ∈ Ê at node E8 in Eq. (6) is across all real or virtual nodes that are connected to E8 via an edge 4̂8 9 . Similar to Gilmer et al. (2017) and Satorras et al. (2021), we make use of
pairwise distances ( x̂8 − x̂ 9 2 and x̂8 − x̂ 9 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 as edge attributes â8 9 , for which an initial edge embedding m̂8 9 is determined via an edge embedding layer. The final node embeddings are used for the predictions via the read-out layers. For aggregation , we use the mean. The main challenges of the implementation of BGNNs are the modification of the graph structure towards a dynamic boundary graph structure to include boundary information via the introduction of modified edges 4̃8 9 , modified node features p̂8 , and coordinates for virtual nodes x̃Ẽ8 . In the following section, we discuss how 4̃8 9 , p̂8 , and x̃Ẽ8 are calculated.

heading: 4 DYNAMIC BOUNDARY GRAPH STRUCTURE
text: A boundary graph structure, which is a dynamic modification of the graph structure is 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 in 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” at 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 particle-boundary interaction that boundary point with the largest contribution is considered. Table 1 shows average numbers of nodes
V , as well as average numbers of boundary edges
Ẽ and the relative increase in edges (ratio of the number of added wall edges to the total number of particle edges). The scalability of BGNNs would suffer if more than one particle per particle-boundary interaction surface was considered.
A dynamic boundary graph structure is obtained in three stages:
• Calculating distances between real particles and triangular boundary surface areas in order to decide if a virtual particle node needs to be inserted into the graph, and subsequently obtain additional edges 4̃8 9
• Obtaining positional coordinates for virtual nodes x̃Ẽ9 as representatives of the relevant fraction of the triangular surface area x̃Ẽ9
• Modifying node features p̃8 by including normal directions to facilitate the learning of geometric relationships between particles and boundary surfaces
Modified edges 4̃8 9 via calculation of distances to boundaries. To obtain modified edges 4̃8 9 , boundary particles are dynamically inserted, but only if a 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 are far apart. To decide whether a 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 D and E ∈ R:
t(D, E) = b + D e0 + E e1 , where D ≥ 0, E ≥ 0, and D + E 6 1, b represents one of the nodes of the triangle, and, e0 and e1 are vectors from b towards the other two nodes (see Fig. 3). The minimal Euclidean squared distance @ of the point p to the triangle is given by the optimization problem:
3 = min D,E
@(D, E) = ‖t(D, E) − p‖2 (7)
s.t. D ≥ 0 , E ≥ 0 , D + E 6 1 . The minimizing arguments D′ and E′ parametrize the closest point t(D′, E′) of the triangle to the point p. As indicated in Fig. 3, seven cases have to be to distinguished: one case (c0) in which t(D′, E′) is located within the (closed) triangle, three cases (c1, c3, c5) in which t(D′, E′) 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(D′, E′) is located on one of two edges (including the triangle corner points as special cases). Whether a virtual particle is inserted is determined by Eq. (4) and the particle-triangle distance 3.
Positional coordinates x̃Ẽ9 for virtual nodes. After having modified the edges 4̃8 9 via the calculation 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 x̃Ẽ9 correspond to t (D′, E′), i.e., the closest points on the triangular boundary surface areas with respect to the interacting real particle.
Modifying node features ?̂8 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 Kipf 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 is able to distinguish the different types of nodes. Second, we adapt the weight initialization of the node feature embedding k̂, 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 particle - 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-independent representation of the normal vectors) of the triangular surface areas (null vectors for real particles).

heading: 5 RELATED WORK
text: There is 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 Sanchez-Gonzalez et al. (2020), Ummenhofer et al. (2019), and, Li et al. (2018), all of which propose methods of learning particle simulations without enforcing constraints. These approaches can be contrasted to works like Ladickỳ et al. (2015) or Schenck & Fox (2018) that utilize strong inductive biases. Ladickỳ et al. (2015) construct features for Random Forest Regression that are influenced by Smooth Particle Hydrodynamics (Gingold & Monaghan, 1977; Lucy, 1977). Schenck & Fox (2018) construct a differentiable fluid dynamics network that is closely related to the Position Based Fluids method (Macklin & Müller, 2013). Importantly, both methods are built on the assumption that the governing equations of the system are known, which is 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 (Pfaff et al., 2020). In contrast, our simulator utilizes the mesh to represent static boundaries in a highly efficient way. We share the opinion of Sanchez-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 Sanchez-Gonzalez et al. (2020) and our work is that Ummenhofer et al. (2019) use static particles as special nodes in the first message passing step only. Consequently, the framework of Sanchez-Gonzalez et al. (2020), which is based on Battaglia et al. (2018), appears to be the most general to us, performing well even without explicit hierarchical clustering as suggested in DPI-Net (Li et al., 2018). Experiments of Sanchez-Gonzalez et al. (2020) further suggest that their simulation of sand particles are superior to the implementation of Ummenhofer et al. (2019). However, Sanchez-Gonzalez et 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 to 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).

heading: 6 EXPERIMENTS
text: 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 to cover a wide range of situations with our available computational resources. The two common geometries are hoppers and rotating drums (see Fig. 1, Fig. 4, 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 I-direction. The upper part of the hopper is delimited along the H-axis by two planes, which are parallel to the G-I plane (see Fig. 4). The G-axis is delimited by two planes, that are inclined at certain angles U, 180◦ − U to the
G-H plane and at corresponding angles U − 90◦, 90◦ − U to the H-I plane. The hopper has an initially closed hole at the bottom, which has an adjustable radius. The rotation axis of the drum is the H-axis (see Fig. 5). 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 particles 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/<3 and 105 J/<3 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 U 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 intermediate 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. C.
Assessment Of Physical Quantities. Granular flow simulations should correctly describe systems on macroscopic scales in terms of particle-averaged positions x̄(C) and particle flows v̄(C) for = particles as a function of time: x̄(C) = 1
= ∑ 8 x8 (C) and v̄(C) = 1= ∑ 8 v8 (C). Hoppers are devices that aim
at adjusting the flow of particles along the direction of gravity, which coincides with the I-axis in our experiments. Rotating drums are commonly utilized as mixing devices for various applications in e.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 a 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 a particle’s initial position x8 (0) is above (below) the median z-coordinate of all particles in the initial state, we assign it to class 2 = +1 (−1). Based on this assignment local entropies B(g:;<, C) at grid cells g:;< are calculated, where the indices :;< identify an individual grid cell. The local entropies B(g:;<, C) are computed from particle counts =2 (g:;<, C), of the respective classes 2 = ±1. The total number of particles in a grid cell is obtained by =(g:;<, C) = =+1 (g:;<, C) + =−1 (g:;<, C). Calculating the particle-number weighted average of the local mixing entropies yields the mixing entropy S(C) of the entire system:
S(C) = −1∑ :;< =(g:;<, C) ∑ :;< ∑ 2=±1 =(g:;<, C) ( 52 (g:;<, C) log 52 (g:;<, C) ) , (8)
where 52 (g:;<, C) denotes the relative fraction of class 2 particles in cell g:;< at time C.
Results. In Fig. 4 and Fig. 5 results 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 noncohesive materials. The right parts of the figures include average position and particle flow plots for hopper, as 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 in 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 to 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 particle-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 6 shows out-ofdistribution (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 ∼ 15°. 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.
Table 2 gives a run-time comparison of the LIGGGHTS simulation versus a forward pass of BGNNs, which only predict every 100 time steps. The highly optimized CPU algorithm (LIGGGHTS) and a non-optimized GPU compatible algorithm (BGNNs) are 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 a 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 of DEM simulations is targeted to 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 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 to 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 [B]
LIGGGHTS CPU AMD EPYC 7H12 1 356 BGNNs GPU NVIDIA A100 100 158
hopper cohesive hopper non-cohesive prediction ground truth prediction ground truth
Ti m
e
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 of 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. C. Simulation data (solid lines) and BGNN predictions (dashed lines) are compared. Simulation uncertainties are due to a change of the particle numbers (±25%) and to different initial conditions. To support the reviewing process, we provide simulation predictions for a hopper with more timesteps in animations at https://bgnn3dsim.bitbucket.io/.

heading: 7 CONCLUSION AND FUTURE DIRECTIONS
text: We have introduced Boundary Graph Neural Networks (BGNNs) in order to achieve an accurate neural network modeling of simulated physical processes with complex geometries. BGNNs dynamically modify graph structures via 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 notably particles stay within the geometric objects without using handcrafted conditions or restrictions. However, it should be mentioned that successful simulations often require precise hyperparameter tuning. So far we have not investigated materials with high Young’s moduli, and thus we have not yet tested the full generalization properties of BGNNs for such scenarios as described in Sec. 6. Another interesting exten-
sion of our work might be to introduce a velocity dependent cut-off radius, and in doing so 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.

heading: 8 REPRODUCIBILITY STATEMENT
text: 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 a hopper with more timesteps in animations at https://bgnn3dsim.bitbucket.io/
We have included error bars and uncertainty estimates wherever we 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 of 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.

heading: 9 ETHICAL STATEMENT
text: 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 is quite clear that a single model will not 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.

heading: A DISCRETE ELEMENT METHOD (DEM) SIMULATOR LIGGGHTS
text: 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(10−5B)) vs. moisture content in particles (O(1B), (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., 2017)), 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.

heading: B NORMAL VECTOR REPRESENTATIONS
text: There is an ambiguity in the representation of planes via normal vectors due to 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 to 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 a given normal vector n = (=1, =2, =3) ∈ R3, we use the following partial order function
5> (n) = 3∑ 8=1 38−1 (sgn (=8) + 1) (B.1)
to retrieve the scalar values 5> (n) and 5> (−n) and sort the two vectors according to their corresponding mapped values. Different sign combinations of the normal vectors are shown in App. Fig. B.1. 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.1 and App. B.2.

heading: B.1 REFLECTION TOY EXAMPLE
text: 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 ('4 5 ) at a plane n as given by
'4 5n (v) = v − 2 v · n n · nn , (B.2)
and try to learn the reflection formula by a simple ReLU network, which takes the 3 components of n and v as input features and predicts the 3 components of '4 5n (v). The training data consists of reflections at four fixed walls: the top, the bottom, the left, and, the right side of a 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).

heading: B.2 SIMULATION EXPERIMENT
text: 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) to the walls is present in the edge features of 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 computing the Earth Movers distances (Bonneel et al., 2011; Flamary & Courty, 2017, EMD) of predicted and simulated trajectories. We use Euclidean distances for the cost matrix, which we compute at time steps 20, 21, · · · , 216 for 5 training trajectories and 5 test trajectories. Table B.1 shows the means (`) and standard deviations (f) of EMD values at different time steps and from 5 different training and test trajectories. A paired Wilcoxon test on 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 on the training data, which might indicate that the usage of orientation-independent features to represent walls, helps to improve generalization performance, while it might not be that helpful for optimization purposes alone.
ground truth: n
ground truth: −n

heading: C EXPERIMENTS

heading: C.1 SIMULATION DETAILS
text: Hopper The initialization consists of two phases. In a first step, particles are randomly inserted into a small cuboid which is positioned at a 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 a certain height with 20,000 particles. In a 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 we only keep particles if they are at the same side of the hyperplane. The inserted particles have a radius of 0.002 m.
Drum For initialization we assume that the direction of gravity is 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 to the gravitation direction during the initialization phase. In this way, we obtain different initial particle distributions within the drum. The inserted particles have a radius of 0.01 m.
C.2 IMPLEMENTATION DETAILS
Graph Neural Network Raw inputs to our 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., 2016) 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 n-parameter as a hyperparameter and set it to 1.0. 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 we assume an increased number of input neurons for fan_in, where we consider the additional neurons as virtual copies of e.g. the wall indication feature in order to be able to upweight the influence of these features. We use the mean-squared error as an objective and train with Adam optimization (Kingma & Ba, 2015). In order to facilitate learning, we provide as hyperparameter options not only
x̂8 − x̂ 9 2, x̂8 − x̂ 9 as features to the network, but also 1‖x̂8−x̂ 9 ‖ , x̂8−x̂ 9‖x̂8−x̂ 9 ‖2 and 1
‖x̂8−x̂ 9 ‖2 , x̂8−x̂ 9 ‖x̂8−x̂ 9 ‖3 reflecting the inverse distance law and the inverse-square law, which are present
in many physical laws. We normalize input and target vectors and use a variant of Kahan summation (Kahan, 1965; Klein, 2006) in order to compute numerically stable statistics across particles of our dataset.
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.

heading: C.3 EXPERIMENTAL RESULTS FOR COHESIVE MATERIAL
text: In the following, experimental results for a cohesive material are shown. The results in the main paper are 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 105 J/<3 corresponds to cohesive granular material, i.e. the particles have a strong tendency to clump together. Figure C.1 shows the corresponding comparison of physical quantities for the cohesive granular material. Like in the non-cohesive case, the predictions for cohesive granular material are widely in agreement with the ground truth simulation.

heading: C.4 A WORD ON HOPPER OOD EXPERIMENTS
text: For hopper geometries (see e.g. Fig. 4), 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, we expect fewer particles to hit the ground for OOD architectures if particle-particle and particle-boundary 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 Mann-Whitney U test which assesses the proportion of in-distribution against the proportion of OOD particles traversing the outlet. The null hypotheses is that the same or a higher proportion of particles traverses the outlet for the case of an OOD trajectory compared to an in-distribution trajectory.
The Mann-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−4 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, we compare the predicted proportion 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-Whitney 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.

heading: C.5 MIXER
text: Additionally to hopper and rotating drum geometries, we apply BGNNs to the geometry of rotating mixers. Interactions within the mixer geometry are especially challenging, since a large number of particles is 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 at different time steps, where BGNN predictions and ground truth data are contrasted. The right parts of the figure include a flow profile and an entropy curve plot.

heading: C.6 ABLATION STUDIES
text: For ablating BGNNs, three design choices are verified:
1. Is our sampling procedure more effective than a dense sampling of the triangularized boundary surface areas?
2. Is a unidirectional particle-wall interaction sufficient to learn the corresponding particle-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.4, C.5, 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 not completely consistent.
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 it more complicated to successfully train such a 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 was however not successful for later stages of non-cohesive particles in the hopper. From a methodological point of view it should be noted that the reason why we skipped bidirectional edges was that unidirectional particle-wall interactions are simpler and easier to train. In any case, the ablation experiments did not show advantages over our default model architecture, but seemed to confirm our choice.
Throughout all ablation experiments max-node aggregation did not work well compared to mean aggregation.
Cohesive Hopper ablations
Non-Cohesive Hopper Ablations
Cohesive Rotating Drum ablations
Non-Cohesive Rotating Drum ablations

Review summaries:
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).

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.  Specifically, 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).  The main claims given the model are:  C1. This approach to describe boundaries is more efficient than those approaches used in prior work.  C2. 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.  C3. The model can be faster than state of the art hardcoded simulation methods.

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.