# Towards Complex Dynamic Physics System Simulation with Graph Neural ODEs

Guangsi Shi, Daokun Zhang, Ming Jin and Shirui Pan Guangsi Shi is with the
Department of Chemical and Biological Engineering, Faculty of Engineering,
Monash University, Clayton, VIC 3800, Australia (email:
guangsi.shi@monash.edu).Daokun Zhang and Ming Jin are with the Department of
Data Science and Artificial Intelligence, Faculty of IT, Monash University,
Clayton, VIC 3800, Australia (emails<EMAIL_ADDRESS>and
ming.jin@monash.edu).Shirui Pan is with the School of Information and
Communication Technology, Griffith University, South Port, QLD 4215, Australia
(email: s.pan@griffith.edu.au).
Corresponding authors: Daokun Zhang and Shirui Pan

###### Abstract

The great learning ability of deep learning models facilitates us to
comprehend the real physical world, making learning to simulate complicated
particle systems a promising endeavour. However, the complex laws of the
physical world pose significant challenges to the learning based simulations,
such as the varying spatial dependencies between interacting particles and
varying temporal dependencies between particle system states in different time
stamps, which dominate particles’ interacting behaviour and the physical
systems’ evolution patterns. Existing learning based simulation methods fail
to fully account for the complexities, making them unable to yield
satisfactory simulations. To better comprehend the complex physical laws, this
paper proposes a novel learning based simulation model – Graph Networks with
Spatial-Temporal neural Ordinary Equations (GNSTODE) – that characterizes the
varying spatial and temporal dependencies in particle systems using a united
end-to-end framework. Through training with real-world particle-particle
interaction observations, GNSTODE is able to simulate any possible particle
systems with high precisions. We empirically evaluate GNSTODE’s simulation
performance on two real-world particle systems, Gravity and Coulomb, with
varying levels of spatial and temporal dependencies. The results show that the
proposed GNSTODE yields significantly better simulations than state-of-the-art
learning based simulation methods, which proves that GNSTODE can serve as an
effective solution to particle simulations in real-world applications.

###### Index Terms:

AI for Physics Science, Learning-based Simulator, Graph Neural Networks,
Neural Ordinary Differential Equations.

##  I Introduction

Simulation has become an indispensable methodology in many scientific and
technological disciplines to comprehend the properties, behaviour, and dynamic
changes of real-world matters. Physical world simulations are often built on
the first principles [1], i.e., simulating according to the established
physics laws without any extra assumptions.

Figure 1: The rollout of a particle system simulation.

Although the conventional Euler and Lagrangian simulation methods [2] have
been employed by many disciplines, including pharmaceutical science and
materials science, their simulation ability is still limited by some domain-
specific challenges. First, the conventional physical simulation methods are
usually developed on some prior expert knowledge describing the objective
systems, like the governing first principles, which makes them hard to
generalize to the generic scenes with no prior knowledge. Second, many
physical systems are described by ordinary differential equations (ODEs) or
partial differential equations (PDEs). As such, many trials are required to
adjust the hyper-parameters of the ODE/PDEs to make them well fit real-world
physical systems, which inevitably involves massive computational costs used
to solve the complex ODE/PDEs. Third, numerous simulators can only work for a
single task and knowledge transfer across tasks is hard to achieve to yield
high-quality simulations. On the other hand, the learning-based simulation
breaks away from the shackles of human knowledge and can expand the knowledge
by uncovering more complex laws through learning from real-world observations.
This makes the learning-based simulation a more promising solution to
comprehending the complex physical world.

Some learning based algorithms have been developed to achieve the complex
physical simulations, including fluid simulation [3], solid simulation [4] and
particle simulation [5, 6]. Among others, graphs provide a straightforward
data structure to describe the complex particle systems, with particles
modelled as nodes and the interacting relations between particles described by
edges. Based on the graph modeling of particle systems, Graph Neural Network
(GNN) [7] has been leveraged as the simulation learner, due to its great
expressive power. The GNN based simulation methods [8, 5, 9] leverage message
passing mechanism to characterize the interactions between particles. Though
competitive performance has been achieved in some cases, the applicability of
GNN based simulation methods are still seriously challenged by the inherent
complexities in the real-world particle systems:

  * •

Challenge I: Varying Spatial Dependencies. In real-world physical systems, the
random particle-particle interacting patterns result in the varying spatial
dependencies between interacting particles. Due to the irregular movements of
particles, different particles can be surrounded by a varying number of
neighbors, making particles have different centralities. With varying
centralities, the interactions imposed by neighboring particles should be
weighed differently to calculate the temporal dynamics of central particles at
each time step. For particles with large centralities, interactions from close
neighboring particles are enough to interpret their movement, while distant
neighboring particles should also be considered for particles with small
centralities. How to automatically determine the range of neighboring
particles used to derive the driving forces of central particles is quite
important to ensure accurate simulations.

  * •

Challenge II: Varying Temporal Dependencies. Particle systems evolve with
varying temporal patterns, i.e., the particle systems might change
significantly in a show time slot or remain stable for a long time. Existing
learning based simulation methods mainly choose a fixed time step to roll out
the simulation, by iteratively predicting the next time’s state from the
current time’s state, as is shown in Fig. 1. However, the rigid fixed-time-
step prediction cannot capture the varying temporal dependencies of particle
systems, without considering the possible dramatic changes occurring in the
fixed time step. As a result, the one-step prediction inevitably involves many
simulation errors, which will be further accumulated along the rollout of
prediction and result in the model collapse dilemma [5]. Therefore, the
ability to adapt to the varying temporal dependencies of particle systems is
critical to a learning based simulator for yielding stable and high-quality
simulations.

To effectively address the particle simulation chanllenges caused by varying
spatial and temporal dependencies, we propose a novel learning based
simulation model – Graph Networks with Spatial-Temporal neural Ordinary
Equations (GNSTODE) – that leverages neural ordinary equations to learn
varying spatial correlations between particles with varying distances and to
model the complex system evolving dynamics. Compared with the existing methods
[5, 10, 11] that use the rigid increment-adding operators (e.g., the Runge-
Kutta integrator [11]) to predict the next time stamps’ system states, the
GNSTODE is endowed with the ability to adapt to the varying spatial and
temporal dependencies by making predictions through performing the integral
operation with regard to the learned spatial and temporal derivatives.

On the spatial domain, a common GNN based strategy [7] to calculate the
temporal dynamics at each time stamp is to use the iterative graph convolution
operation to increasingly account for the effects of distant neighboring
particles. To better adapt to varying spatial dependencies, we infer the
temporal dynamics at each time step by numerically solving a neural Ordinary
Differential Equation (ODE) on the spatial domain, where we use a GNN to model
the derivative of the temporal dynamics with regard to the neighborhood
radius. As the solution to the neural ODE, the final temporal dynamics are
obtained by integrating the GNN based derivative along the neighborhood
radius. Through learning from real-world data, the modelled derivative can
automatically explore the varying importance of neighboring particles with
varying distances to the central particles, by adaptively assigning large
derivative values to small neighborhood radiuses while making the large
neighborhood radiuses have small derivative values.

On the temporal domain, to make the model adaptive to the varying temporal
dependencies within a fixed time step, we use a neural network to parameterize
a continuous temporal dynamics function, that interpolates the known temporal
dynamics modeled by spatial neural ODE at finite time points. Based on the
fact that the modeled continuous temporal dynamics function is the derivative
of the system state with regard to time, we predict the next time’s system
state by solving another neural ODE on the temporal domain, by integrating the
continuous temporal dynamics from the last time stamp to the current time
stamp. Through using continuous temporal dynamics to model any possible state
change trends within a fixed time step, the proposed GNSTODE can make an
accurate prediction for the next time’s state. The accurate time point-wise
prediction makes no simulation errors accumulated along the simulation
rollout, which guarantees the stability and accuracy of the overall
simulation.

In the GNSTODE model, to make the spatial and temporal ODEs well collaborate
with each other, we integrate them into a united end-to-end learning
framework. We also develop an algorithm for training the proposed GNSTODE
model efficiently, in which numerical ODE solvers are leveraged to achieve the
fast feedforward prediction, while the adjoint sensitivity method [12] is
employed for instant error backpropagation. We conduct extensive experiments
on simulating the particle systems with varying spatial and temporal
dependencies. The results show that the proposed GNSTODE model achieves
significantly better simulation performance than state-of-the-art learning
based simulation methods in all cases. In summary, the contribution of this
paper is threefold:

  * •

To the best of our knowledge, we are the first to analyze the two key
challenges for learning to simulate particle systems: the varying spatial and
temporal dependencies. This not only well motivates the proposed approach but
also provides a high-level instruction for the follow-up development of the
learning based simulation models.

  * •

We propose a novel learning based simulation model GNSTODE that is able to
adapt to varying spatial and temporal dependencies in particle systems, by
creatively leveraging two coupled ODEs to model the varying dependencies in
spatial and temporal domains. An efficient training algorithm is also designed
to make sure GNSTODE is practical to real-world simulations.

  * •

To verify the effectiveness of the proposed GNSTODE model, we conduct
extensive simulation experiments on real-world particle systems with different
spatial and temporal dependencies. The experimental results show that the
proposed GNSTODE model consistently outperforms state-of-the-art learning
based simulation models in all cases by large margins.

The rest of this paper is organized as follows. In Section II, we first review
the related work. Then, the problem of learning to simulate is formally
defined in in Section III. Section IV then reviews the required preliminaries.
After that, Section V details the proposed GNSTODE model. Extensive
experiments are then presented in Section VI. Finally, we conclude this paper
in Section VII.

##  II Related Work

In this section, we review three streams of related work: physics system
simulation, Graph Neural Networks (GNNs) and Neural Ordinary Equations (Neural
ODEs).

###  II-A Physics System Simulation

Traditional physics system simulation is achieved through modeling the system
with basic physics laws. It mainly uses two rationales: 1) the Euler rationale
that treats the physics systems as continuums, with representative models of
Computational Fluid Dynamics (CFD) [13] and Finite Element Method (FEM) [14];
2) the Lagrangian rationale that models the physics systems as discrete
matters (particles), with examples of Discrete Element Method (DEM) [15],
Lattice Boltzmann Method (LBM) [16] and Smooth Particle Hydrodynamics (SPH)
[17], etc. As many real-world systems are described by particles, in this
paper, we consider the simulation under the Lagrangian rationale, i.e.,
particle system simulation.

As the underlying physics laws are hard to fully be uncovered, traditional
particle system simulation methods fail to yield satisfactory simulations for
many real-world systems. As a more advanced simulation paradigm, learning
based particle system simulation has been proposed to perceive the underlying
physics laws from real-world data with machine learning. Learning based
particle system simulation mainly leverages neural networks to model the
interactions between particles, with the purpose of capturing the relational
inductive biases [18]. Graph Interaction Networks (GINs) [19] is the pioneer
of learning to simulate particle systems. Neural Relational Inference (NRI)
[20] is then proposed to consider different interacting types. HOGN [11]
imports the Hamiltonian mechanics [21] as physics informed inductive biases
into INs for more accurate particle system simulation. To better characterize
the complex particle systems, DPI-Nets [8] try to capture the hierarchical and
long-range interactions between particles with the proposed Propagation
Networks. GNS [5] models the hierarchical and long-range particle-particle
interactions by the iterative massage-passing mechanism [22]. GNS variants are
then developed by improving its scalability [10] and adapting it to the mesh-
based simulation [4]. To make neural models reason like humans, VGPL [23]
reinforces the learning based particle system simulation with visual and
dynamics priors.

However, the existing learning based particle system simulation does not
account for the varying spatial and temporal dependencies in real-world
particle systems, limiting their performance. In this paper, we propose to use
Graph Neural ODEs to characterize the varying spatial and temporal
dependencies, proving a more effective method for real-world particle system
simulation.

###  II-B Graph Neural Networks (GNNs)

GNNs are one of the most powerful neural network architectures for learning
from graph structured data [7, 24], which are first proposed by [25]. After
that, a number of spectral GNNs [26] are then proposed, by performing feature
filtering on graph spectral domains. Many spatial variants have also been
developed for GNNs, like Graph Convolution Network (GCN) [27] that uses
neighborhood aggregation to update node representations, and Graph Attention
Network (GAT) [28] that considers the attention weights between neighboring
nodes when performing neighborhood aggregation. Among others, GNS [5] provides
a general GNNs based framework for particle system simulation, by modeling
particles as nodes and characterizing the interactions between particles by
edges. Recent research on GNNs includes improving GNNs’ scalability for
handling large graphs [29, 30], developing more powerful neighborhood
aggregation operators [31, 32] and learning informative graph representations
with self-supervised learning [33, 34], etc. GNNs have been applied to solve
many real-world problems in various disciplines, including time-series
prediction [35, 36], anomaly detection [37], hyperspectral image
classification [38], knowledge graph completion [39], molecular and material
property prediction [22, 40] and protein structure prediction [41], etc.

In this paper, we propose a novel framework to make the best of GNNs’
expressive power for accurate particle system simulation, by making the GNNs
based particle-particle interaction modeling well interplay with Neural ODEs
to capture the complex temporal and spatial dependencies in real-world
particle systems.

###  II-C Neural ODEs

Neural ODEs [42] are a new class of deep neural networks, which model deep
neural network transformations as the numerical solutions of ODEs. They have
been applied to construct continuous-depth residual networks [42], with much
lower memory cost than deep residual networks [43]. Neural ODEs also provide
an elegant way to model the continuous-time latent states for time series
[42]. As an example, ODE-RNNs [44] augment Recurrent Neural Networks (RNNs)
with Neural ODEs to accurate model hidden states for irregularly-sampled time
series. In addition to time series, Neural ODEs have also been applied to
model other temporal systems. LG-ODE [45] uses Neural ODEs to model the multi-
agent dynamic systems with known graph structures and irregularly-sampled
partial observations. TrajODE [46] applies Neural ODEs to model the
continuous-time trajectory dynamics in the noisy spatial trajectory data.
Besides, Neural ODEs have also been leveraged to model spatio-temporal point
processes [47], learn event functions [48], and dynamic graph[49, 50].

The success of Neural ODEs on modeling temporal dynamic systems inspires us to
leverage them to capture the complex continuous-time temporal dynamics in the
real-world particle systems. To fully unleash the expressive power of Neural
ODEs, we also apply Neural ODEs to characterize the dynamics in the spatial
domain, i.e., the varying spatial dependencies between particles with varying
distances.

##  III Problem Formulation

For the problem of learning to simulate, we are given a collection of observed
particle state trajectories $\\{\mathcal{X}\in\mathbb{R}^{n\times d\times
T}\\}$ as training samples, with $\mathcal{X}=\\{X_{1},X_{2},\cdots,X_{T}\\}$,
where $n$ is the number of particles in the particle system, $d$ is the number
of particle features (e.g., coordinates, velocities, accelerations, masses and
electric charges) for describing particles’ states, $T$ is the number of time
stamps of the trajectories, and $X_{t}\in\mathbb{R}^{n\times d}$ is the
particle system state at time stamp $t$, with its $i$-th column
$\bm{x}_{t,i}\in\mathbb{R}^{d}$ representing the feature vector of the $i$-th
particle. Using training samples, we aim to learn a simulator,
$\mathcal{S}:\mathbb{R}^{n\times d}\rightarrow\mathbb{R}^{n\times d}$, which
is able to predict the particle system state at time stamp $t$ with its
previous state at time stamp $t-1$, i.e., $\hat{X}_{t+1}=\mathcal{S}(X_{t})$.
The training objective is to minimize the difference between the predicted
state $\hat{X}_{t+1}$ and the ground-truth state $X_{t+1}$.

To construct the simulator $\mathcal{S}$, we first need to learn a temporal
dynamics function $\mathcal{D}:\mathbb{R}^{n\times
d}\rightarrow\mathbb{R}^{n\times d}$, for predicting the temporal dynamics at
each time stamp: $D_{t}=\mathcal{D}(X_{t})$. The temporal dynamics
$D_{t}\in\mathbb{R}^{n\times d}$ describes the trend by which the particle
system state $X_{t}$ evolves to the next time stamp’s state $X_{t+1}$. The
simulator $\mathcal{S}$ is then obtained by updating the particle system state
$X_{t}$ with its temporal dynamics $D_{t}$:
$\mathcal{S}(X_{t})=\mathcal{F}(X_{t},D_{t})=\hat{X}_{t+1}$, where
$\mathcal{F}:\mathbb{R}^{n\times d}\times\mathbb{R}^{n\times
d}\rightarrow\mathbb{R}^{n\times d}$ is the updating function.

At each time stamp $t$, we design a spatial Neural ODE to formulate the
temporal dynamics function $\mathcal{D}$, by considering the varying spatial
dependencies between particles. For characterizing the spatial dependencies
between particles, a spatial graph is constructed according to the particle
coordinates, $G_{t}=(V_{t},E_{t},X_{t})$, where $V_{t}$ is the set of nodes
with each node corresponding to a particle, $E_{t}\subseteq V_{t}\times V_{t}$
is the set of edges with each edge connecting two particles having a distance
smaller than a predefined threshold, and $X_{t}\in\mathbb{R}^{n\times d}$
describes the particle features. In the spatial Neural ODE based temporal
dynamics function formulation, the Graph Interaction Network (GIN) [19] will
be used for effectively capturing the interactions between particles.

To capture the varying temporal dependencies between consecutive time stamps,
we design a temporal Neural ODE to formulate the updating function
$\mathcal{F}$. Instead of using a simple one-step updating operation (e.g.,
addition), the temporal Neural ODE extrapolates the continuous temporal
dynamics between consecutive time stamps, so that the accumulative updating
along a continuous time slot can be achieved to accurately predict the next
time stamp’s state.

With the well-trained simulator $\mathcal{S}$, given the initial state of a
particle system $X_{1}$, we can simulate its state evolving trajectory
$\hat{\mathcal{X}}=\\{X_{1},\hat{X}_{2},\cdots,\hat{X}_{T}\\}$, by
sequentially applying the simulator $\mathcal{S}$ at each time stamp, with
$\hat{X}_{2}=\mathcal{S}(X_{1})$ and $\hat{X}_{t+1}=\mathcal{S}(\hat{X}_{t})$
for $t\geq 2$. As the simulator $\mathcal{S}$ is designed for effectively
capturing the particle systems’ complex evolving patterns, i.e., the varying
spatial and temporal dependencies, we can expect that the simulated particle
state trajectory $\hat{\mathcal{X}}$ will be consistent with the ground-truth
state trajectory $\mathcal{X}$.

##  IV Preliminaries

Figure 2: The overall framework of GNSTODE. The framework is composed of two
key components: the neural spatial ODE component is used to predict current
temporal dynamics, while the neural temporal ODE component is employed to
predict the next time stamp’s system state.

In this section, we review the preliminaries on Graph Interaction Network
(GIN) [18] and Neural ODEs [42].

###  IV-A Graph Interaction Network (GIN)

At each time stamp $t$, by formulating the particle-particle interaction
relations with a particle interaction network $G_{t}=(V_{t},E_{t},X_{t})$, we
can use the Graph Interaction Network (GIN) [18] to learn informative particle
representations that well encode the particle-particle interaction patterns.
The particle temporal dynamics can then be predicted from particle
representations by performing a non-linear transformation. GIN is a stacked
multi-layer neural network, with each layer performing a message-passing
operation. By denoting the $l$-th-layer input representation of the $i$-th
particle as $\bm{h}_{l,i}\in\mathbb{R}^{d_{l}}$, in the $l$-th layer, the
representation is updated as

$\displaystyle\bm{m}_{l,i}=\sum_{j\in\mathcal{N}(i)}\mathcal{M}_{l}(\bm{h}_{l,i},\bm{h}_{l,j},\bm{e}_{i,j}),$
(1) $\displaystyle\bm{h}_{l+1,i}=\mathcal{U}_{l}(\bm{h}_{l,i},\bm{m}_{l,i}),$

where $\mathcal{N}(i)$ is the set of particle $i$’s neighboring particles;
$\bm{e}_{i,j}\in\mathbb{R}^{d_{e}}$ represents the features of the edge
connecting particles $i$ and $j$, which can be the distance or interacting
forces between them; $\mathcal{M}_{l}(\cdot,\cdot,\cdot)$ is the $l$-th-layer
function for calculating the messages between interacting particles;
$\bm{m}_{l,i}$ represents the aggregated messages passed to particle $i$ in
the $l$-th layer; and $\mathcal{U}_{l}(\cdot,\cdot)$ is the function for
updating particle representations with aggregated messages in the $l$-th
layer.

With the message-passing mechanism serving as an effective relational
inductive bias, GIN can accurately predict particle temporal dynamics from the
particle-particle interaction patterns. For mathematical convenience, we
denote input particle representations of the $l$-th layer as
$H_{l}\in\mathbb{R}^{n\times d_{l}}$ with its $i$-th column being the $i$-th
particle’s representation $\bm{h}_{l,i}$, and set $H_{0}=X_{t}$, then we
reformulate the message passing operation in Eq. (1) as

$H_{l+1}=\mathrm{GIN}_{l}(H_{l}).$ (2)

###  IV-B Neural ODEs

Many deep neural networks (e.g., residual networks and recurrent neural
networks) work through iteratively operating a transformation on the networks’
hidden states:

$\bm{z}_{t+1}=\bm{z}_{t}+f_{t}(\bm{z}_{t}),$ (3)

where $t\in\\{0,1,\cdots,T\\}$ is the layer number,
$\bm{z}_{t}\in\mathbb{R}^{d}$ is the hidden state at the $t$-th layer and
$f_{t}(\cdot)$ is the $t$-th layer’s transformation function.

Neural ODEs [42] reformulate the discrete updating as a continuous integral
operation, by formulating the relation between hidden state $\bm{z}_{t}$ and
layer number $t$ with an ODE:

$\frac{d\bm{z}_{t}}{dt}=f(\bm{z}_{t},t),$ (4)

where $f(\cdot,\cdot)$ is the layer-dependent transformation function
parameterized by a neural network. Given the input state $\bm{z}_{0}$, the
hidden state at the $t$-th layer $\bm{z}_{t}$ is calculated by

$\bm{z}_{t}=\bm{z}_{0}+\int_{0}^{t}f(\bm{z}_{\tau},\tau)d\tau,$ (5)

which can be solved by numerical ODE solvers, such as Euler and Runge-Kutta
[42]. Similar to traditional deep neural networks, the learnable parameters in
$f(\cdot,\cdot)$ are updated by error backpropagation, where their gradients
can be fast computed by the adjoint sensitivity method [12].

##  V Methodology

In this section, we detail the proposed GNSTODE framework to learn the
simulator function $\mathcal{S}$ for predicting the next time stamp’s system
state according to the current state. To yield high-quality simulation, the
proposed GNSTODE framework effectively models the varying spatial dependencies
via a neural spatial ODE component and captures the varying temporal
dependencies with a neural temporal ODE component.

###  V-A Overall Framework

Fig. 2 shows the overall GNSTODE learning framework. At each time stampt $t$,
a neural spatial ODE is used to construct the temporal dynamics function
$\mathcal{D}:\mathbb{R}^{n\times d}\rightarrow\mathbb{R}^{n\times d}$, for
inferring the current temporal dynamics $D_{t}$ from the current system state
$X_{t}$. The neural spatial ODE is solved through an integral along the
varying distances between neighboring particles and central particles, to
characterize the varying interacting effects between particles with varying
distances. On the other hand, a neural temporal ODE is used to construct the
updating function $\mathcal{F}:\mathbb{R}^{n\times d}\times\mathbb{R}^{n\times
d}\rightarrow\mathbb{R}^{n\times d}$, for predicting the current system state
$\hat{X}_{t}$ from the last time stamp’s state $X_{t-1}$ and the last time
stamp’s temporal dynamics $D_{t-1}$. Through extrapolating the continuous
temporal dynamics function between the time stamps $t-1$ and $t$, the temporal
neural ODE can effectively capture the complex accumulative updating effects
between consecutive time stamps. The reconstruction loss at the time stamp
$\mathcal{L}_{t}$ is then calculated through comparing the predicted system
state $\hat{X}_{t}$ with the ground truth state $X_{t}$. The model is trained
by minimizing the overall construction loss $\mathcal{L}$, which is obtained
by summing up the reconstruction losses at all time stamps.

###  V-B Neural Spatial ODE Component

At each time stamp $t$, by representing the particle system state $X_{t}$ as a
spatial graph $G_{t}=(V_{t},E_{t},X_{t})$, we use the Graph Interaction
Network (GIN) [18] to model the interactions between particles. We aim to make
the best of GIN’s expressive power to predict the temporal dynamics $D_{t}$
accurately. GIN embeds the complex interactions between neighboring particles
into particle latent representations with the iterative message-passing
operation:

$H_{l+1}=\mathrm{GIN}_{l}(H_{l}),$ (6)

where $\mathrm{GIN}_{l}(\cdot)$ is the $l$-th GIN layer’s message passing
operation defined by Eq. (1) and $H_{l}$ is the input particle representations
of the $l$-th GIN layer with $H_{0}=X_{t}$.

To avoid over-smoothing and make the hidden representations more indicative to
particle original features, by augmenting a skip-connection operation, we
reformulate the GIN’s layer-wise updating in Eq. (6) as

$H_{l+1}=H_{l}+\mathrm{GIN}_{l}(H_{l}).$ (7)

By regarding the operation in Eq. (7) as the discrete approximation of a
continuous-layer updating, we can write the derivative of the hidden
representations $H_{l}$ with regard to the layer number $l$, $H^{\prime}_{l}$,
as,

$H^{\prime}_{l}=\frac{dH_{l}}{dl}=\mathrm{GIN}_{l}(H_{l}),$ (8)

where the GIN layer-wise operation $\mathrm{GIN}_{l}(\cdot)$ is used to
formulate the hidden representation derivative. Assuming the upper bound of
layer number is set to $L$, the final-layer particle representations can be
obtained by solving the ODE problem in Eq. (8) as:

$H_{L}=H_{0}+\int_{0}^{L}H^{\prime}_{l}dl=X_{t}+\int_{0}^{L}H^{\prime}_{l}dl.$
(9)

In implementation, we set $L=1$, which warrants $L$ to be automatically re-
scaled to arbitrary values with the learnable parameters in $H^{\prime}_{l}$.

The temporal dynamics at time stamp $t$, $D_{t}$, is finally obtained by
performing a multi-layer neural network based transformation
$\mathrm{NN}(\cdot)$ on $H_{L}$:

$D_{t}=\mathrm{NN}(H_{L}).$ (10)

###  V-C Neural Temporal ODE Component

At each time stamp $t$, given the system state $X_{t}$ and the inferred
temporal dynamics $D_{t}$, we can predict the next time stamp’s state with the
updating function $\mathcal{F}$:

$\hat{X}_{t+1}=\mathcal{F}(X_{t},D_{t}).$ (11)

A common implementation of $\mathcal{F}$ is to use the simple addition
operator:

$\hat{X}_{t+1}=X_{t}+D_{t}.$ (12)

However, such a naive updating operation cannot account for the complex
accumulative effects between time stamps $t$ and $t+1$. To achieve continuous-
time updating, we model the relation between the particle system state and
time as an ODE:

$\left\\{\begin{aligned}
&\frac{d\hat{X}_{\tau}}{d\tau}=D_{\tau},\;\;\tau\in[t,t+1],\\\
&D_{\tau}=D_{t}\;\;\mathrm{for}\;\tau=t.\end{aligned}\right.$ (13)

Here, we use a neural network to parameterize the temporal dynamics function
$D_{\tau}$ with input time variable $\tau\in[t,t+1]$. The next time stamp’s
state is then predicted by solving Eq. (13) as

$\hat{X}_{t+1}=X_{t}+\int_{t}^{t+1}D_{\tau}d\tau.$ (14)

###  V-D Loss Function

For each time stamp $t\geq 2$, by comparing the predicted particle system
state $\hat{X}_{t}$ with the ground-truth state $X_{t}$, we calculate the
reconstruction loss as

$\mathcal{L}_{t}=\|\hat{X}_{t}-X_{t}\|_{\mathrm{F}}^{2},$ (15)

where $\|\cdot\|_{\mathrm{F}}$ is the Frobenius norm. The overall
reconstruction loss of the whole state trajectory is obtained by summing up
the reconstruction losses at all time stamps:

$\mathcal{L}=\sum_{t=2}^{T}\mathcal{L}_{t}.$ (16)

The parameters of the proposed GNSTODE is optimized by minimizing the
reconstruction loss $\mathcal{L}$ of each training state trajectory. We use
error propagation to update model parameters, where the Runge-Kutta method
[42] is used to numerically calculate the ODE solutions in Eq. (9) and (14) in
the feedforward process, and the gradients of the ODE parameters are fast
calculated by the adjoint sensitivity method [12] in the backpropagation
process.

Algorithm Description and Time Complexity. Algorithm 1 describes the overall
workflow for training the proposed GNSTODE model. We first initialize model
parameters with random numbers. The model parameters are then updated by
stochastic gradient descent: randomly sample a batch of training samples at
each iteration, and update model parameters with gradient descent by reducing
the corresponding reconstruction loss. Taking the maximum number of epochs
max_epoch as a constant, the time complexity of Algorithm 1 is
$O((nd^{2}+nds_{1}+nds_{2})T|\mathcal{X}|)$, where $s_{1}$ and $s_{2}$ are
respectively the numbers of steps adopted by the Runge-Kutta method for
calculating the ODE solutions in Eq. (9) and (14), and $|\mathcal{X}|$ is the
number of training samples in the training set $\mathcal{X}$.

Algorithm 1 The algorithm for training GNSTODE

Input: The particle state trajectory training set
$\\{\mathcal{X}\in\mathbb{R}^{n\times d\times T}\\}$.
Output: The trained GNSTGODE model.

1: Initialize model parameters with random numbers.

2: for epoch in 1, 2, $\cdots$, max_epoch do

3: $\mathcal{B}\leftarrow$ Randomly split $\mathcal{X}$ into batches of size
$B$;

4: for each batch in $\mathcal{B}$ do

5: Calculate temporal dynamics at time stamps $t=2,3,\cdots,T$ with Eq. (10);

6: Predict particle system states at time stamps $t=2,3,\cdots,T$ with Eq.
(14);

7: Update model parameters by reducing the reconstruction loss in Eq. (16)
with gradient descent;

8: end for

9: end for

10: return the trained GNSTGODE model.

##  VI Experiments

In this section, we conduct experiments to verify the effectiveness of the
proposed GNSTODE model for particle system simulation. We first introduce the
datasets used for experiments and the experimental setup. Then, we demonstrate
the experimental results, including the performance comparison with baselines,
ablation study and parameter sensitivity study.

TABLE I: Dataset statistics.

Dataset | #Particles | #Particle Features
---|---|---
Gravity | 20 | 5
100 | 5
500 | 5
Coulomb | 20 | 6
100 | 6
500 | 6

###  VI-A Benchmark Datasets

We evaluate the proposed GNSTODE model on six widely used particle simulation
benchmark datasets [10], including two types of particle systems (Gravity and
Coulomb) and three different particle scales (20, 100, 500) for each type.

In the Gravity particle system, each particle $i$ is descried by a five-
dimensional feature vector $[m_{i},x_{i},y_{i},\dot{x}_{i},\dot{y}_{i}]$,
where $m_{i}$ is the mass of the particle and $\bm{x}_{i}=[x_{i},y_{i}]$ is
the coordinate of the particle in the two-dimensional Euclidean space and
$\dot{\bm{x}}_{i}=[\dot{x}_{i},\dot{y}_{i}]$ is the velocity of the particle.
Differently, each particle $i$ is described by a 6-dimensional feature vector
$[m_{i},c_{i},x_{i},y_{i},\dot{x}_{i},\dot{y}_{i}]$ in the Coulomb system,
with an additional feature dimension – the electric charge $c_{i}$. The
statistics of the particle systems is summarized in Table I.

For each particle system with varying types and varying scales, we
respectively simulate 100 training, 20 validation and 20 testing trajectories
with 200 time stamps for each, using the Leapfrog simulator [10]. The Leapfrog
simulator leverages Newton’s second law to update particle coordinates and
velocities, and respectively utilize the Gravity and Coulomb forces to
formulate particle accelerations in the Gravity and Coulomb systems. In the
Gravity system, the particle acceleration is calculated as

$\bm{a}_{i}=-G\sum_{j\neq i}\frac{m_{j}}{\|\bm{x}_{i}-\bm{x}_{j}\|_{2}^{2}},$
(17)

where $G$ is the gravity constant. In the Coulomb system, the acceleration is
instead computed as

$\bm{a}_{i}=k\frac{1}{m_{i}}\sum_{j\neq i}\frac{c_{i}\cdot
c_{j}}{\|\bm{x}_{i}-\bm{x}_{j}\|_{2}^{2}},$ (18)

where $k$ is the Coulomb constant.

For each simulation, particles’ initial coordinates are determined by the
random sampling from a two-dimensional uniform random distribution with
intensity 0.42 (the average particle number at per $1\times 1$ square);
particle masses are uniformly set to 1; particle velocities are initialised by
sampling over the uniform distribution on $(-1,1)$; particle charges are
initialised by the random sampling from the uniform distribution on
$(0.5,1.5)$ for the Coulomb system; the Gravity and Coulomb constants are
respectively set as $G=2$ and $k=2$; and the time step is set as $0.01$.

###  VI-B Experimental Settings

To provide a comprehensive evaluation on the simulation performance of the
proposed GNSTODE model, we configure three different experimental settings:

  * •

Varying Scales. This setting evaluates the robustness of the GNSTODE model
against varying particle scales. As is described, the three different particle
scales (20, 100, 500) are respectively used for the Gravity and Coulomb
systems for the evaluation purpose.

  * •

Varying Intensities. This setting aims to mimic the varying spatial
dependencies between particles. Fixing the particle system scale to 20, we
respectively simulate three different system variants with initial particle
intensities 0.083 0.42 and 1.63 for the Gravity and Coulomb systems.

  * •

Varying Time Steps. This setting mimics the varying temporal dependencies. On
the Gravity and Coulomb particle trajectories with scale 20 and initial
particle intensity 0.42, we respectively perform three downsampling operations
by keeping only every 5th, 10th and 20th time stamps and benchmark on the
downsampled trajectories.

TABLE II: Simulation performance comparison on the Gravity and Coulomb systems
under the setting of varying scales. The best and second best performers are
respectively highlighted by the boldface and underline.

Dataset |  Metrics $\downarrow$ | Scale | GNS | HOGN | HDGN | HHOGN | GNSTODE
---|---|---|---|---|---|---|---
Gravity | RMSE | 20 | 0.8654 | 1.1053 | 0.8114 | 1.1211 | 0.7053
100 | 1.7098 | 1.7120 | 1.0471 | 1.5190 | 0.9937
500 | 1.6813 | 1.7943 | 1.0033 | 1.5533 | 0.7447
Energy Error | 20 | 0.0659 | 0.4227 | 0.0584 | 0.4237 | 0.0167
100 | 0.4116 | 0.2160 | 0.0466 | 0.1918 | 0.0199
500 | 0.1268 | 0.2201 | 0.0134 | 0.0858 | 0.0029
Coulomb | RMSE | 20 | 0.8471 | 1.1022 | 0.7770 | 0.7475 | 0.6797
100 | 0.6519 | 0.7459 | 0.4968 | 0.7857 | 0.4505
500 | 0.6921 | 0.7857 | 0.5562 | 0.5442 | 0.3334
Energy Error | 20 | 0.0642 | 0.4071 | 0.0248 | 0.4301 | 0.0116
100 | 0.7375 | 0.3810 | 0.8405 | 0.4075 | 0.3273
500 | 0.7373 | 0.8901 | 0.2175 | 0.8614 | 0.0423

TABLE III: Simulation performance comparison on the Gravity and Coulomb
systems under the setting of varying intensities. The best and second best
performers are respectively highlighted by the boldface and underline.

Dataset |  Metrics $\downarrow$ | Intensity | GNS | HOGN | HDGN | HHOGN | GNSTODE
---|---|---|---|---|---|---|---
Gravity | RMSE | 0.083 | 0.1971 | 0.5036 | 0.3124 | 0.5055 | 0.1932
0.42 | 0.8654 | 1.1053 | 0.8114 | 1.1210 | 0.6936
1.63 | 1.5903 | 1.5109 | 1.2862 | 1.5668 | 1.2595
Energy Error | 0.083 | 0.0094 | 0.1804 | 0.0220 | 0.1762 | 0.0090
0.42 | 0.0658 | 0.4226 | 0.0584 | 0.4237 | 0.0149
1.63 | 0.5787 | 0.7136 | 0.0484 | 0.8233 | 0.0464
Coulomb | RMSE | 0.083 | 0.2139 | 0.2049 | 0.2083 | 0.2556 | 0.2044
0.42 | 0.9040 | 1.4009 | 0.9442 | 1.3949 | 0.7468
1.63 | 0.4660 | 0.4574 | 0.4568 | 0.4851 | 0.2009
Energy Error | 0.083 | 1.1377 | 0.7508 | 0.9783 | 1.5492 | 0.9576
0.42 | 0.0752 | 0.0933 | 0.0499 | 0.8810 | 0.0313
1.63 | 3.7563 | 3.7467 | 3.6365 | 4.7369 | 1.1669

###  VI-C Baseline Methods

We compare the proposed GNSTODE model with the following four competitive
learning based particle system simulation baselines:

  * •

GNS [5] employs multi-step message passing to capture the complex particle-
particle interactions.

  * •

HOGN [11] utilizes the GIN [18] to model the Hamiltonian of particle systems
and applies the Hamiltonian to capture the complex interactions between
particles.

  * •

HDGN [10] constructs a hierarchical graph to model particles’ multi-level
neighbors and leverages the GIN [18] to capture the cross-level interactions.

  * •

HHOGN [10] extends the HOGN model [11] to the hierarchical particle
interaction graph for achieving more accurate particle system simulation.

TABLE IV: Simulation performance comparison on the Gravity and Coulomb systems
under the setting of varying time steps. The best and second best performers
are respectively highlighted by the boldface and underline.

Dataset |  Metrics $\downarrow$ | Step | GNS | HOGN | HDGN | HHOGN | GNSTODE
---|---|---|---|---|---|---|---
Gravity | RMSE | 5 | 0.9557 | 1.4083 | 0.8817 | 1.4068 | 0.7517
10 | 1.0739 | 1.5242 | 1.0213 | 1.5292 | 0.7965
20 | 1.3126 | 1.6155 | 1.5292 | 1.7314 | 0.8833
Energy Error | 5 | 0.0971 | 0.9396 | 0.9787 | 0.9835 | 0.0197
10 | 0.1571 | 1.5965 | 0.1444 | 1.4450 | 0.0414
20 | 0.8050 | 1.7440 | 0.2874 | 2.5928 | 0.0512
Coulomb | RMSE | 5 | 0.9040 | 1.4006 | 0.9442 | 1.3949 | 0.7468
10 | 0.9903 | 1.5084 | 1.0656 | 1.5188 | 0.7436
20 | 1.1939 | 1.5937 | 1.2679 | 1.6505 | 0.8487
Energy Error | 5 | 0.0751 | 0.9325 | 0.0849 | 0.8810 | 0.0313
10 | 0.0663 | 1.4966 | 0.1284 | 1.4442 | 0.0225
20 | 0.2994 | 1.7826 | 0.3376 | 2.2266 | 0.0489

###  VI-D Evaluation Metrics

Following [10], we evaluate the particle system simulation performance using
the following two metrics:

  * •

RMSE (Root Mean Square Error). We expect the simulated particle trajectories
to be as close to the ground truth as possible. RMSE measures the discrepancy
between the predicted and ground-truth trajectories:

$\text{RMSE}=\sqrt{\frac{1}{N}\sum_{i=1}^{N}\sum_{t=2}^{T}\|X^{c}_{i,t}-\hat{X}^{c}_{i,t}\|_{\mathrm{F}}^{2}},$
(19)

where $N$ is the number of test particle system trajectories, $T$ is the
number of time stamps of the test trajectories, and $X^{c}_{i,t}$ and
$\hat{X}^{c}_{i,t}$ are respectively the ground-truth and predicted particle
coordinates at the time stamp $t$.

  * •

Energy Error. Both the Gravity and Coulomb systems follow the principle of
energy conservation, meaning that the total energy of the whole system should
be constant over time. The Energy Error measures the extent to which the
simulation violates the principle of energy conservation:

$\text{Energy
Error}=\frac{1}{N}\sum_{i=1}^{N}\frac{H_{i,1}-\hat{H}_{i,T}}{H_{i,1}},$ (20)

where $H_{i,1}$ and $\hat{H}_{i,T}$ are respectively the Hamiltonian (total
energy) of the ground-truth system state at time stamp $1$ and the predicted
system state at time stamp $T$ of the $i$-th test trajectory.

For both the RMSE and Energy Error metrics, a lower score indicates a better
simulation performance.

###  VI-E Implementation Details

All experiments are conducted on a personal computer with the Ubuntu 18.04 OS,
an NVIDIA GeForce RTX 2080Ti (12GB memory) GPU, an Intel Core i9-10900X (3.70
GHz) CPU and a 32 GB RAM. Default parameter configurations are used for
implementing baseline models.

When implementing the GNSTODE model, 15-nearest neighbor graph is used to
construct the spatial graph and the latent particle embeddings’ dimension is
set same as the original particle features’ dimension (5 for the Gravity
system and 6 for the Coulomb system). The maximum number of epochs is set to
200. At each epoch, we collect the paired particle system states at every two
consecutive time stamp from the training trajectories, shuffle them and split
them into a number of batches, then iterate each batch of particle state pairs
to minimize the corresponding reconstruction errors. The batch size is set to
50. We optimize the GNSTODE model with the Adam optimizer [51].

For the proposed GNSTODE model and baseline models, the validation set is
employed to select the best epoch, and the test set is used to evaluate the
simulation performance. On each dataset and under each setting, we run the
proposed GNSTODE model and baseline models for five times with different
random parameter initializations and report the averaged RMSE and Energy Error
scores as models’ simulation performance.

###  VI-F Performance Comparison

Tables II-IV compare the simulation performance of different models under the
three different settings, i.e., varying scales, varying intensities and
varying time steps. For each comparison, the best and second best performers
are respectively highlighted by the boldface and underline. From Tables II-IV,
we can find that the proposed GNSTODE model consistently achieves the best
RMSE and Energy Error scores on each dataset and under each setting, except
for the second best Energy Error on the Couloub system with intensity 0.083.

Table II shows that the proposed GNSTODE model significantly outperforms the
compared baseline methods on different scales of particle systems. This
verifies that the proposed GNSTODE model is robust to the varying scales of
particle systems. As is shown in Table II, models tend to achieve better
simulation performance on the particle systems with larger scales, which
provide more supervision on particle movement patterns.

From Table III, we can see that the proposed GNSTODE model achieves superior
RMSE and Energy Error scores than baseline models on particle systems with
varying intensities. This is mainly contributed by the effectiveness of the
proposed neural spatial ODE component in capturing the varying spatial
dependencies between particles. The slight performance drop on the particle
systems with larger intensities indicates that the larger spatial dependency
variations bring more difficulties to the learning based particle system
simulation.

Table IV demonstrates the best simulation performance of the proposed GNSTODE
model on the particle systems with different time steps. This is consistent
with the fact that the designed neural temporal ODE makes the proposed GNSTODE
able to adapt the varying temporal dependencies between every two consecutive
time stamps. As large time steps incorporate more temporal dependency
variations, long-time-step simulation usually yields higher RMSE and Energy
Error scores.

###  VI-G Simulation Visualization

In Fig. 3 and 4, we visualize the simulations produced by the GNS, HDGN and
the proposed GNSTODE model on a Gravity system and a Coulomb system with scale
20 and intensity 0.42, and compare them with the ground-truth trajectories.
Through the comparison, we can find that the simulations produced by the GNS
and HDGN baseline models have significant discrepancies with the ground-truth
trajectories, while the GNSTODE simulations can be almost perfectly aligned to
the ground-truth trajectories. This further verifies the superior simulation
performance of the proposed GNSTODE model over the baseline models.

(a) Ground Truth

(b) GNS

(c) HDGN

(d) GNSTODE

Figure 3:  The ground-truth trajectory and simulations of a Gravity system
with scale 20, intensity 0.42 and time step 1.

(a) Ground Truth

(b) GNS

(c) HDGN

(d) GNSTODE

Figure 4: The ground-truth trajectory and simulations of a Coulomb system with
scale 20, intensity 0.42 and time step 1.

###  VI-H Ablation Study

To further verify the effectiveness of the proposed neural spatial and
temporal ODE components, we also conduct an ablation study by comparing the
full GNSTODE model with its two ablated variants:

  * •

w/o Spatial ODE. For this variant, we replace the neural spatial ODE component
in Eq. (9) with the GIN message-passing operation in Eq.(7) to model the
interactions between particles.

  * •

w/o Temporal ODE. For this variant, we replace the neural temporal ODE
component in Eq. (14) with the simple addition operation in Eq. (12) to
predict the next time stamp’s system state.

Tables V-VI respectively compare the RMSE and Energy Error scores of the full
GNSTODE model and its ablated variants on the Gravity system with different
scale, intensity and time step configurations. The best performer is
highlighted by the boldface for each comparison. From Tables V-VI, we can find
that the full GNSTODE consistently outperforms its ablated variants – w/o
Spatial ODE and w/o Temporal ODE – by large margins under all configurations.
This proves that the neural spatial and temporal ODE components are both
important to the high-quality GNSTODE simulations. Ablating either of them
will result in inferior simulation performance.

TABLE V:  The RMSE scores of the full GNSTODE model and its ablated variants on the Gravity system with different scale, intensity and time step configurations. The best performers are highlighted by the boldface.  | Scale | Intensity | Time Step
---|---|---|---
Method |  20 |  100 |  0.083 |  1.63 |  5 |  10
Full GNSTODE |  0.7053 |  0.9937 |  0.1934 |  1.2595 |  0.7351 |  0.7965
w/o Spatial ODE |  0.8414 |  1.5549 |  0.2142 |  1.3644 |  0.8825 |  0.9901
w/o Temporal ODE |  0.7437 |  1.0134 |  0.1954 |  1.2839 |  0.7422 |  0.8434
TABLE VI:  The Energy Error scores of the full GNSTODE model and its ablated variants on the Gravity system with different scale, intensity and time step configurations. The best performers are highlighted by the boldface.  | Scale | Intensity | Time Step
---|---|---|---
Method |  20 |  100 |  0.083 |  1.63 |  5 |  10
Full GNSTODE |  0.0167 |  0.0199 |  0.0090 |  0.0464 |  0.0196 |  0.0414
w/o Spatial ODE |  0.0405 |  0.0308 |  0.0096 |  0.0495 |  0.0210 |  0.1304
w/o Temporal ODE |  0.0446 |  0.1557 |  0.0117 |  0.0945 |  0.0378 |  0.0912

###  VI-I Parameter Sensitivity Study

Figure 5: Parameter sensitivity with regard to the batch size, epoch number
and embedding size.

In this section, we study the sensitivity of the proposed GNSTODE model with
regard to the three main model parameters, the batch size, the maximum number
of epochs (epoch number) and the latent embeddings’ dimension, on the Gravity
dataset with scale, intensity and time step respectively set to 20, 0.42 and
1. In turn, we fix any two of the three parameters as default values, and
study the performance change of the GNSTODE model when the value of the
remainder parameter varies in a proper range. Fig. 5 shows the model’s
sensitivity with regard to the three parameters.

From Fig. 5, we can find that the model’s performance degrades with the
increase of batch size. The large batch size makes the gradient descent hard
to capture the uniqueness of each particle system’s updating with varying
spatial and temporal dependencies, resulting in inferior simulation
performance. Fig. 5 shows that the GNSTODE’s performance gradually increases
and then remains at a stable level with the increase of the epoch number. The
study indicates that the epoch number is required to be large enough to make
sure the model is sufficiently trained, while an overlarge epoch number is not
necessary when the model training is saturated. Fig. 5 demonstrates that the
GNSTODE’s performance is not sensitive to the dimension of particle latent
embeddings. This is mainly contributed by the strong learning ability of the
proposed neural spatial and temporal ODE components for modeling the
complicated particle system state updates, making GNSTODE still able to
achieve satisfactory performance even with a small particle embedding
dimension.

##  VII Conclusion and Future Work

In this paper, we employed the GNNs and dual neural ODEs in the both spatial
and temporal domains to tackle the problem of learning to simulate particle
systems for the first time, established a new paradigm for physics learning
and outlined a comprehensive solution procedure. Moreover, to verify the
efficacy of our model, we designed three experimental settings on the Gravity
and Coulomb particle systems, i.e., varying scales, varying intensities and
varying time steps, which will also inspire the follow-up research to evaluate
the learning based particle system simulation models comprehensively. The
experimental results show that our model consistently achieves better
simulation performance than the state-of-the-art learning based simulation
models in all cases. More importantly, the results confirm our model’s
robustness to varying scales of particle systems and the ability to
automatically adapt to particle systems with varying spatial and temporal
dependencies.

The changes in the actual physical world are significantly more complex than
we believe. Although our model provides a solution to modeling the complex
spatial and temporal dependencies in particle systems, it still has a long way
to go before we can fully comprehend the laws of the real world. In addition
to addressing the studied simulation problem, we can make further explorations
by studying a series of related tasks, such as modeling the environment
interaction (i.e., modeling the interactions between particles and boundaries)
and modeling the interactions between particles with different material
compositions, which are great value to process engineering. Another
interesting future work is to study how to effectively inject physical
inductive bias into the learning based simulation models to make them more
consistent with the conventional Euler and Lagrangian simulation methods [2]
in rationale. Overall, the research on AI for physics science needs to be
further explored thoroughly.

## References

  * [1] C. Herfeld and M. Ivanova, “Introduction: first principles in science—their status and justification,” _Synthese_ , vol. 198, no. 14, pp. 3297–3308, 2021.
  * [2] S. Subramaniam, “Lagrangian–eulerian methods for multiphase flows,” _Progress in Energy and Combustion Science_ , vol. 39, no. 2-3, pp. 215–245, 2013.
  * [3] S. L. Brunton, B. R. Noack, and P. Koumoutsakos, “Machine learning for fluid mechanics,” _Annual Review of Fluid Mechanics_ , vol. 52, pp. 477–508, 2020\.
  * [4] T. Pfaff, M. Fortunato, A. Sanchez-Gonzalez, and P. W. Battaglia, “Learning mesh-based simulation with graph networks,” _arXiv preprint arXiv:2010.03409_ , 2020.
  * [5] A. Sanchez-Gonzalez, J. Godwin, T. Pfaff, R. Ying, J. Leskovec, and P. Battaglia, “Learning to simulate complex physics with graph networks,” in _International Conference on Machine Learning_. PMLR, 2020, pp. 8459–8468.
  * [6] J. Shlomi, P. Battaglia, and J.-R. Vlimant, “Graph neural networks in particle physics,” _Machine Learning: Science and Technology_ , vol. 2, no. 2, p. 021001, jan 2021. [Online]. Available: https://doi.org/10.1088/2632-2153/abbf9a
  * [7] Z. Wu, S. Pan, F. Chen, G. Long, C. Zhang, and P. S. Yu, “A comprehensive survey on graph neural networks,” _IEEE Transactions on Neural Networks and Learning Systems_ , vol. 32, no. 1, pp. 4–24, 2021.
  * [8] Y. Li, J. Wu, R. Tedrake, J. B. Tenenbaum, and A. Torralba, “Learning particle dynamics for manipulating rigid bodies, deformable objects, and fluids,” in _International Conference on Learning Representations_ , 2019.
  * [9] A. Sanchez-Gonzalez, N. Heess, J. T. Springenberg, J. Merel, M. Riedmiller, R. Hadsell, and P. Battaglia, “Graph networks as learnable physics engines for inference and control,” _arXiv preprint arXiv:1806.01242_ , 2018.
  * [10] K. Martinkus, A. Lucchi, and N. Perraudin, “Scalable graph networks for particle simulations,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 35, no. 10, 2021, pp. 8912–8920.
  * [11] A. Sanchez-Gonzalez, V. Bapst, K. Cranmer, and P. Battaglia, “Hamiltonian graph networks with ode integrators,” _arXiv preprint arXiv:1909.12790_ , 2019.
  * [12] L. S. Pontryagin, _Mathematical theory of optimal processes_. CRC press, 1987.
  * [13] J. H. Ferziger, M. Perić, and R. L. Street, _Computational methods for fluid dynamics_. Springer, 2002, vol. 3.
  * [14] O. C. Zienkiewicz, R. L. Taylor, and J. Z. Zhu, _The finite element method: its basis and fundamentals_. Elsevier, 2005.
  * [15] H. Zhu, Z. Zhou, R. Yang, and A. Yu, “Discrete particle simulation of particulate systems: theoretical developments,” _Chemical Engineering Science_ , vol. 62, no. 13, pp. 3378–3396, 2007.
  * [16] S. Chen and G. D. Doolen, “Lattice boltzmann method for fluid flows,” _Annual review of fluid mechanics_ , vol. 30, no. 1, pp. 329–364, 1998.
  * [17] J. J. Monaghan, “Smoothed particle hydrodynamics,” _Annual review of astronomy and astrophysics_ , vol. 30, no. 1, pp. 543–574, 1992.
  * [18] P. W. Battaglia, J. B. Hamrick, V. Bapst, A. Sanchez-Gonzalez, V. Zambaldi, M. Malinowski, A. Tacchetti, D. Raposo, A. Santoro, R. Faulkner _et al._ , “Relational inductive biases, deep learning, and graph networks,” _arXiv preprint arXiv:1806.01261_ , 2018.
  * [19] P. Battaglia, R. Pascanu, M. Lai, D. J. Rezende _et al._ , “Interaction networks for learning about objects, relations and physics,” in _Advances in neural information processing systems_ , 2016, pp. 4502–4510.
  * [20] T. Kipf, E. Fetaya, K.-C. Wang, M. Welling, and R. Zemel, “Neural relational inference for interacting systems,” in _International Conference on Machine Learning_. PMLR, 2018, pp. 2688–2697.
  * [21] S. Greydanus, M. Dzamba, and J. Yosinski, “Hamiltonian neural networks,” in _Advances in neural information processing systems_ , vol. 32, 2019.
  * [22] J. Gilmer, S. S. Schoenholz, P. F. Riley, O. Vinyals, and G. E. Dahl, “Neural message passing for quantum chemistry,” _arXiv preprint arXiv:1704.01212_ , 2017.
  * [23] Y. Li, T. Lin, K. Yi, D. Bear, D. L. Yamins, J. Wu, J. B. Tenenbaum, and A. Torralba, “Visual grounding of learned physical models,” in _ICML_ , 2020\.
  * [24] D. Zhang, J. Yin, X. Zhu, and C. Zhang, “Network representation learning: A survey,” _IEEE transactions on Big Data_ , vol. 6, no. 1, pp. 3–28, 2020\.
  * [25] M. Gori, G. Monfardini, and F. Scarselli, “A new model for learning in graph domains,” in _Proceedings. 2005 IEEE International Joint Conference on Neural Networks, 2005._ , vol. 2. IEEE, 2005, pp. 729–734.
  * [26] M. Defferrard, X. Bresson, and P. Vandergheynst, “Convolutional neural networks on graphs with fast localized spectral filtering,” in _Advances in neural information processing systems_ , 2016, pp. 3844–3852.
  * [27] T. N. Kipf and M. Welling, “Semi-supervised classification with graph convolutional networks,” in _International Conference on Learning Representations_ , 2017, pp. 1–14.
  * [28] P. Veličković, G. Cucurull, A. Casanova, A. Romero, P. Liò, and Y. Bengio, “Graph attention networks,” in _International Conference on Learning Representations_ , 2018, pp. 1–12.
  * [29] F. Wu, A. Souza, T. Zhang, C. Fifty, T. Yu, and K. Weinberger, “Simplifying graph convolutional networks,” in _International conference on machine learning_. PMLR, 2019, pp. 6861–6871.
  * [30] W. L. Hamilton, R. Ying, and J. Leskovec, “Inductive representation learning on large graphs,” in _Proceedings of the 31st International Conference on Neural Information Processing Systems_ , 2017, pp. 1025–1035.
  * [31] K. Xu, W. Hu, J. Leskovec, and S. Jegelka, “How powerful are graph neural networks?” _arXiv preprint arXiv:1810.00826_ , 2018.
  * [32] J. Klicpera, S. Weißenberger, and S. Günnemann, “Diffusion improves graph learning,” _Advances in Neural Information Processing Systems_ , vol. 32, pp. 13 354–13 366, 2019.
  * [33] S. Pan, R. Hu, S.-f. Fung, G. Long, J. Jiang, and C. Zhang, “Learning graph embedding with adversarial training methods,” _IEEE transactions on cybernetics_ , vol. 50, no. 6, pp. 2475–2487, 2019.
  * [34] Y. Liu, S. Pan, M. Jin, C. Zhou, F. Xia, and P. S. Yu, “Graph self-supervised learning: A survey,” _arXiv preprint arXiv:2103.00111_ , 2021.
  * [35] Z. Wu, S. Pan, G. Long, J. Jiang, and C. Zhang, “Graph wavenet for deep spatial-temporal graph modeling,” _arXiv preprint arXiv:1906.00121_ , 2019\.
  * [36] Z. Wu, S. Pan, G. Long, J. Jiang, X. Chang, and C. Zhang, “Connecting the dots: Multivariate time series forecasting with graph neural networks,” in _Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, 2020, pp. 753–763.
  * [37] Y. Liu, Z. Li, S. Pan, C. Gong, C. Zhou, and G. Karypis, “Anomaly detection on attributed networks via contrastive self-supervised learning,” _IEEE Transactions on Neural Networks and Learning Systems_ , 2021.
  * [38] S. Wan, C. Gong, P. Zhong, S. Pan, G. Li, and J. Yang, “Hyperspectral image classification with context-aware dynamic graph convolutional network,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 59, no. 1, pp. 597–612, 2020.
  * [39] S. Ji, S. Pan, E. Cambria, P. Marttinen, and P. S. Yu, “A survey on knowledge graphs: Representation, acquisition, and applications,” _IEEE Transactions on Neural Networks and Learning Systems_ , pp. 1–21, 2021.
  * [40] M. Karamad, R. Magar, Y. Shi, S. Siahrostami, I. D. Gates, and A. B. Farimani, “Orbital graph convolutional neural network for material property prediction,” _Physical Review Materials_ , vol. 4, no. 9, p. 093801, 2020\.
  * [41] J. Jumper, R. Evans, A. Pritzel, T. Green, M. Figurnov, O. Ronneberger, K. Tunyasuvunakool, R. Bates, A. Žídek, A. Potapenko _et al._ , “Highly accurate protein structure prediction with alphafold,” _Nature_ , vol. 596, no. 7873, pp. 583–589, 2021.
  * [42] R. T. Chen, Y. Rubanova, J. Bettencourt, and D. Duvenaud, “Neural ordinary differential equations,” _arXiv preprint arXiv:1806.07366_ , 2018.
  * [43] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 770–778.
  * [44] Y. Rubanova, R. T. Chen, and D. Duvenaud, “Latent odes for irregularly-sampled time series,” _arXiv preprint arXiv:1907.03907_ , 2019.
  * [45] Z. Huang, Y. Sun, and W. Wang, “Learning continuous system dynamics from irregularly-sampled partial observations,” _arXiv preprint arXiv:2011.03880_ , 2020.
  * [46] Y. Liang, K. Ouyang, H. Yan, Y. Wang, Z. Tong, and R. Zimmermann, “Modeling trajectories with neural ordinary differential equations.” in _International Joint Conference on Artificial Intelligence_ , 2021, pp. 1498–1504.
  * [47] R. T. Chen, B. Amos, and M. Nickel, “Neural spatio-temporal point processes,” _arXiv preprint arXiv:2011.04583_ , 2020.
  * [48] ——, “Learning neural event functions for ordinary differential equations,” _arXiv preprint arXiv:2011.03902_ , 2020.
  * [49] M. Jin, Y. Zheng, Y.-F. Li, S. Chen, B. Yang, and S. Pan, “Multivariate time series forecasting with dynamic graph neural odes,” _IEEE Transactions on Knowledge and Data Engineering_ , 2022.
  * [50] M. Jin, Y.-F. Li, and S. Pan, “Neural temporal walks: Motif-aware representation learning on continuous-time dynamic graphs,” in _Advances in Neural Information Processing Systems_ , 2022.
  * [51] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _arXiv preprint arXiv:1412.6980_ , 2014.