Graph Neural Network-based surrogate model for granular flows
Accurate simulation of granular flow dynamics is crucial for assessing
various geotechnical risks, including landslides and debris flows. Granular
flows involve a dynamic rearrangement of particles exhibiting complex
transitions from solid-like to fluid-like responses. Traditional continuum and
discrete numerical methods are limited by their computational cost in
simulating large-scale systems. Statistical or machine learning-based models
offer an alternative. Still, they are largely empirical, based on a limited set
of parameters. Due to their permutation-dependent learning, traditional machine
learning-based models require huge training data to generalize. To resolve
these problems, we use a graph neural network, a state-of-the-art machine
learning architecture that learns local interactions. Graphs represent the
state of dynamically changing granular flows and the interaction laws, such as
energy and momentum exchange between grains. We develop a graph neural
network-based simulator (GNS) that takes the current state of granular flow and
predicts the next state using Euler explicit integration by learning the local
interaction laws. We train GNS on different granular trajectories. We then
assess the performance of GNS by predicting granular column collapse. GNS
accurately predicts flow dynamics for column collapses with different aspect
ratios unseen during training. GNS is hundreds of times faster than
high-fidelity numerical simulators. The model also generalizes to domains much
larger than the training data, handling more than twice the number of particles
than it was trained on.
Yongjin Choi, Krishna Kumar
# Graph Neural Network-based surrogate model for granular flows

###### Abstract

Accurate simulation of granular flow dynamics is crucial for assessing various geotechnical risks, including landslides and debris flows. Granular flows involve a dynamic rearrangement of particles exhibiting complex transitions from solid-like to fluid-like responses. Traditional continuum and discrete numerical methods are limited by their computational cost in simulating large-scale systems. Statistical or machine learning-based models offer an alternative. Still, they are largely empirical, based on a limited set of parameters. Due to their permutation-dependent learning, traditional machine learning-based models require huge training data to generalize. To resolve these problems, we use graph neural network, a state-of-the-art machine learning architecture that learns local interactions. Graphs represent the state of dynamically changing granular flows and the interaction laws, such as energy and momentum exchange between grains. We develop a graph neural network-based simulator (GNS) that takes the current state of granular flow and predicts the next state using Euler explicit integration by learning the local interaction laws. We train GNS on different granular trajectories. We then assess the performance of GNS by predicting granular column collapse. GNS accurately predicts flow dynamics for column collapses with different aspect ratios unseen during training. GNS is hundreds of times faster than high-fidelity numerical simulators. The model also generalizes to domains much larger than the training data, handling more than twice the number of particles than it was trained on.

keywords: graph neural network, learned physics simulator, granular column collapse, surrogate model +
Footnote †: journal: Computers and Geotechnics

## 1 Introduction

Landslides cause extensive material displacement and significant infrastructure damage. Accurate modeling of granular flow runout is crucial to understanding the impact of landslides. Numerical methods, such as particle-based and continuum approaches, are often employed to assess landslide runouts. Particle-based approaches, like the Discrete Element Method (DEM) (Staron and Hinch, 2005; Kermani et al., 2015; Kumar et al., 2017), can model grain-grain interactions but are limited to representative elemental volumes. Traditional continuum approaches, such as the Finite Element Method, can predict the initiation of such failures but suffer from mesh distortions when capturing runout dynamics. Hybrid Eulerian-Lagrangian approaches like the Material Point Method (MPM) (Mast et al., 2014; Kumar et al., 2017) can simulate large-deformation flows without undergoing mesh distortions. However, the hybrid nature of MPM requires tracking both the grid and the material
points, which is computationally expensive. Multiple full-scale simulations are necessary for a comprehensive evaluation of runout hazard scenarios. Similarly, a back analysis to estimate material parameters requires a broad parametric sweep involving hundreds to thousands of simulations. However, current state-of-the-art numerical methods are restricted to, at most, a few full-scale simulations, limiting our ability in scenario testing or back analysis.

An alternative to numerical simulations is the development of statistical or machine learning models to evaluate landslide risks. These surrogate models build correlations between landslide risks and their influencing factors through simple empirical correlation without considering the complex granular flow dynamics. Several studies adopt probabilistic approaches, such as Monte Carlo simulation and Bayesian analysis, to evaluate the landslide runout distance based on factors including topology and geology (Gao et al., 2021; Zeng et al., 2021; Sun et al., 2021; Zhao et al., 2022). Machine learning models can predict the travel distance and potential path of granular flows based on the geometry and ground properties (Durante and Rathje, 2021; Ju et al., 2022; Yang and Hambleton, 2021). Although researchers have been able to correlate the runout of granular flow based on statistical or data-driven techniques, these techniques do not explicitly consider granular flow dynamics--the actual physics governing the flow behavior. Thus, due to a lack of physics, these statistical models do not generalize outside their training range in modeling other boundary conditions or geometry.

Building surrogate models that replicate the entire granular flow dynamics is challenging. The surrogate model must capture complex behaviors involving highly non-linear, static, collisional, and frictional dissipation regimes (Soga et al., 2016). Learning fundamental interaction laws is crucial for generalizing beyond the training datasets. Techniques like max-pooling in convolutional neural networks learn spatially invariant behavior, i.e., they learn features irrespective of their spatial location. However, CNNs are primarily limited to mesh-based systems with fixed neighbors.

Granular flow is a dynamic system where neighbor interactions evolve throughout the runout (Lajeunesse et al., 2005; Zhang et al., 2016; Soga et al., 2016). A traditional Multi-Layer Perceptron (MLP) could model such a dynamic system. However, generalizing MLPs requires an exhaustive dataset to overcome combinatorial dependence, i.e., the outputs of the models depend on the order of the inputs (Battaglia et al., 2018; Haeri and Skoniczny, 2022). Unreasonably large training datasets are needed to map the entire parameter space of particle arrangements and dynamics.

To address these limitations, we utilize graph neural networks (GNNs), a state-of-the-art machine learning architecture that enables permutation invariant learning (Battaglia et al., 2016, 2018; Sanchez-Gonzalez et al., 2020), to develop a data-driven surrogate model for granular flow dynamics. At any given time, the physical state of the granular system is represented as a graph. We develop a GNN-based Simulator (GNS) that operates on the graph to learn the fundamental interaction law. We demonstrate the capability of GNS in replicating the granular flow dynamics by studying the collapse of a granular column. Granular column collapse is a simple physical experiment that captures the overall dynamics of large-scale runouts. GNS, trained on granular flow trajectories, successfully predicts the runout dynamics of column collapse outside its training range and generalizes to upscaled domain sizes.
## 2 Method

This section describes the individual components of GNS: graphs, graph neural networks (GNNs), and message passing.

### Graph Neural Networks and Message Passing

#### 2.1.1 Graphs

Graphs can represent interactions in physical systems (Battaglia et al., 2016; Sanchez-Gonzalez et al., 2020). We represent the granular media as a graph \(G=(\mathbf{V},\mathbf{E})\) consisting of a set of vertices (\(\mathbf{v}_{i}\in\mathbf{V}\)) representing the soil grains or aggregation of grains and edges (\(\mathbf{e}_{i,j}\in\mathbf{E}\)) connecting a pair of vertices (\(\mathbf{v}_{i}\) and \(\mathbf{v}_{j}\)) representing the interaction between the grains. Consider an example involving interaction between grains in a box (see fig. 1). We encode the state of the physical system, such as the kinematics of grains and their interaction (fig. 1a and fig. 1d), as a graph (fig. 1b and fig. 1c). The vertices describe the position and velocity of the grains, and the edges describe the directional interaction between them, shown as arrows in fig. 1b and fig. 1c. The state of the grain \(i\) is represented as a vertex feature vector \(\mathbf{v}_{i}\). The vertex feature vector includes velocities, mass, and distance to the boundary. The edge feature vector \(\mathbf{e}_{i,j}\) includes information about the interaction between grains \(i\) and \(j\) such as the relative distance between the grains. Thus, we can store and process the state of granular bodies and their interactions as graphs.

Graphs offer a permutation-invariant form of encoding data, where the interaction between vertices is independent of the order of vertices or their position in Euclidean space. As graphs represent the interactions between grains as edge connections, graphs are permutation invariants. For example, by storing the relative positional information in \(\mathbf{e}_{i,j}\), rather than the absolute position, machine learning models operating on these networks learn the interaction behavior of different relative distances between grains. Therefore, graphs can efficiently represent the physical state of granular flow involving multi-grain interactions.

#### 2.1.2 Graph neural networks (GNNs)

GNN takes a graph \(G=(\mathbf{V},\mathbf{E})\) as an input, computes properties and updates the graph, and outputs an updated graph \(G^{\prime}=(\mathbf{V}^{\prime},\mathbf{E}^{\prime})\) with an identical structure, where \(\mathbf{V}^{\prime}\) and \(\mathbf{E}^{\prime}\) are the set of updated vertex and edge features (\(\mathbf{v}^{\prime}_{i}\) and \(\mathbf{e}^{\prime}_{i,j}\)). Message passing is the process of updating the graph by propagating information through it.

In the grains-in-a-box example, the GNN first takes the original graph \(G=(\mathbf{V},\mathbf{E})\) (fig. 1b) that describes the current state of the physical system (\(\mathbf{X}_{t}\)). The GNN then updates the state of the physical system through message passing, which models the exchange of energy and momentum between the grains, and returns an updated graph \(G^{\prime}=(\mathbf{V}^{\prime},\mathbf{E}^{\prime})\) (fig. 1c). We decode \(G^{\prime}\), the output of GNN, to extract information related to the future state of the physical system (\(\mathbf{X}_{t+1}\)), such as the next position or acceleration of the grains (fig. 1d).

#### 2.1.3 Message passing

Message passing consists of three operations: message construction (eq. (1)), message aggregation (eq. (2)), and the vertex update function (eq. (3)).

\[\mathbf{e}^{\prime}_{i,j}=\phi_{\mathbf{\Theta}_{\phi}}\left(\mathbf{v}_{i},\mathbf{v}_{j}, \mathbf{e}_{i,\ j}\right) \tag{1}\]
\[\bar{\mathbf{v}}_{i}=\Sigma_{j\in N(i)}\ \mathbf{e}^{\prime}_{i,j} \tag{2}\]

\[\mathbf{v}^{\prime}_{i}=\gamma_{\mathbf{\Theta}_{\gamma}}\left(\mathbf{v}_{i},\bar{\mathbf{v}}_{ i}\right) \tag{3}\]

The subscript \(\mathbf{\Theta}_{\phi}\) and \(\mathbf{\Theta}_{\gamma}\) represent a set of learnable parameters in each computation. The message construction function \(\phi_{\mathbf{\Theta}_{\phi}}\) (eq. (1)) takes the feature vectors of the receiver and sender vertices (\(\mathbf{v}_{i}\) and \(\mathbf{v}_{j}\)) and the feature vector of the edge connecting them (\(\mathbf{e}_{i,\ j}\)) and returns an updated edge feature vector \(\mathbf{e}^{\prime}_{i,j}\) as the output. \(\phi_{\mathbf{\Theta}_{\phi}}\) is a matrix operation including the learnable parameter \(\mathbf{\Theta}_{\phi}\). The updated edge feature vector \(\mathbf{e}^{\prime}_{i,j}\) is the message sent from vertex \(j\) to \(i\). Figure 1(a) shows an example of constructing messages on edges directed to vertex \(0\) originating from vertices \(1\), \(2\), and \(3\) (\(\mathbf{e}^{\prime}_{0,\ 1}\), \(\mathbf{e}^{\prime}_{0,\ 2}\), \(\mathbf{e}^{\prime}_{0,\ 3}\)). Here, we define the message construction function \(\phi_{\mathbf{\Theta}_{\phi}}\) as \(((\mathbf{v}_{i}+\mathbf{v}_{j})\times\mathbf{e}_{i,\ j})\times\mathbf{\Theta}_{\phi}\). The updated feature vector \(\mathbf{e}^{\prime}_{0,\ 1}\) is computed as \(((\mathbf{v}_{0}+\mathbf{v}_{1})\times\mathbf{e}_{0,\ 1})\times\mathbf{\Theta}_{\phi}\), where \(\mathbf{v}_{0}\) and \(\mathbf{v}_{1}\) are the receiver and sender vertex feature vectors, and \(\mathbf{e}_{0,\ 1}\) is their edge feature vector. Suppose we assume all values of \(\mathbf{\Theta}_{\phi}\) are \(1.0\) for simplicity, we obtain \(\mathbf{e}^{\prime}_{0,\ 1}=(([1,\ 0,\ 2])+[1,\ 3,\ 2])\times[2,\ 1,\ 0]^{T})\times 1=[4,\ 3,\ 0]\). Similarly, we can compute the messages \(\mathbf{e}^{\prime}_{0,\ 2}=[0,\ 3,\ 9]\) and \(\mathbf{e}^{\prime}_{0,\ 3}=[3,\ 4,\ 9]\).

The next step in message passing is the message aggregation \(\Sigma_{j\in N(i)}\) (eq. (2)), where \(N(i)\) is the set of sender vertices \(j\) related to vertex \(i\)). It collects all the messages directing to vertex \(i\) and aggregates those into a single vector with the same dimension as the aggregated message (\(\bar{\mathbf{v}}_{i}\)). The aggregation rule can be element-wise vector summation or averaging; hence it is a permutation invariant computation. In fig. 1(a), the aggregated message \(\bar{\mathbf{v}}_{0}=[7,\ 10,\ 18]\) is the element-wise summation of the messages directing to vertex \(0\) as \(\bar{\mathbf{v}}_{0}=\mathbf{e}^{\prime}_{0,\ 1}+\ \mathbf{e}^{\prime}_{0,\ 2}+\ \mathbf{e}^{ \prime}_{0,\ 3}\).

The final step of the message passing is updating vertex features using eq. (3). It takes the aggregated message (\(\bar{\mathbf{v}}_{i}\)) and the current vertex feature vector \(\mathbf{v}_{i}\), and returns an updated vertex feature vector \(\mathbf{v}^{\prime}_{i}\), using predefined vector operations including the learnable parameter \(\mathbf{\Theta}_{\gamma}\). Figure 1(b) shows an example of the update at vertex \(0\). Here, we define the update function \(\gamma_{\mathbf{\Theta}_{\gamma}}\) as \(\mathbf{\Theta}_{\gamma}\left(\mathbf{v}_{i}+\bar{\mathbf{v}}_{i}\right)\). The updated feature vector \(\mathbf{v}^{\prime}_{0}\) is computed as \(\mathbf{\Theta}_{\gamma}\left(\mathbf{v}_{0}+\bar{\mathbf{v}}_{0}\right)\).

Figure 1: An example of a graph and graph neural network (GNN) that process the graph (modified from Battaglia et al. (2018)): (a) A state of the current physical system (\(\mathbf{X}_{t}\)) where the grains are bouncing in a box boundary; (b) Graph representation of the physical system (\(G\)). There are three vertices representing grains and six edges representing their directional interaction shown as arrows; (c) The updated graph (\(G^{\prime}\)) that GNN outputs through message passing; (d) The predicted future state of the physical system (\(\mathbf{X}_{t+1}\)) (i.e., the positions of the grains at the next timestep) decoded from the updated graph.
Assuming all parameters in \(\mathbf{\Theta}_{\gamma}\) are \(1.0\) for simplicity, we obtain \(\mathbf{v}_{0}^{\prime}=\mathbf{v}_{0}+\bar{\mathbf{v}}_{0}=[1,\ 0,\ 2]+[7,\ 10,\ 18]=[8,\ 10,\ 20]\). Similarly, we update the other vertex features \((\mathbf{v}_{1}^{\prime},\ \mathbf{v}_{2}^{\prime},\ \mathbf{v}_{3}^{\prime})\).

After message passing, the graph vertex and edge features \((\mathbf{v}_{i}\) and \(\mathbf{e}_{i,\ j})\) are updated to \(\mathbf{v}_{i}^{\prime}\) and \(\mathbf{e}_{i,\ j}^{\prime}\). The GNN may include multiple message passing steps to propagate the information further through the network.

Unlike the example shown above, where we assume a constant value of \(1.0\) for the learnable parameters, in a supervised learning environment, the optimization algorithm will find a set of the best learnable parameters \((\mathbf{\Theta}_{\phi},\mathbf{\Theta}_{\gamma})\) in the message passing operation.

### Graph Neural Network-based Simulator (GNS)

In this study, we use GNN as a surrogate simulator to model granular flow behavior. Figure 3 shows an overview of the general concepts and structure of the GNN-based simulator (GNS) proposed by Sanchez-Gonzalez et al. (2020). Consider a granular flow domain represented as material points (fig. 3a), which represent the collection of grains. In GNS, we represent the physical state of the granular domain at time \(t\) with a set of \(\mathbf{x}_{i}^{t}\) describing the state and properties of each material point. The GNS takes the current state of the granular flow \(\mathbf{x}_{t}^{i}\in\mathbf{X}_{t}\) and predicts its next state \(\mathbf{x}_{i+1}^{i}\in\mathbf{X}_{t+1}\) (fig. 3a). The GNS consists of two components: a parameterized function approximator \(d_{\Theta}\) and an updater function (fig. 3b). The function approximator \(d_{\Theta}\) takes \(\mathbf{X}_{t}\) as an input and outputs dynamics information \(\mathbf{y}_{i}^{t}\in\mathbf{Y}_{t}\). The updater then computes \(\mathbf{X}_{t+1}\) using \(\mathbf{Y}_{t}\) and \(\mathbf{X}_{t}\). Figure 3c shows the details of \(d_{\Theta}\) which consists of an encoder, a processor, and a decoder. The encoder (fig. 3-c1) takes the state of the system \(\mathbf{X}^{t}\) and embeds it into a latent graph \(G_{0}=(\mathbf{V}_{0},\ \mathbf{E}_{0})\) to represent the relationships between material points. The vertices \(\mathbf{v}_{i}^{t}\in\mathbf{V}_{0}\) contain latent information of the current state of the material point, and the edges \(\mathbf{e}_{i,j}^{t}\in\mathbf{E}_{0}\) contain latent information of the pair-wise relationship between material points. Next, the processer (fig. 3-c2) converts the input graph \(G_{0}\) to the output graphs \(G_{M}\) through \(M\) stacks of message-passing GNN (\(G_{0}\to G_{1}\to\cdots\to G_{M}\)). The message passing computes the interaction between vertices. Finally, the decoder (fig. 3-c3) extracts the dynamics of the points \((\mathbf{Y}^{t})\) from \(G_{M}\), such as the acceleration of the physical system. The entire simulation (fig. 3a) involves running GNS surrogate model through \(K\) timesteps predicting from the initial state
to \(\mathbf{X}_{K}\left(\mathbf{X}_{0},\ \mathbf{X}_{1},\ \dots,\ \mathbf{X}_{K}\right)\), updating at each step (\(\mathbf{X}_{t}\rightarrow\mathbf{X}_{t+1}\)). We call this successive prediction from GNS the "rollout".

In the following sections, we explain the details of our input \(\mathbf{X}^{t}\) (fig. 3a), the encoder, processor, and decoder in \(d_{\Theta}\) (fig. 3c), and how we compute \(\mathbf{X}^{t+1}\) from \(\mathbf{X}^{t}\) using the GNS updater function (fig. 3b).

#### 2.2.1 Input

The input to the GNS, \(\mathbf{x}_{i}^{t}\in\mathbf{X}^{t}\) (eq. (4)), is a vector consisting of the current material point position \(\mathbf{p}_{i}^{t}\), the material point velocity context \(\mathbf{\hat{p}}_{i}^{\leq t}\), information on boundaries \(\mathbf{b}_{i}^{t}\), and material point type embedding \(\mathbf{f}\). The current state \(\mathbf{x}_{i}^{t}\) will be used to construct vertex feature (\(\mathbf{v}_{i}^{t}\)) (eq. (6)).

\[\mathbf{x}_{i}^{t}=\left[\mathbf{p}_{i}^{t},\ \mathbf{p}_{i}^{\leq t},\ \mathbf{b}_{i}^{t},\ \mathbf{f}\right] \tag{4}\]

The velocity context \(\mathbf{\hat{p}}_{i}^{\leq t}\) includes the current and previous material point velocities for \(n\) timesteps \(\left[\mathbf{\hat{p}}_{i}^{t-n},\cdots,\ \mathbf{\hat{p}}_{i}^{t}\right]\) with \(n+1\) velocities. We use \(n=4\) to include sufficient velocity context in the vertex feature \(\mathbf{v}_{i}^{t}\). Sanchez-Gonzalez et al. (2020) show that having \(n>1\) significantly improves the model performance. We compute the velocities using the finite difference of the position sequence (i.e., \(\mathbf{\hat{p}}_{i}^{t}=\left(\mathbf{p}_{i}^{t}-\mathbf{p}_{i}^{t-1}\right)/\Delta t\)). \(\mathbf{b}_{i}^{t}\) is boundary information. For a 2D problem, \(\mathbf{b}_{i}^{t}\) has four components, each indicating the distance between material points and the four walls. We normalize \(\mathbf{b}_{i}^{t}\) by the connectivity radius \(R\) which defines the interaction zone, explained in the next section, and restrict it between -1.0 to 1.0. \(\mathbf{b}_{i}^{t}\) is used to evaluate boundary interaction for a material point. \(\mathbf{f}\) is a vector embedding describing a material point type.

We define the interaction between material points \(i\) and \(j\) as \(\mathbf{r}_{i,\ j}^{t}\) using the distance and displacement of the material points in the current timestep (see eq. (4)). The former reflects

Figure 3: The structure of the graph neural network (GNN)-based physics simulator (GNS) for granular flow (modified from Sanchez-Gonzalez et al. (2020)): (a) The entire simulation procedure using the GNS, (b) The computation procedure of GNS and its composition, (c) The computation procedure of the parameterized function approximator \(d_{\Theta}\) and its composition.
the level of interaction, and the latter reflects its spatial direction. \(\mathbf{r}_{i,\ j}^{t}\) will be used to construct edge features (\(\mathbf{e}_{i,j}^{t}\)).

\[\mathbf{r}_{i,j}^{t}=\left[\left(\mathbf{p}_{i}^{t}-\mathbf{p}_{j}^{t}\right),\ \|\mathbf{p}_{i}^{t}-\mathbf{p}_{j}^{t}\|\right] \tag{5}\]

#### 2.2.2 Encoder

The vertex and edge encoders (\(\varepsilon_{\mathbf{\Theta}}^{v}\) and \(\varepsilon_{\mathbf{\Theta}}^{e}\)) convert \(\mathbf{x}_{i}^{t}\) and \(\mathbf{r}_{i,\ j}^{t}\) into the vertex and edge feature vectors (\(\mathbf{v}_{i}^{t}\) and \(\mathbf{e}_{i,j}^{t}\)) (eq. (6)) and embed them into a latent graph \(G_{0}=(\mathbf{V}_{0},\ \mathbf{E}_{0})\), \(\mathbf{v}_{i}^{t}\in\ \mathbf{V}_{0}\), \(\mathbf{e}_{i,j}^{t}\in\ \mathbf{E}_{0}\).

\[\mathbf{v}_{i}^{t}=\varepsilon_{\mathbf{\Theta}}^{v}\left(\mathbf{x}_{i}^{t}\right),\ \mathbf{e}_{i,j}^{t}=\varepsilon_{\mathbf{\Theta}}^{e}\left(\mathbf{r}_{i,j}^{t}\right) \tag{6}\]

We use a two-layered 128-dimensional multi-layer perceptron (MLP) for the \(\varepsilon_{\mathbf{\Theta}}^{v}\) and \(\varepsilon_{\mathbf{\Theta}}^{e}\). The MLP and optimization algorithm search for the best candidate for the parameter set \(\mathbf{\Theta}\) that estimates a proper way of representing the physical state of the material points and their relationship which will be embedded into \(G_{0}\).

The edge encoder \(\varepsilon_{\mathbf{\Theta}}^{v}\) uses \(\mathbf{x}_{i}^{t}\) (eq. (4)) without the current position of the material point (\(\mathbf{p}_{i}^{t}\)), but with its velocities (\(\mathbf{\dot{\mathbf{p}}}_{i}^{\pm t}\)), as velocity governs the momentum, and the interaction dynamics is independent of the absolute position of the material points. Rubanova et al. (2022) confirmed that including position causes poorer model performance. We only use \(\mathbf{p}_{i}^{t}\) to predict the next position \(\mathbf{p}_{i}^{t+1}\) based on the predicted velocity \(\mathbf{\dot{\mathbf{p}}}_{i}^{t+1}\) using Explicit Euler integration.

We consider the interaction between two material points by constructing edges between them all pairs of vertices located within a certain distance called connectivity radius \(R\) (see the shaded circular area in fig. 3b). The connectivity radius is a critical hyperparameter that governs how effectively the model learns the local interaction. \(R\) should be sufficiently large to include the local interaction between material points and capture the simulation domain's global dynamics.

#### 2.2.3 Processor

The processor performs message passing (based on eq. (1) to eq. (3)) on the initial latent graph (\(G_{0}\)) from the encoder for \(M\) times (\(G_{0}\to G_{1}\rightarrow\cdots\to G_{M}\)) and returns a final updated graph \(G_{M}\). We use two-layered 128-dimensional MLPs for both the message construction function \(\phi_{\mathbf{\Theta}_{\phi}}\) and vertex update function \(\gamma_{\mathbf{\Theta}_{r}}\), and element-wise summation for the message aggregation function \(\Sigma_{j\in N(i)}\) in eq. (1) to eq. (3). We set \(M=10\) to ensure sufficient message propagation through the network. These stacks of message passing model information propagation through the network of material points.

#### 2.2.4 Decoder

The decoder \(\delta_{\mathbf{\Theta}}^{v}\) extracts the dynamics \(\mathbf{y}_{i}^{t}\in\mathbf{Y}^{t}\) of the material points from the vertices \(\mathbf{v}_{i}^{t\prime}\) (eq. (7)) using the final graph \(G_{M}\). We use a two-layered 128-dimensional MLP for \(\delta_{\mathbf{\Theta}}^{v}\), which learns to extract the relevant dynamics for material points from \(G_{M}\).

\[\mathbf{y}_{i}^{t}=\delta_{\mathbf{\Theta}}^{v}\left(\mathbf{v}_{i}^{t\prime}\right) \tag{7}\]
#### 2.2.5 Update

We use the dynamics \(\mathbf{y}_{i}^{t}\) to predict the velocity and position of the material points at the next timestep (\(\mathbf{\dot{p}}_{i}^{t+1}\) and \(\mathbf{p}_{i}^{t+1}\)) based on Euler integration (eq. (8) and eq. (9)), which makes \(\mathbf{y}_{i}^{t}\) analogous to acceleration \(\mathbf{\ddot{p}}_{i}^{t}\).

\[\mathbf{\dot{p}}_{i}^{t+1}=\mathbf{\dot{p}}_{i}^{t}+\mathbf{y}_{i}^{t}\Delta\mathrm{t} \tag{8}\]

\[\mathbf{p}_{i}^{t+1}=\mathbf{p}_{i}^{t}+\mathbf{\dot{p}}_{i}^{t+1}\Delta\mathrm{t} \tag{9}\]

Based on the new position and velocity of the material points, we update \(\mathbf{x}_{i}^{t}\in\mathbf{X}^{t}\) (eq. (4)) to \(\mathbf{x}_{i}^{t+1}\in\mathbf{X}^{t+1}\). The updated physical state \(\mathbf{X}^{t+1}\) is then used to predict the position and velocity for the next timestep.

The updater imposes inductive biases, such as an inertial frame, on GNS to force it only to learn the interaction dynamics, improving learning efficiency. A traditional neural network learns both the update scheme and the interaction dynamics:

\[p^{t+1}=NN(p^{t},v^{t})\,. \tag{10}\]

Whereas, using an inertial prior, we force the GNS only to learn the interaction dynamics, by hardcoding the update function:

\[p^{t+1}=p^{t}+v^{t}\cdot\Delta t+NN(p^{t},v^{t})\,. \tag{11}\]

Nevertheless, GNS does not directly predict the next position from the current position and velocity (i.e., \(\mathbf{p}_{i}^{t+1}=GNS\left(\mathbf{p}_{i}^{t},\ \mathbf{\dot{p}}_{i}^{t}\right)\)) which has to learn the static motion and inertial motion. Instead, it uses (1) the inertial prior (eq. (8)) where the prediction of the next velocity \(\mathbf{\dot{p}}_{i}^{t+1}\) should be based on the current velocity \(\mathbf{\dot{p}}_{i}^{t}\) and (2) the static prior (eq. (9)) where the prediction of the next position \(\mathbf{p}_{i}^{t+1}\) should be based on the current position \(\mathbf{p}_{i}^{t}\). These make GNS focus on learning unknown dynamics by hardcoding known physics. Since GNS learns the dynamics of material points through interactions independent of absolute position, GNS is generalizable to other geometric conditions.

## 3 Training and Evaluation

We now train the GNS to predict granular column collapse. This section explains how we generate training data, details of the training process, and how we evaluate the performance of the GNS.

### Material Point Method

We utilize the Material Point Method (MPM) to generate the GNS training dataset of granular flow simulations. The MPM is a hybrid Eulerian-Lagrangian approach designed for modeling large-deformation flows (Soga et al., 2016). In the MPM, a continuum body is represented by individual material points that traverse a static background grid. The governing equation is solved at the nodes, and the updated velocity field is subsequently mapped back to the material points. We employ the position information stored in these
material points to construct the current state \(\mathbf{X}^{t}\) in the GNS. For more information on MPM refer to Soga et al. (2016).

### Datasets

The training datasets include 26 granular flow trajectories of square-shaped granular mass in a two-dimensional box boundary simulated by the MPM explicit time integration method using the CB-Geo MPM code (Kumar et al., 2019). Each simulation has a different initial configuration regarding the size of the square granular mass, position, and velocity. Table 1 presents the details of the training dataset generated using the MPM simulation. The datasets are published on DesignSafe (Kumar and Choi, 2023). A shows all the training trajectories with different initial configurations and initial velocities.

We also create the validation datasets to check if the model experiences overfitting. The datasets include seven trajectories of randomly picked rectangular-shaped granular mass with different initial configurations not included in the training datasets.

### Training

Our GNS has a learnable parameter set \(\Theta\). We train \(\Theta\) to minimize the loss calculated as the mean squared error (MSE) between \(\mathbf{y}^{t}_{i}\) (predicted proxy-acceleration) and the ground truth acceleration \(\mathbf{\vec{p}}^{t}_{i}\) for all material points \(i=1,\ 2,\ \dots,\ N\) as shown in eq. (12) based on gradient (\(\nabla loss_{\Theta}\))-based optimizer, Adam (Kingma and Ba, 2014).

\[loss_{\Theta}=\frac{1}{n}\sum_{i=1}^{N}\left(\mathbf{y}^{t}_{i}-\mathbf{\vec{p}}^{t}_ {i}\right)^{2} \tag{12}\]

For training the GNS, we have to set hyperparameters to learn the flow behavior from the training trajectories properly. The first key hyperparameter is the connectivity radius \(R\)

\begin{table}
\begin{tabular}{l l l} \hline \hline \multicolumn{2}{c}{Property} & \multicolumn{1}{c}{Description} \\ \hline Simulation boundary & \multicolumn{1}{c}{1.0\(\times\)1.0 m} \\ Mesh size & \multicolumn{1}{c}{0.025\(\times\)0.025 m} \\ Material points per cell & \multicolumn{1}{c}{16} \\ Granular mass geometry & \multicolumn{1}{c}{0.2\(\times\)0.2 m and 0.3\(\times\)0.3 m} \\ Simulation duration (timesteps) & \multicolumn{1}{c}{(each includes 1,024 and 2,304 particles)} \\ \hline \multirow{6}{*}{Material property} & Model & Mohr-Coulomb \\  & Density & 1,800 \(kg/m^{3}\) \\  & Youngs modulus & 2 GPa \\  & Poisson ratio & 0.3 \\  & Friction angle & 30\({}^{\circ}\) \\  & Cohesion & 0.1 kPa \\  & Tension cutoff & 0.05 kPa \\ \hline \hline \end{tabular}
\end{table}
Table 1: Details of the Material Point Method (MPM) simulation geomaterials and properties used for generating the training datasets.
which governs the model's capacity to learn the interactions of material points. We set \(R=0.030\) m which includes about 9 to 10 material points along the diameter. The circular area defined by \(R\) can incorporate approximately 70 material points inside. Another important hyperparameter is the Gaussian noise value for perturbing the ground truth position in the training trajectories. Since every predicted position for each timestep is based on the previous prediction, which includes a prediction error, the simulation over the large timesteps is subjected to an exponential error accumulation. To avoid this issue, we train the model on input positions with Gaussian noise that emulates the prediction error made by a one-step prediction (\(\mathbf{X}_{t}\rightarrow\mathbf{X}_{t+1}\)). The inclusion of noise in training leads to more rigorous long-rollout predictions.

We use the learning rate (\(lr\)) decay with the initial value of \(10^{-4}\) and decay rate of 0.1 (\(lr=10^{-4}\times 0.1^{step/5\times 10^{6}}\)) for more stable convergence. We use the batch size of two, i.e., \(\mathbf{X}_{t}\) from two different trajectories are used simultaneously in updating the learnable parameters. For information on the scalability of the GNS algorithm, refer to Kumar and Vantassel (2022).

We investigate if the model experiences overfitting by plotting the loss history (fig. 4) for the training and validation datasets evaluated for every 10K training steps. The training and validation losses show a drastic decrease until 2M steps. After that, the validation loss tends to remain slightly larger than the training loss. Figure 4 shows no overfitting during the training.

### GNS prediction of granular flows

We trained the GNS to predict the collapse of a granular column (as studied by Lajeunesse et al. (2004); Lube et al. (2005)). Figure 5 shows the granular column collapse experiments to evaluate its ability to replicate granular flow dynamics. Granular column collapse is a simple physical experiment that captures the transient response of granular flow dynamics. The experiment involves the collapse of a granular column of initial height \(H_{0}\) and length \(L_{0}\) on a flat surface due to gravity. As the gate holding the column is removed, the granular

Figure 4: Evolution of GNS loss in training and validation with epochs.
material destabilizes, resulting in a runout. We measure the final runout deposit with the final height \(H_{f}\) and runout \(L_{f}\).

The runout of the granular column is governed by the initial aspect ratio (\(a=H_{0}/L_{0}\)) (Staron and Hinch, 2005; Kumar, 2015). For short columns (\(a\lesssim 2\)) (fig. 5a), the soil mass fails along the flanks of the column above a well-defined failure surface (dashed line). The soil mass beneath the failure surface remains in static equilibrium throughout the collapse forming a truncated conical shape. With the increase in aspect ratio, the portion of the sliding mass above the failure surface increase, and the static part becomes smaller, forming a conical shape. For tall columns (\(a\gtrsim 2\)) (fig. 5b), the majority of the soil mass is involved in the collapse, and it initially experiences a free fall due to gravitational acceleration. As the falling mass reaches the failure surface, the vertical kinetic energy is converted to horizontal acceleration, resulting in a longer runout distance than the short column (fig. 5a). In addition, researchers (Kumar, 2015; Staron and Hinch, 2005; Kermani et al., 2015; Utili et al., 2015) observed a transition zone where the flow dynamics change from short to tall columns. The normalized runout (\(\left(L_{f}-L_{0}\right)/L_{0}\)) of a granular column is only a function of its aspect ratio (\(a\)). The normalized runout represents how far the granular mass runs out before reaching the final deposit state compared to the initial length of the column. Short columns show a linear relationship with the initial aspect ratio. In contrast, tall columns have a power-law relationship with the initial aspect ratio.

The GNS was trained only on the aspect ratio of 1.0. However, we evaluate its performance in predicting the runout dynamics of other aspect ratios by comparing the GNS predictions with the MPM simulations. Table 2 presents the test cases for evaluating GNS performance.

## 4 Results and Discussions

We evaluate the GNS predictions of granular column collapse against the MPM simulations in terms of the (1) geometry of sliding mass, (2) evolution of runout and height with time, and (3) energy evolution during the collapse. Figure 6 shows the normalized runout (\(\left(L_{f}-L_{0}\right)/L_{0}\)) predictions of GNS for different aspect ratios in comparison with MPM. \(L_{f}\) is the distance from the left wall to the material point that runs out the farthest, as shown in fig. 5. Previous research observed a transition zone for the relationship between the normalized runout and aspect ratio that distinguishes short-column from tall-column dynamics.

Figure 5: Schematics for the granular column collapse configuration and its behavior on the aspect ratio.
For both GNS and MPM, we observe the transition around an initial aspect ratio \(a=1.2\) (fig. 6). Table 3 summarizes the errors between GNS predictions and MPM simulations for different aspect ratios. In general, the GNS runout prediction is within 5% of the MPM runout estimate. Figure 6 suggests that the GNS successfully captures the dependence of the final runout with the initial aspect ratio, including the transition from the short to the tall column.

### GNS Predictions of Granular Flow Dynamics

#### 4.1.1 Short Column

We now evaluate the GNS rollout (prediction) of the granular flow dynamics with time for a short column (\(a=0.8\)). Figure 7 shows the time evolution of granular flow for the

\begin{table}
\begin{tabular}{c c c c c} \hline \hline \multirow{2}{*}{Test case} & \multirow{2}{*}{\(H_{0}\times L_{0}\)} & Duration & Simulation & Number of \\  & & & (timesteps) & boundary & material points \\ \hline Short & & & & X: 0 to 1.0 & \\ columns & \(a=0.5\) & 0.2 \(\times\) 0.4 & 400 & \begin{tabular}{c} X: 0 to 1.0 \\ Y: 0 to 0.5 \\ Y: 0 to 1.0 \\ Y: 0 to 0.5 \\ Y: 0 to 0.5 \\ Y: 0 to 0.5 \\ \end{tabular} & 1956 \\  & \(a=0.8\) & \(0.24\times 0.30\) & 400 & \begin{tabular}{c} X: 0 to 1.0 \\ Y: 0 to 1.0 \\ Y: 0 to 0.5 \\ Y: 0 to 0.5 \\ \end{tabular} & 1824 \\  & \(a=1.0\) & 0.30 \(\times\) 0.30 & 400 & \begin{tabular}{c} X: 0 to 1.0 \\ Y: 0 to 0.5 \\ Y: 0 to 0.5 \\ \end{tabular} & 2304 \\ \hline Tall & & & & X: 0 to 1.0 & \\ columns & \(a=2.0\) & 0.30 \(\times\) 0.15 & 400 & \begin{tabular}{c} X: 0 to 1.0 \\ Y: 0 to 0.5 \\ Y: 0 to 1.0 \\ Y: 0 to 0.5 \\ \end{tabular} & 1152 \\  & \(a=3.0\) & 0.36 \(\times\) 0.12 & 400 & \begin{tabular}{c} X: 0 to 1.0 \\ Y: 0 to 0.5 \\ \end{tabular} & 1106 \\  & \(a=4.0\) & 0.35 \(\times\) 0.075 & 400 & \begin{tabular}{c} X: 0 to 1.0 \\ Y: 0 to 0.5 \\ \end{tabular} & 576 \\ \hline Up-scaled & \(a=0.8\) & 0.36 \(\times\) 0.45 & 500 & 
\begin{tabular}{c} X: 0 to 1.5 \\ Y: 0 to 1.0 \\ \end{tabular} & 5120 \\ \hline \hline \end{tabular}
\end{table}
Table 2: Granular column collapse simulation cases for testing GNS.

\begin{table}
\begin{tabular}{c c c c} \hline \hline \multirow{2}{*}{Aspect ratio, \(a\)} & \multicolumn{2}{c}{Normalized runout} & \multirow{2}{*}{Runout error (\%)} \\ \cline{2-3}  & MPM & GNS & \\ \hline
0.5 & 0.831 & 0.811 & 2.48 \\
0.8 & 1.444 & 1.445 & 0.06 \\
1.0 & 2.071 & 2.152 & 3.78 \\
2.0 & 3.892 & 3.682 & 5.70 \\
3.0 & 5.620 & 5.341 & 5.23 \\
4.0 & 5.753 & 6.070 & 5.21 \\ \hline \hline \end{tabular}
\end{table}
Table 3: Normalized runout from MPM and GNS depending on aspect ratios and corresponding prediction error.
short column collapse. We use a normalized time (\(t/\tau_{c}\)) to compare the flow evolution, where \(t\) is physical time, and \(\tau_{c}\) is the critical time defined as the time required for the flow to fully mobilize. \(\tau_{c}\) is defined as \(\sqrt{H_{0}/g}\), where \(g\) is the gravitational acceleration. In fig. 7, the collapse shows three stages. First, the flow is mobilized by the failure of the flank and reaches full mobilization around \(t/\tau_{c}=1.0\). The majority of the runout occurs until \(t/\tau_{c}=2.5\). Beyond \(t/\tau_{c}>2.5\), the spreading decelerates due to the basal friction and finally stops at around \(t/\tau_{c}=4.0\) for both MPM and GNS rollout (prediction). As seen in fig. 7, although the GNS has only seen an aspect ratio \(a=1.0\) during training, GNS successfully captures the overall time-evolution of granular flows for a short column (\(a=0.8\)).

In addition to the visual comparison of profiles, we quantitatively investigate the flow dynamics of the GNS rollout of the short column by comparing the normalized runout and height evolution with the MPM. Figure 8a shows the evolution of normalized runout and height with time. The normalized runout of the MPM (see the gray line corresponding to the left axis in fig. 8a) shows the three stages of collapse. The collapse of the granular column starts with the failure of the flank and evolves slowly until the runout is fully mobilized by \(t/\tau_{c}=1.0\). As the collapse proceeds, the runout acceleration increases (\(t/\tau_{c}=1.0\) to \(2.5\)). After this time, the runout deaccelerates due to basal friction, and finally stops at \(t/\tau_{c}\approx 4.0\). Both GNS and MPM show a constant normalized height (see the gray line corresponding to the right axis in fig. 8a) as only the flank of the column collapse, leaving a static truncated-conical core. GNS predicts an almost identical evolution of runout as the MPM simulation, which is noteworthy as only a small portion of the training trajectories (\(5\) out of \(26\)) includes the deacceleration behavior leading to the flow coming to rest due to the basal friction before hitting the walls. Overall, the quantitative comparison shown in fig. 8a confirms that the GNS can accurately model the granular flow dynamics for the short column.

Figure 8b shows the energy evolutions from GNS rollout and MPM simulation. Based on the principle of energy conservation, the granular flow must satisfy \(E_{0}=E_{p}+E_{k}+E_{d}\)

Figure 6: Normalized runout distance (\(\left(L_{f}-L_{0}\right)/L_{0}\)) with different aspect ratios (\(a\)).
Figure 7: Evolution of flow with normalized time for GNS and MPM for the short column with \(a=0.8\). Units are in \(m\). The color represents the magnitude of the displacement. Subfigure (e) shows the final deposit at the last timestep.
where \(E_{0}\) is the potential energy of the column before material points start to mobilize, \(E_{p}\) is the potential energy, \(E_{k}\) is the kinetic energy, and \(E_{d}\) is the dissipation energy due to friction along the boundary and material. In fig. 7(b), both GNS rollout and MPM simulation show identical energy evolutions. A significant fraction of the stored potential energy is converted to kinetic energy in the initial stages of the failure, reaching a peak value of kinetic energy at \(t/\tau_{c}=1\). The kinetic energy dissipates due to the basal friction and flow ceases at \(t/\tau_{c}=4.0\) when \(E_{k}\) is fully dissipated.

#### 4.1.2 Tall column

Tall columns exhibit different runout dynamics than the short column. GNS was only trained on granular mass with an aspect ratio of 1.0 and has not seen the dynamics of a tall column during training. As an example, we demonstrate the GNS prediction for a tall column with \(a=2.0\). Figure 9 shows the GNS rollout and MPM simulation of the runout evolution for this case. GNS rollout predicts an identical runout profile with normalized time as the MPM simulation. Similar to the short column, the tall column also shows the three stages: the initial mobilization of the flow (\(t/\tau_{c}\) to 1.0), runout (\(t/\tau_{c}=1.0\) to 2.5) along the failure surface, deacceleration (\(t/\tau_{c}=2.5\) to 4.0). In the tall column, however, a larger volume of sliding mass above the failure plane is mobilized. During the initial stages of the collapse, the granular mass experiences free fall due to gravity dominated by collisional dissipation. As the granular mass reaches the failure surface, the vertical kinetic energy is converted to horizontal acceleration, resulting in longer runout distances. GNS rollout shows similar behavior to the MPM runout simulation.

Figure 9(a) shows the normalized runout and height evolution for the tall column. Although the runout evolution remains identical in the initial phase of the collapse, MPM

Figure 8: (a) Normalized runout and height evolution with normalized time and (b) normalized energy evolution with normalized time for the short column \(a=0.8\). \(H_{t}\) is the height from the bottom corner of the boundary to the highest part of the column at \(t\). \(E_{p}=\sum_{i=1}^{n}m_{i}gh_{i}\) is the potential energy of the column, and \(E_{k}=\frac{1}{2}\sum_{i}^{n}m_{i}v_{i}^{2}\) is the kinetic energy of the column, where \(m_{i}\), \(h_{i}\), and \(v_{i}\) is the mass, height, and velocity of the material point \(i\), and \(n\) is the total number of material points. \(E_{d}=E_{0}-E_{p}-E_{k}\) is the dissipation energy where \(E_{0}\) is the potential energy before material points start to move.
Figure 9: Evolution of flow with normalized time for GNS and MPM for the tall column with \(a=2.0\). Units are in \(m\). The color represents the magnitude of the displacement. Subfigure (e) shows the final deposit at the last timestep.
shows a slightly larger normalized runout compared to the GNS. The final height in both GNS and MPM remains the same.

Figure 10 presents the normalized energy evolution of the GNS rollout and the MPM simulation. During the initial stages of the collapse (\(t/\tau_{c}\) to \(1.0\)), a large amount of initial potential energy is converted to kinetic energy due to the free fall of mass under gravity. Both GNS and MPM show almost identical energy profiles. GNS shows a larger potential energy loss as the flow accelerates with an almost similar gain in kinetic energy. It indicates that GNS predicts larger frictional dissipation in tall columns, which could be from the training data focused only on short columns having higher frictional dissipation than tall columns. At the final stage, MPM shows less dissipation due to the basal boundary friction, resulting in a slightly longer runout than GNS rollout. Generally, energy dissipation behavior in GNS is consistent with MPM showing a more significant potential drop and increased dissipation energy accumulation.

Overall, the GNS rollout is consistent with the MPM simulation with a runout error of \(5.7\) % for the tall column with \(a=2.0\), implying that the GNS can capture the dynamics of granular flows in collision-dominated tall columns despite only being trained on \(a=1.0\).

#### 4.1.3 Upscaled domain

GNS is generalizable to different initial configurations of the flow simulation owing to the strong inductive bias of the GNN(Battaglia et al., 2018). The strengths of GNS surrogate models would be to train them on small-scale experiments and then predict large-scale dynamic scenarios with complex boundary conditions. We now evaluate the scalability of GNS to a larger domain, including more material points than the training dataset. Figure 11 shows the GNS rollout of a short column \(a=0.8\) with \(5120\) material points (up to \(5\times\) more material points than the training dataset) for a larger simulation domain and longer rollout duration than the training dataset.

GNS successfully predicts the flow dynamics for an upscaled domain size showing a similar runout profile with the MPM simulation. The GNS rollout predicts a normalized runout of

Figure 10: (a) Normalized runout and height evolution with normalized time and (b) normalized energy evolution with normalized time for the tall column with \(a=2.0\).
1.74 while the MPM simulation shows 1.76, showing an error of 1.30%. Figure 12 shows that GNS rollout successfully replicate energy evolution observed in an upscaled domain compared to the MPM simulation. Hence, GNS can reproduce the flow dynamics even for the upscaled geometries beyond the training dataset.

The primary source of GNS rollout error is not from the simulation scale but from the portion of material points that shows a large amount of displacement during column collapse. Figure 13 shows the evolution of mean squared error (MSE) of displacement over all material points (\(N\)) with time \(t\) computed as \(\frac{1}{n}\sum_{i}^{N}\left(\boldsymbol{p}_{i}^{t}-\boldsymbol{p}_{MPMi}^{t} \right)^{2}\), where \(\boldsymbol{p}_{MPMi}^{t}\) is material point position from MPM. When we compare the MSE for \(a=0.8\) with 1,824 of material points and its upscaled domain (2.22x material points), upscaling does not alter the MSE significantly. Figure 14 shows the evolution of the squared error of displacement of individual material points for the upscaled domain (\(a=0.8\)). The squared error shows larger values for those material points which run out further, i.e., the larger final displacements, but the proportion of the error with respect to the final runout is small so that GNS could simulate the upscaled case without significant error.
Figure 11: Evolution of flow with normalized time for GNS and MPM for the upscaled case of short column with \(a=0.8\). Units are in \(m\). The color represents the magnitude of the displacement. Subfigure (e) shows the final deposit at the last timestep.
Figure 12: Normalized runout and height evolution with normalized time and (b) normalized energy evolution with normalized time for the upscaled case of the short column with \(a=0.8\). Note that the data after \(t/\tau_{c}>5.0\) is abbreviated since the flow reaches a static state.

Figure 13: Evolution of mean squared displacement error over all material points with time.
## 6 Limitations

The GPU memory limits the current implementation of the GNS surrogate model. A GPU with 40 GB memory can simulate up to around 50K material points (approximately 3M edge connections). However, this shortcoming can be improved by optimizing the size of the connectivity radius \(R\). We use \(R\) of 0.030 m, which includes more interaction between neighbors. Multi-GPU GNS rollouts will enable the scalability of GNS to larger and more complex domains.

## 7 Conclusion

Traditional numerical methods are computationally intensive when simulating large-scale granular flows. Statistical or conventional machine learning-based surrogate models are not generalizable since they do not explicitly consider the underlying physics. In this study, we develop a graph neural network (GNN)-based simulator (GNS) as a generalizable granular flow surrogate simulator. The use of graphs efficiently represents the physical state of interacting material points. At the same time, the message passing operation of GNN encourages the neural network to learn the interaction between material points. The expressive power of graphs and message passing that models the interaction between material points allows GNS to accurately predict granular flow dynamics for various conditions, including those not seen during training. We demonstrate the performance of GNS on granular column collapse. GNS precisely simulates different flow dynamics involved in columns for different initial aspect ratios and can also be applied to the upscaled domain with more than 2 to \(5\times\) material points with a longer simulation duration than the data provided for training. GNS also shows a remarkable speed-up of \(150\times\) computation speed compared to the parallelized CPU version of MPM. The computational efficiency and generalizability of the GNS surrogate can expedite evaluating runout hazards requiring numerous scenarios.

Figure 14: Evolution of squared displacement error for each material point with normalized time in upscaled case of \(a=0.8\). The line color represents the final displacement of each material point.