# Reduced-Order Neural Operators: Learning Lagrangian Dynamics on Highly
Sparse Graphs

Hrishikesh Viswanath
Purdue University
&Yue Chang
University
of Toronto
&Julius Berner
California Institute
of Technology
&Peter Yichen Chen
Massachusetts Institute
of Technology
&Aniket Bera
Purdue University

###### Abstract

We present a neural operator architecture to simulate Lagrangian dynamics,
such as fluid flow, granular flows, and elastoplasticity. Traditional
numerical methods, such as the finite element method (FEM), suffer from long
run times and large memory consumption. On the other hand, approaches based on
graph neural networks are faster but still suffer from long computation times
on dense graphs, which are often required for high-fidelity simulations. Our
model, GIOROM or Graph Interaction Operator for Reduced-Order Modeling, learns
temporal dynamics within a reduced-order setting, capturing spatial features
from a highly sparse graph representation of the input and generalizing to
arbitrary spatial locations during inference. The model is geometry-aware and
discretization-agnostic and can generalize to different initial conditions,
velocities, and geometries after training. We show that point clouds of the
order of 100,000 points can be inferred from sparse graphs with $\sim$1000
points, with negligible change in computation time. We empirically evaluate
our model on elastic solids, Newtonian fluids, Non-Newtonian fluids, Drucker-
Prager granular flows, and von Mises elastoplasticity. On these benchmarks,
our approach results in a 25$\times$ speedup compared to other neural network-
based physics simulators while delivering high-fidelity predictions of complex
physical systems and showing better performance on most benchmarks. The code
and the demos are provided at https://github.com/HrishikeshVish/GIOROM.

##  1 Introduction

Understanding the dynamics of physical systems is crucial in fields like
computational fluid mechanics, digital twins, graphics, and robotics. Accurate
simulations of temporal dynamics are essential for predicting turbulent ocean
behaviors, forecasting weather, and designing efficient physics-informed
robotic control algorithms, serving a broad spectrum of scientific and
engineering applications.

To numerically solve a system of temporal partial differential equations
(PDEs) governed by temporal dynamics, the system is discretized both
temporally and spatially. Temporal discretization breaks down the system’s
continuous evolution into discrete time steps. Similarly, spatial
discretization partitions the physical domain into a discrete set of spatial
points.

Traditional numerical methods, such as the material point method (MPM) [35],
the finite element method (FEM) [32], and smooth particle hydrodynamics (SPH)
[63] can accurately simulate temporal dynamics. However, these methods suffer
from high computational complexity and large memory requirements at finer
spatial-temporal resolutions, which are required for high-stake, high-fidelity
simulations [22, 57].

Machine learning-based methods, on the other hand, offer a faster alternative
to these traditional methods [3, 96, 18, 39]. However, data-driven approaches
are often trained on a particular discretization of the input and cannot
generalize to other discretizations during inference [42]. Models based on
graph neural networks [79] can, in theory, be applied to general meshes
(represented as graphs) but still face issues when the training and inference
resolutions differ significantly. Moreover, the computation cost scales poorly
with the size of the graph and aggregating the information between all
neighbors can be computationally prohibitive on dense graphs. One can, in
principle, use a sparse graph to speed up the computations. However, coarse
neighborhoods often lead to unphysical volume collapse111Volume collapse
happens when the average inter-particle distance in the predicted output is
significantly lower than in the reduced-order ground truth. during inference.
Neural operators [45, 56, 42], which learn the mappings between infinite-
dimensional function spaces, offer a principled way of dealing with input and
output data at arbitrary resolution. In particular, graph neural operators
[46, 47, 48] can be applied to graphs of any resolution by taking special care
of the construction and aggregation of neighborhoods.

We propose a new parameterization to the graph neural operator that allows it
to better capture spatial interactions of the particles within their
neighborhood. Integrating this parameterization into the geometry-informed
neural operator [48], we can efficiently learn the mapping between the
velocity and acceleration fields of 3D Lagrangian dynamics. Combining this
with state-of-the-art reduced-order techniques using neural fields allows us
to operate on sparse graphs while predicting at arbitrary high resolution
during inference and preserving the volume of the systems. By operating on
sparse graphs, our proposed model, GIOROM or Graph Interaction Operator for
Reduced-Order Modeling, is at least 25$\times$ faster than other neural
solvers. Moreover, using our proposed neural field, we can still predict the
full-order dynamics (with orders of magnitude higher resolution) with similar
or even better performance than the baselines. Our key contributions can be
summarized as follows:

  1. 1.

A framework to learn Lagrangian dynamics on highly sparse graphs: We propose a
neural operator-based model that can accurately and efficiently learn temporal
dynamics on very sparse graphs, achieving 25$\times$ speedup over other neural
physics solvers, while also delivering better performance on most benchmarks.

  2. 2.

Novel graph neural operator: Our model is based on an interaction operator
(IO) that significantly improves the performance of graph neural operators for
modeling spatial interactions.

  3. 3.

Arbitrary spatial evaluation using neural fields: We significantly accelerate
graph-based approaches by predicting the dynamics on reduced-order
discretizations and inferring the full system state using neural
representations of reduced-order models. We show results on full-order
representations up to two orders of magnitude larger than the reduced-order
representations.

##  2 Related Work

#### Neural physics solvers

Neural network-based solvers have shown great success in accelerating physical
simulations, including problems in fluid dynamics [79, 41, 86, 61, 82, 29],
solid mechanics [23, 14, 37], climate modeling [67], and robotics [65, 38].
Such approaches can be purely data-driven or physics-informed, i.e.,
leveraging an underlying PDE [77, 83, 78, 64]. Moreover, different
architectures have been proposed for the neural networks. Approaches based on
convolutional neural networks can be used to numerically solve systems of PDE
on fixed regular grids [43, 62, 84, 4]. For more general meshes, graph neural
networks (GNNs) have been proposed, e.g., in the context of mesh-based physics
[13, 71, 28, 21], Lagrangian dynamics [79], parametric PDEs [72], and rigid
body physics [40]. GNNs can efficiently capture spatial interactions between
particles. However, their time complexity scales with the size of the graph
since they require message-passing operations on every node. For finer
resolutions, this can be computationally prohibitive. Moreover, in their
standard formulation, they do not generalize to graphs that have significantly
different sizes than the ones seen during training.

#### Neural operator models

Neural operators are a class of discretization agnostic neural network
architectures that can generalize to arbitrary discretization of input data.
These architectures have been used in solving parametric PDEs [56, 47, 44, 3,
74, 55, 42, 73, 50, 88, 81, 25], fluid dynamics [20, 89, 70], protein
interactions [53, 52, 19], 3D physics [94, 91, 12, 75, 30, 75], weather
modeling [11, 67], robotics [10, 69], and computer vision [27, 76, 87]. In our
work, we propose a new parameterization for the graph neural operator to
efficiently capture spatial dynamics in Lagrangian systems and generalize to
different discretizations of reduced-order inputs.

#### Reduced-order model

Reduced-order models (ROMs) simplify high-dimensional dynamic systems by
projecting them onto a lower-dimensional manifold, resulting in faster and
less expensive computations [9, 31, 68]. These methods gain computational
efficiency by simulating a subset of the original spatial samples [1]. Recent
neural field-based ROMs [66, 95, 90, 17] demonstrate the ability to train a
discretization-agnostic low-dimensional representation, allowing the trained
model to generalize over various geometric discretizations. We extend these
ROMs by using a neural operator to compute the dynamics of the spatial
samples, such that the model is a non-intrusive and discretization-agnostic
ROM system. This combination creates a machine-learning model that generalizes
across geometries while being significantly faster than traditional models.

##  3 Method: Problem setting and reduced-order modeling

###  3.1 Problem setting

In this work, we are concerned with a spatial-temporal dynamical system that
can be described in the form

$\displaystyle{\boldsymbol{\mathcal{J}}}({\boldsymbol{\phi}},{\boldsymbol{\nabla}}{\boldsymbol{\phi}},{\boldsymbol{\nabla}}^{2}{\boldsymbol{\phi}},\ldots,\dot{{\boldsymbol{\phi}}},\ddot{{\boldsymbol{\phi}}},\ldots)={\boldsymbol{0}},\quad{\boldsymbol{\phi}}({\boldsymbol{X}},t):\Omega\times\mathcal{T}\to\mathbb{R}^{d}\
,$ (1)

where ${\boldsymbol{\phi}}$ represents a multidimensional vector field that
depends on both space and time. The symbols ${\boldsymbol{\nabla}}$ and
$\dot{(\cdot)}$ signify the spatial gradient and time derivative,
respectively. Here, $\Omega\subset\mathbb{R}^{d}$ and
$\mathcal{T}\subset\mathbb{R}$ denote the spatial and temporal domains,
respectively. While ${\boldsymbol{\phi}}$ can generally be any spatial-
temporal vector field, we focus on the deformation map arising from continuum
mechanics [24], i.e.,
${\boldsymbol{x}}={\boldsymbol{\phi}}({\boldsymbol{X}},t)$. Its first and
second temporal derivatives denote “acceleration”
${\boldsymbol{a}}=\ddot{{\boldsymbol{\phi}}}$ and “velocity”
${\boldsymbol{v}}=\dot{{\boldsymbol{\phi}}}$, respectively. After introducing
a temporal discretization $\\{t_{n}\\}_{n=0}^{{N}}$, we aim to solve for a
sequence of spatial functions
$\\{{\boldsymbol{\phi}}_{t_{n}}({\boldsymbol{X}})\\}_{n=0}^{{N}}$, where
${\boldsymbol{\phi}}_{t_{n}}({\boldsymbol{X}})={\boldsymbol{\phi}}({\boldsymbol{X}},t_{n})$.

###  3.2 Discretization-agnostic reduced-order modeling

To represent these continuous spatial functions in a numerical simulation, one
needs to introduce spatial discretizations. Let
$\\{{\boldsymbol{X}}^{j}\\}_{j=1}^{{P}}$ be the ${P}$-point discretization of
the spatial domain $\Omega$, where ${P}$ is the number of full-order spatial
samples. Traditional full-order numerical PDE solvers directly operate on
these spatial discretizations [32]. Due to the fact that ${P}$ needs to be
very large for high-fidelity simulations, these full-order solvers are
extremely expensive and are not suitable for time-critical applications [57].

To address this computation challenge, reduced-order modeling (ROM) approaches
leverage data priors to evolve the spatial functions over time using only a
small number of spatial samples [15]. Concretely, ROM techniques leverage a
${Q}$-point discretization of the spatial domain
$\\{{\boldsymbol{X}}^{k}\\}_{k=1}^{{Q}}$, where ${Q}\ll{P}$. In particular, by
leveraging neural fields and projection-based ROM, the work by [16] proposes a
technique that can infer the continuous spatial function at arbitrary spatial
locations from just a few spatial samples, i.e., the ${Q}$-point
discretization. Section 4.3 lists the details of the upsampling algorithm. To
evolve these ${Q}$-point discretizations over time, we seek a mapping between
$\\{{\boldsymbol{\phi}}_{t_{n}}^{k}\\}_{k=1}^{{Q}}$ and
$\\{{\boldsymbol{\phi}}_{t_{n+1}}^{k}\\}_{k=1}^{{Q}}$, where
${\boldsymbol{\phi}}_{t_{n}}^{k}=\mathbf{X}_{t_{n}}^{k}={\boldsymbol{\phi}}_{t_{n}}(\mathbf{X}^{k})$.

###  3.3 Time integration

To further break down this mapping, we first leverage an explicit Euler time
integrator [2] with step-size $\Delta t$,

$\displaystyle{\boldsymbol{\phi}}_{t_{n+1}}^{j}={\boldsymbol{\phi}}_{t_{n}}^{j}+\Delta
t\,\dot{{\boldsymbol{\phi}}}_{t_{n}}^{j}$ (2)
$\displaystyle\dot{{\boldsymbol{\phi}}}_{t_{n+1}}^{j}=\dot{{\boldsymbol{\phi}}}_{t_{n}}^{j}+\Delta
t\,\ddot{{\boldsymbol{\phi}}}_{t_{n}}^{j}$ (3)

As such, the one and only unknown in the equation above is the acceleration
$\mathbf{A}_{t_{n}}^{j}=\ddot{{\boldsymbol{\phi}}}_{t_{n}}^{j}$, which is
necessary for computing the velocity
$\mathbf{V}_{t_{n+1}}^{j}=\dot{{\boldsymbol{\phi}}}_{t_{n+1}}^{j}$. We propose
to predict the acceleration field from the current and past velocity fields
via neural operators.

##  4 Method: Operator learning

Neural operators are a class of machine learning models that learn to map
functions in infinite-dimensional function spaces using a finite collection of
discretized input-output pairs. In particular, we want to learn the mapping
from the current and past velocity fields
$(\dot{{\boldsymbol{\phi}}}_{t_{i}})_{i=n-w}^{n}$ to the current acceleration
field $\ddot{{\boldsymbol{\phi}}}_{t_{n}}$. The hyperparameter $w$ defines the
time window given by the past $w$ time steps. In practice, we use a finite
collection of spatially discretized input-output pairs, as, e.g., provided by
the $Q$-point discretization.

Our proposed neural operator leverages the geometry-informed neural operator
(GINO) [48] as a backbone. Its encoder and decoder are given by integral
transforms, parametrized as graph neural operators. This allows to map from an
arbitrary discretization of the input and output fields to a uniform
discretization in the latent space. We can then leverage a Fourier neural
operator to efficiently learn global interactions via spectral convolutions in
the latent space. See Appendix B.1 for further details. To efficiently learn
local interactions, we propose residual connections combined with a novel
interaction operator to pre- and postprocess the input and output fields. This
operator represents a discretization-agnostic adaptation of message-passing
architectures, such as the Interaction Network proposed in [5]. As described
in Section 3.3, we use the predicted acceleration field to obtain the
deformation field ${\boldsymbol{\phi}}_{t_{n+1}}^{k}$ at the $Q$-point
discretization $\\{{\boldsymbol{X}}^{k}\\}_{k=1}^{{Q}}$ via Euler steps.
Finally, we leverage neural representations of reduced-order models [16] to
obtain the full deformation field ${\boldsymbol{\phi}}_{t_{n+1}}$. We refer to
Figure 1 for an overview and provide further details on the aforementioned
building blocks in the following.

Figure 1: The overall architecture of GIOROM. The neural operator
$\mathcal{G}_{\theta}$ predicts the acceleration of a Lagrangian system
$\mathbf{A}_{t_{k}}$ at time $t_{k}$ from the past $w$ velocity instances
$\mathbf{V}_{t_{k-w:k}}$. The positions are derived through Euler integration.
The neural field is used to evaluate the deformation field at arbitrary
locations.

###  4.1 Geometry-informed neural operator

Inter-particle spatial interactions are typically encoded using graph neural
networks in many neural physics simulators. However, graph neural networks
(GNNs) that use nearest-neighbor connectivity are not discretization-agnostic
and degenerate to a pointwise operator [42]. In contrast, a graph neural
operators (GNO) consider a radius graph, where a point $x$ is connected to all
points within a ball $B_{r}(x)$ of a certain radius $r$ [46, 47, 48]. This can
be understood as a discretization of an integral transform

$\operatorname{GNO}(v)(x)=\int_{B_{r}(x)}\kappa_{\theta}(x,y,v(y))dy,$ (4)

where $v$ denotes a suitable input function and $\kappa$ is a learnable
kernel, parametrized by a neural network. If the input function $v$ is
discretized at points $y_{i}\in B_{r}(x)$, the integral transform in (4) can
be approximated by

$\operatorname{GNO}(v)(x)\approx\sum_{y_{i}\in
B_{r}(x)}\kappa_{\theta}(x,y_{i},v(y_{i}))\Delta y_{i},$ (5)

where $\Delta y_{i}$ are suitable integration weights.

While geometries in the physical domain are complex and irregular, we follow
[48] and efficiently learn the temporal dynamics on a uniform grid using a
Fourier neural operator (FNO) [45]. To switch between the given discretization
and a uniform grid, we note that the output function $\operatorname{GNO}(v)$
in (5) can be evaluated on points $x$ different from the discretization points
$\\{y_{i}\\}_{i}$. Leveraging this flexibility, the GINO architecture is thus
given by

$\operatorname{GINO}=\operatorname{GNO}_{\mathrm{dec}}\circ\operatorname{FNO}_{L}\circ\dots\circ\operatorname{FNO}_{1}\circ\operatorname{GNO}_{\mathrm{enc}},$
(6)

where $L$ denotes the number of FNO layers and each layer has its own
learnable parameters.

Each FNO layer uses spectral convolutions to efficiently model global spatial
dependencies on the regular latent grid. More precisely, we use the
convolution theorem to define the spectral convolution as

$\mathcal{K}(v)=\mathcal{F}^{-1}(R_{\theta}\,\mathcal{F}(v)),$ (7)

where $\mathcal{F}$ and $\mathcal{F}^{-1}$ denote the mapping to the Fourier
coefficients and its inverse, and $R_{\theta}$ are learnable parameters
(corresponding to a fixed number of Fourier coefficients of the kernel). Since
the input function $v$ is given on a regular grid, we can efficiently
approximate $\mathcal{F}$ and $\mathcal{F}^{-1}$ using the fast Fourier
transform (FFT). The FNO is then defined as

$\operatorname{FNO}(v)(x)=\varrho(W_{\theta}(v(x))+\mathcal{K}(v)(x)),$ (8)

where $\varrho$ is a pointwise non-linearity and $W_{\theta}$ is a learnable
linear layer.

###  4.2 Interaction neural operator

To improve the performance of GINO, we add additional layers inspired by
successful message-passing GNNs, such as the interaction network proposed in
[5]; see Appendix C for further details. However, crucially, we adhere to the
operator framework, such that our resulting interaction operator is agnostic
to the discretization of the input function.

In particular, we define a general interaction operator (IO) as

$\operatorname{IO}_{k}(v)(x)=f_{\theta}\Big{(}v(x),\int_{B_{r}(x)}\kappa_{\theta}(k(x,y),v(x),v(y))dy\Big{)}.$
(9)

In contrast to existing GNOs as in (4), we let the kernel $\kappa_{\theta}$ in
(9) depend on $v(x)$ and an additional function $k$ representing edge
features. Moreover, we allow for residual connections through $f_{\theta}$,
which is parametrized by a neural network. Our proposed operator then takes
the form

$\mathcal{G}_{\theta}(v)=(\operatorname{IO}_{k,\mathrm{dec}}\circ\operatorname{GINO}\circ\operatorname{IO}_{h,\mathrm{enc}})(v).$
(10)

We choose $h(x,y)=g_{\theta}((x-y)/r)$ and
$k(x,y)=h(x,y)+\kappa_{\theta,\mathrm{enc}}(h(x,y),v(x),v(y))$, where
$g_{\theta}$ is parametrized by a neural network and
$\kappa_{\theta,\mathrm{enc}}$ is the kernel of
$\operatorname{IO}_{h,\mathrm{enc}}$. This effectively creates a residual
connection between the interaction operators. We can discretize (9) similar to
(5) but require the evaluation point $x$ to be included in the set of
discretization points $\\{y_{i}\\}_{i}$ at which we know the value of the
input function $v$ and the edge features $k$. This is possible since we use
the same discretization for the input and output functions of the interaction
operators. We empirically observe that this parametrization helps to model
spatial information about the interactions between the graph nodes,
significantly improving performance on complex geometries; see Figure 6.

Figure 2: Visual comparison of GIOROM prediction with other neural solvers.
GIOROM achieves higher accuracy than other neural solvers across various
materials and problem dimensions. In addition to the accuracy advantage,
GIOROM is also significantly faster, as shown in table 4.

The discretized version of $\mathcal{G}_{\theta}$ takes a matrix
$\mathbf{V}_{t_{n-w:n}}\in\mathbb{R}^{w\times{Q}\times d}$ of $w$ consecutive
velocities defined over ${Q}$ samples $\\{{\boldsymbol{X}}^{k}\\}_{k=1}^{{Q}}$
in $d$-dimensional space for $w$ time steps as inputs222In practice, we use an
embedding of the velocities as well as the distances of each particle to the
boundaries and the particle type as input to the model to improve
performance.. It then outputs the pointwise accelerations
$\mathbf{A}_{t_{n}}\in\mathbb{R}^{{Q}\times d}$ of the samples at the current
time step. Since $\mathcal{G}_{\theta}$ is a neural operator, the output is
consistent across different numbers of particles ${Q}$. Using the predicted
accelerations $\mathbf{A}_{t_{n}}$, we can leverage (2) and (3) to compute the
next time-step positions
${\boldsymbol{\phi}}_{t_{n+1}}({\boldsymbol{X}}^{k})$.

###  4.3 Full-order inference using neural fields

Equipped with the next time-step positions
${\boldsymbol{\phi}}_{t_{n+1}}({\boldsymbol{X}}^{k})$ at the reduced-order
${Q}$-point discretizations $\\{{\boldsymbol{X}}^{k}\\}_{k=1}^{{Q}}$, we will
compute the next time-step positions
${\boldsymbol{\phi}}_{t_{n+1}}({\boldsymbol{X}}^{j})$ at the full-order
${P}$-point discretizations $\\{{\boldsymbol{X}}^{j}\\}_{j=1}^{{P}}$. To do
so, we leverage a neural representation of projection-based reduced-order
models (ROM) [8]. ROM assumes that
${\boldsymbol{\phi}}_{t_{n+1}}({\boldsymbol{X}})$ can be represented as a
weighted sum of a small number of basis functions $\\{U^{i}\\}_{i=1}^{r}$ with
weights $\\{q_{t_{n+1}}^{i}\\}_{i=1}^{r}$:

$\displaystyle{\boldsymbol{\phi}}_{t_{n+1}}({\boldsymbol{X}})={\boldsymbol{X}}+\mathbf{U}({\boldsymbol{X}})\mathbf{q}_{t_{n+1}}.$
(11)

Here, the matrix $\mathbf{U}({\boldsymbol{X}})$ and the vector
$\mathbf{q}_{t_{n+1}}$ are vectorized versions of
$\\{U^{i}({\boldsymbol{X}})\\}_{i=1}^{r}$ and
$\\{q_{t_{n+1}}^{i}\\}_{i=1}^{r}$, respectively.

The basis functions are implemented using neural fields whose weights are
learnable. We follow the same training procedure as described in [16]. After
training, the spatially-dependent $\mathbf{U}({\boldsymbol{X}})$ stays fixed
over time while the weights $\\{q_{t_{n+1}}^{k}\\}_{k=1}^{r}$ for each basis
change at each time step.

When we have the function
${\boldsymbol{\phi}}_{t_{n+1}}({\boldsymbol{X}}^{k})$ for the sub-sampling
${\boldsymbol{X}}^{k}\in\\{{\boldsymbol{X}}^{k}\\}_{k=1}^{{Q}}$, we can
calculate $\mathbf{q}_{t_{n+1}}$ by solving the least squares problem:

$\displaystyle\min_{\mathbf{q}_{t_{n+1}}}\sum_{k=1}^{{Q}}\|{\boldsymbol{\phi}}_{t_{n+1}}({\boldsymbol{X}}^{k})-({\boldsymbol{X}}^{k}+\mathbf{U}({\boldsymbol{X}}^{k})\mathbf{q}_{t_{n+1}})\|_{2}^{2}.$
(12)

After $\mathbf{q}_{t_{n+1}}$ is obtained, we are able to calculate
${\boldsymbol{\phi}}_{t_{n+1}}({\boldsymbol{X}}^{j})$ by:

$\displaystyle{\boldsymbol{\phi}}_{t_{n+1}}({\boldsymbol{X}}^{j})={\boldsymbol{X}}^{j}+\mathbf{U}({\boldsymbol{X}}^{j})\mathbf{q}_{t_{n+1}},\quad\forall{\boldsymbol{X}}^{j}\in\\{{\boldsymbol{X}}^{j}\\}_{j=1}^{{P}}.$
(13)

##  5 Experiments

###  5.1 Dataset

We trained our model on simulations generated by SPlisHSPlasH [7] and NCLAW
[58] simulators. The NCLAW dataset contains multiple materials (jelly,
plasticine, sand, and water) and 15 trajectories with different initial
conditions, each with 2000 time steps per material and 4 geometries per
material. We also include 2D simulations generated by [71], with 1000
trajectories per material - sand, jelly, water, and multi-material. Additional
details are given in Appendix E. The dataset contains multiple physical
materials - elasticity, plasticine (von Mises), granular media (Drucker Prager
sand flows), Newtonian fluids, multi-material, and non-Newtonian fluids. The
datasets are generated using 3 types of simulators - MPM, SPH, and FEM.

###  5.2 Model setup and hyperparameters

#### Sparse graph construction

To model the physics in reduced-order space, it is very important to ensure
that the volume is preserved. This helps restore the full-graph post inference
without any loss in accuracy. The key factor is to ensure that there is a
single connected component for every object on which the model learns. As
shown in Figure 3, if this property is violated, the model treats each
component as a separate object, thus leading to volume collapse. The graph is
constructed by increasing the neighborhood radius until there is a single
connected component.

Figure 3: Our method operates on a reduced-order graph. (a) depicts the graph
with all the points, each connected to all points within a radius of 0.015.
Operating on this full-order graph is expensive. (b) We obtain several
disconnected components when subsampling the graph and picking the radius too
small, which can lead to volume collapse in simulations. (c) We improve the
prediction of the dynamics by choosing a radius that induces a single
connected component.

#### Loss function

The model predicts per-particle acceleration from a sequence of past
velocities of the samples. The loss is defined as the mean squared error
between the predicted acceleration and the ground-truth acceleration in the
simulation sequence. For a consecutive pair of positions $\mathbf{X}_{t_{n}}$
and $\mathbf{X}_{t_{n+1}}$, with corresponding velocities $\mathbf{V}_{t_{n}}$
and $\mathbf{V}_{t_{n+1}}$, the corresponding acceleration is defined as
$\hat{\mathbf{A}}_{t_{n}}=(\mathbf{V}_{t_{n+1}}-\mathbf{V}_{t_{n}})/\Delta t$.
The loss is thus given as

$L(\theta)=\|\mathcal{G}_{\theta}(\mathbf{V}_{t_{n-w:n}})-\hat{\mathbf{A}}_{t_{n}}\|_{2}^{2}$
(14)

Additional model hyperparameters can be found in Appendix C.2.

Table 1: GIOROM model outperforms various graph-based and neural operator-
based models for different types of materials.

Material | Noise | Model | 1-step Loss | Rollout-Loss
---|---|---|---|---
Plasticine | 1e-4 | GIOROM | 1.2e-8 | 0.0094
|  | GNS | 6.8e-8 | 0.0098
|  | GINO | 3.2e-7 | 0.088
|  | GAT | 9.7e-8 | 0.024
Elasticity | 3e-4 | GIOROM | 1.07e-8 | 0.0034
|  | GNS | 8.73e-9 | 0.0069
|  | GINO | 1.54e-8 | 0.687
|  | GAT | 8.61e-9 | 0.011
Water | 3e-4 | GIOROM | 7.05e-9 | 0.0074
|  | GNS | 5.21e-9 | 0.0079
|  | GINO | 3.3e-7 | 0.10
|  | GAT | 1.0e-8 | 0.098
Sand | 3e-4 | GIOROM | 7.03e-9 | 0.0063
|  | GNS | 6.34e-9 | 0.0071
|  | GINO | 2.1e-7 | 0.583
|  | GAT | 9.49e-9 | 0.062
Jelly | 3e-4 | GIOROM | 1.5e-8 | 0.0019
|  | GNS | 4.37e-9 | 0.0023
|  | GINO | 3.0e-7 | 0.11
|  | GAT | 7.9e-8 | 0.038
Multi-material | 3e-4 | GIOROM | 7.0e-9 | 0.0056
|  | GNS | 4.37e-9 | 0.013
|  | GINO | 9.7e-8 | 0.749
|  | GAT | 3.9e-8 | 0.024

##  6 Results

#### Baselines

We compare our model against three graph neural network-based models and a
neural operator-based model. For the graph neural networks, we consider GNS
[79], which uses a series of message-passing layers in an encoder-decoder
fashion, MeshGraphNet [71], a slightly modified version of GNS, and GAT [51].
For the neural operators, we compare against GINO [48], a combination of graph
neural operators and a Fourier neural operator [45].

#### Scale of complex simulations

Our model is able to simulate various materials across different geometries
and initial conditions, as shown in Table 7 and Figure 7 (2D) and Figure 8
(3D) in Appendix G. In all of these cases, the physics was simulated in a
reduced-order setting on approximately 500-2000 particles, and the rest of the
particles were recovered using neural fields with nearly no additional cost in
terms of computation and time. The full-order settings were as dense as 95K
particles, showing that the model was able to simulate the dynamics accurately
using only 2% of the full-order degrees of freedom.

#### Comparison to other neural solvers

Figure 2 highlights the performance of our model against other neural solvers.
We show in tables 3 and 4 that our model has a comparable parameter count with
GAT and GNS models but is significantly faster during inference on graphs with
the same connectivity. This speedup is due to fewer layers of message-passing
and more layers of spectral convolution. It can be seen in the table that on
denser graphs, the speedup of GIOROM is at least 3x compared to GNS and at
least 8x compared to GAT, for inputs of the same size. During inference,
GIOROM predicts the physics by operating on highly sparse graphs. To render a
point cloud of 7056 degrees of freedom, it requires only 1776 points, thus
achieving a 25$\times$ speedup compared to GNS. The up-scaled point clouds are
rendered using neural fields and have negligible computation time for
upscaling, which is shown in Table 2. From Table 1 and Table 6 in Appendix G,
it can be observed that GIOROM has a lower rollout-loss compared to all other
neural solvers in the full-order setting.

Figure 4: Discretization agnosticism. Our approach yields consistent results
across discretizations. Table 2: Contrasting the inference times (in seconds)
for highly dense point clouds up-sampled from highly sparse graphs (1776
points).

Time Steps | Rollout Size | Rollout Time (s) | Full Order Size | Upscale Time (s)
---|---|---|---|---
200 | 1776 | 4.4 | 7,000 | 5e-5
200 | 1776 | 4.4 | 40,000 | 3e-4
200 | 1776 | 4.4 | 60,000 | 8e-4
200 | 1776 | 4.4 | 100,000 | 9e-4

Table 3: Contrasting the change in computation time with the increase in
connectivity radius for a graph with 7056 points. The times shown represent
the overall time needed to infer all 200 time steps.

Model | Time Steps | Number of Spatial Points | Connectivity Radius
---|---|---|---
0.040 | 0.050 | 0.060 | 0.070 | 0.080 | 0.090 | 0.100
GIOROM | 200 | 7056 points | 28.1s | 40.4s | 52.7s | 75.3s | 95.8s | 107.5s | 115.3s
GNS | 200 | 7056 points | 43.5s | 73.5s | 111.6s | 162s | 226.2s | 305.9s | 386.0s
GAT | 200 | 7056 points | 146.5s | 236.5s | 394.2s | 532.8s | 645.2s | 733.8s | 812.5s

Table 4: Contrasting the change in computation time with an increase in graph
size at a fixed radius of 0.060. The times shown represent the overall time
needed to infer all 200 time steps.

Model | Parameters | Connectivity | Material | Time Steps | Graph Size
---|---|---|---|---|---
1776 points | 4143 points | 5608 points | 7056 points
GIOROM | 1,721,374 | 0.060 | Plasticine | 200 | 4.4s | 17.7s | 32.5s | 52.7s
GNS | 1,592,987 | 0.060 | Plasticine | 200 | 7.8 s | 38.3s | 68.7s | 111.6s
GAT | 1,999,003 | 0.060 | Plasticine | 200 | 71.1s | 153.4s | 295.3s | 394.2s

##  7 Ablation studies

The following section discusses all the major ablations we performed on our
model. For additional ablations, refer to Appendix H.

###  7.1 Generalization to unseen and complex geometries

We test the model in two settings, firstly, on an upsampled denser graph of
the same geometry. In the second case, we test the model on upsampled graphs
of previously unseen geometries, with greater asymmetry and mean surface
curvature compared to the graphs it was trained on, as shown in Table 8. For
3D water and sand trajectories, the model was trained on armadillo and blub
and was shown to generalize to bunny and spot, as shown in Figure 8. For
elasticity, as seen in Figure 8, the armadillo was the previously unseen
geometry.

###  7.2 Key architecture choices

#### Cost-accuracy tradeoff on graph connectivity radius

We experiment with different radii during training. It can be observed that
increasing the connectivity radius improves the accuracy of the rollout in all
of the simulations. Furthermore, on sub-sampled graphs, a larger radius is
necessary to ensure that the network learns on a single connected component.
The presence of stray components within the graph causes a volume collapse
during rollout.

#### Graph radius and viscosity

During inference, the connectivity radius plays an important role in modifying
the material properties of the fluid. Increasing the connectivity introduces
rigidity into the graph structure, increasing the viscosity of the fluid. This
is experimentally shown in Table 11, in Appendix H, which compares the mean
inter-particle distance and the mean particle velocity against the
connectivity radius. Figure 5 highlights the changes in viscosity during
inference.

Figure 5: The figure highlights the increase in viscosity as the radius
increases, due to a larger neighborhood, (a) with a radius of 0.010, (b)
0.015, (c) 0.025, (d) 0.040, (e) 0.055.

###  7.3 Discretization invariance and generalizability

We perform studies to show the discretization invariance of the model in
different settings. The three settings we consider are the effects of latent
grid dimensions, sub-sampling during training, and super-sampling during
inference. Figure 4 highlights how the model can accurately infer the
trajectories on different discretizations of the same input.

#### Effects of sub-sampling on training

The model was trained on sub-sampled point clouds at different sampling rates
(0.5$\times$, 0.3$\times$, 0.1$\times$, 0.01$\times$). The reason for doing
this is that sparse graphs have lower computational overhead compared to full-
order settings. Table 5 shows the training performance of the model on
different sampling rates with both fixed and variable radii.

Table 5: The tables denote the invariance to sampling rate during training, as
long as appropriate radii are chosen. When the radii are fixed, the model
fails to accurately model the acceleration.

Sampling Rate | Radius | Loss
---|---|---
1.0$\times$ | 0.015 | 0.0090
0.5$\times$ | 0.035 | 0.0098
0.25$\times$ | 0.070 | 0.0086
0.10$\times$ | 0.095 | 0.0087
0.01$\times$ | 0.120 | 0.0085

Sampling Rate | Radius | Loss
---|---|---
1.0$\times$ | 0.015 | 0.0090
0.5$\times$ | 0.015 | 0.010
0.25$\times$ | 0.015 | 0.050
0.10$\times$ | 0.015 | 0.175
0.01$\times$ | 0.015 | 1.248

#### Zero-shot super-resolution during inference

Zero-shot super-resolution refers to the model’s ability to generalize to
super-sampled dense graphs after being trained on sparse graphs. We show in
Table 9 that the model maintains the same rollout error during inference on
previously unseen super-sampled graphs.

##  8 Conclusion

In conclusion, our proposed model, GIOROM, marks a substantial advance in
simulating Lagrangian dynamics across various fields. We showcase results on
2D and 3D simulations of several materials such as Newtonian fluids, non-
Newtonian fluids, elastic materials, plasticine, and granular material.
Utilizing a reduced-order modeling approach on sparse graphs, GIOROM is up to
25× faster than previous neural network-based physics solvers while achieving
similar or better performance. Moreover, our neural-operator-based model
generalizes well across different initial conditions, velocities,
discretizations, and geometries.

Despite its promising performance, GIOROM has several limitations that warrant
further exploration. For instance, the choice of noise parameters can
significantly alter the rollout accuracy of the model. Additionally, while
GIOROM is capable of generalizing across different settings, like many machine
learning and reduced-order methods, it struggles with extreme out-of-
distribution scenarios [46, 17]. Moreover, while GIOROM is primarily designed
for continuous systems, future research might explore mechanisms to explicitly
handle discontinuities [6, 26].

## References

  * [1] S. S. An, T. Kim, and D. L. James.  Optimizing cubature for efficient integration of subspace deformations.  ACM transactions on graphics (TOG), 27(5):1–10, 2008.
  * [2] U. M. Ascher and L. R. Petzold.  Computer methods for ordinary differential equations and differential-algebraic equations.  SIAM, 1998.
  * [3] K. Azizzadenesheli, N. Kovachki, Z. Li, M. Liu-Schiaffini, J. Kossaifi, and A. Anandkumar.  Neural operators for accelerating scientific simulations and design.  Nature Reviews Physics, pages 1–9, 2024.
  * [4] F. Bamer, D. Thaler, M. Stoffel, and B. Markert.  A monte carlo simulation approach in non-linear structural dynamics using convolutional neural networks.  Frontiers in Built Environment, 7, 2021.
  * [5] P. Battaglia, R. Pascanu, M. Lai, D. Jimenez Rezende, et al.  Interaction networks for learning about objects, relations and physics.  Advances in neural information processing systems, 29, 2016.
  * [6] Y. Belhe, M. Gharbi, M. Fisher, I. Georgiev, R. Ramamoorthi, and T.-M. Li.  Discontinuity-aware 2d neural fields.  ACM Transactions on Graphics (TOG), 42(6):1–11, 2023.
  * [7] J. Bender and D. Koschier.  Divergence-free smoothed particle hydrodynamics.  In Proceedings of the 14th ACM SIGGRAPH/Eurographics symposium on computer animation, pages 147–155, 2015.
  * [8] P. Benner, S. Gugercin, and K. Willcox.  A survey of projection-based model reduction methods for parametric dynamical systems.  SIAM review, 57(4):483–531, 2015.
  * [9] G. Berkooz, P. Holmes, and J. L. Lumley.  The proper orthogonal decomposition in the analysis of turbulent flows.  Annual review of fluid mechanics, 25(1):539–575, 1993.
  * [10] R. Bhaskara, H. Viswanath, and A. Bera.  Trajectory prediction for robot navigation using flow-guided markov neural operator.  arXiv preprint arXiv:2309.09137, 2023.
  * [11] S. Bire, B. Lütjens, K. Azizzadenesheli, A. Anandkumar, and C. N. Hill.  Ocean emulation with fourier neural operators: Double gyre.  Authorea Preprints, 2023.
  * [12] B. Bonev, T. Kurth, C. Hundt, J. Pathak, M. Baust, K. Kashinath, and A. Anandkumar.  Spherical fourier neural operators: Learning stable dynamics on the sphere.  In International conference on machine learning, pages 2806–2823. PMLR, 2023.
  * [13] Y. Cao, M. Chai, M. Li, and C. Jiang.  Efficient learning of mesh-based physical simulation with bsms-gnn.  arXiv preprint arXiv:2210.02573, 2022.
  * [14] G. Capuano and J. J. Rimoli.  Smart finite elements: A novel machine learning application.  Computer Methods in Applied Mechanics and Engineering, 345:363–381, 2019.
  * [15] K. Carlberg, C. Farhat, J. Cortial, and D. Amsallem.  The gnat method for nonlinear model reduction: effective implementation and application to computational fluid dynamics and turbulent flows.  Journal of Computational Physics, 242:623–647, 2013.
  * [16] Y. Chang, P. Y. Chen, Z. Wang, M. M. Chiaramonte, K. Carlberg, and E. Grinspun.  Licrom: Linear-subspace continuous reduced order modeling with neural fields.  In SIGGRAPH Asia 2023 Conference Papers, SA ’23, New York, NY, USA, 2023. Association for Computing Machinery.
  * [17] P. Y. Chen, J. Xiang, D. H. Cho, Y. Chang, G. A. Pershing, H. T. Maia, M. M. Chiaramonte, K. T. Carlberg, and E. Grinspun.  CROM: Continuous reduced-order modeling of PDEs using implicit neural representations.  In The Eleventh International Conference on Learning Representations, 2023.
  * [18] S. Cuomo, V. S. Di Cola, F. Giampaolo, G. Rozza, M. Raissi, and F. Piccialli.  Scientific machine learning through physics–informed neural networks: Where we are and what’s next.  Journal of Scientific Computing, 92(3):88, 2022.
  * [19] G. Dharuman, L. Ward, H. Ma, P. V. Setty, O. Gokdemir, S. Foreman, M. Emani, K. Hippe, A. Brace, K. Keipert, et al.  Protein generation via genome-scale language models with bio-physical scoring.  In Proceedings of the SC’23 Workshops of The International Conference on High Performance Computing, Network, Storage, and Analysis, pages 95–101, 2023.
  * [20] P. C. Di Leoni, L. Lu, C. Meneveau, G. E. Karniadakis, and T. A. Zaki.  Neural operator prediction of linear instability waves in high-speed boundary layers.  Journal of Computational Physics, 474:111793, 2023.
  * [21] M. Fortunato, T. Pfaff, P. Wirnsberger, A. Pritzel, and P. Battaglia.  Multiscale meshgraphnets.  arXiv preprint arXiv:2210.00612, 2022.
  * [22] M. Gao, X. Wang, K. Wu, A. Pradhana, E. Sifakis, C. Yuksel, and C. Jiang.  GPU optimization of material point methods.  ACM Transactions on Graphics (TOG), 37(6):1–12, 2018.
  * [23] A. R. Geist and S. Trimpe.  Structured learning of rigid-body dynamics: A survey and unified view from a robotics perspective.  GAMM-Mitteilungen, 44(2):e202100009, 2021.
  * [24] O. Gonzalez and A. M. Stuart.  A first course in continuum mechanics, volume 42.  Cambridge University Press, 2008.
  * [25] S. Goswami, A. Bora, Y. Yu, and G. E. Karniadakis.  Physics-informed deep neural operator networks.  In Machine Learning in Modeling and Simulation: Methods and Applications, pages 219–254. Springer, 2023.
  * [26] S. Goswami, M. Yin, Y. Yu, and G. E. Karniadakis.  A physics-informed variational deeponet for predicting crack path in quasi-brittle materials.  Computer Methods in Applied Mechanics and Engineering, 391:114587, 2022.
  * [27] J. Guibas, M. Mardani, Z. Li, A. Tao, A. Anandkumar, and B. Catanzaro.  Adaptive fourier neural operators: Efficient token mixers for transformers.  arXiv preprint arXiv:2111.13587, 2021.
  * [28] X. Han, H. Gao, T. Pfaff, J.-X. Wang, and L.-P. Liu.  Predicting physics in mesh-reduced space with temporal attention.  arXiv preprint arXiv:2201.09113, 2022.
  * [29] Z. Hao, C. Su, S. Liu, J. Berner, C. Ying, H. Su, A. Anandkumar, J. Song, and J. Zhu.  DPOT: Auto-regressive denoising operator transformer for large-scale pde pre-training.  arXiv preprint arXiv:2403.03542, 2024.
  * [30] J. He, S. Koric, D. Abueidda, A. Najafi, and I. Jasiuk.  Geom-deeponet: A point-cloud-based deep operator network for field predictions on 3d parameterized geometries.  arXiv preprint arXiv:2403.14788, 2024.
  * [31] P. Holmes, J. L. Lumley, G. Berkooz, and C. W. Rowley.  Turbulence, coherent structures, dynamical systems and symmetry.  Cambridge university press, 2012.
  * [32] T. J. Hughes.  The finite element method: linear static and dynamic finite element analysis.  Courier Corporation, 2012.
  * [33] S. Janny, M. Nadri, J. Digne, and C. Wolf.  Space and time continuous physics simulation from partial observations.  In International Conference on Learning Representation, Vienna, Austria, May 2024.
  * [34] J. Jeon, J. Lee, R. Vinuesa, and S. J. Kim.  Residual-based physics-informed transfer learning: A hybrid method for accelerating long-term cfd simulations via deep learning.  International Journal of Heat and Mass Transfer, 220:124900, 2024.
  * [35] C. Jiang, C. Schroeder, J. Teran, A. Stomakhin, and A. Selle.  The material point method for simulating continuum materials.  ACM SIGGRAPH 2016 Courses, pages 1–52, 2016.
  * [36] R. Jiang, P. Y. Lu, E. Orlova, and R. Willett.  Training neural operators to preserve invariant measures of chaotic attractors.  Advances in Neural Information Processing Systems, 36, 2024.
  * [37] H. Jin, E. Zhang, and H. D. Espinosa.  Recent advances and applications of machine learning in experimental solid mechanics: A review.  Applied Mechanics Reviews, 75(6):061001, 2023.
  * [38] B. Kaczmarski, A. Goriely, E. Kuhl, and D. E. Moulton.  A simulation tool for physics-informed control of biomimetic soft robotic arms.  IEEE Robotics and Automation Letters, 8(2):936–943, 2023.
  * [39] G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang.  Physics-informed machine learning.  Nature Reviews Physics, 3(6):422–440, 2021.
  * [40] J. Kneifl, J. Fehr, S. L. Brunton, and J. N. Kutz.  Multi-hierarchical surrogate learning for structural dynamics of automotive crashworthiness using graph convolutional neural networks.  arXiv preprint arXiv:2402.09234, 2024.
  * [41] D. Kochkov, J. A. Smith, A. Alieva, Q. Wang, M. P. Brenner, and S. Hoyer.  Machine learning–accelerated computational fluid dynamics.  Proceedings of the National Academy of Sciences, 118(21):e2101784118, 2021.
  * [42] N. Kovachki, Z. Li, B. Liu, K. Azizzadenesheli, K. Bhattacharya, A. Stuart, and A. Anandkumar.  Neural operator: Learning maps between function spaces with applications to pdes.  Journal of Machine Learning Research, 24(89):1–97, 2023.
  * [43] K. Lee and K. T. Carlberg.  Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders.  Journal of Computational Physics, 404:108973, 2020.
  * [44] Z. Li, D. Z. Huang, B. Liu, and A. Anandkumar.  Fourier neural operator with learned deformations for pdes on general geometries.  Journal of Machine Learning Research, 24(388):1–26, 2023.
  * [45] Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar.  Fourier neural operator for parametric partial differential equations.  arXiv preprint arXiv:2010.08895, 2020.
  * [46] Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar.  Neural operator: Graph kernel network for partial differential equations.  arXiv preprint arXiv:2003.03485, 2020.
  * [47] Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, A. Stuart, K. Bhattacharya, and A. Anandkumar.  Multipole graph neural operator for parametric partial differential equations.  Advances in Neural Information Processing Systems, 33:6755–6766, 2020.
  * [48] Z. Li, N. Kovachki, C. Choy, B. Li, J. Kossaifi, S. Otta, M. A. Nabian, M. Stadler, C. Hundt, K. Azizzadenesheli, et al.  Geometry-informed neural operator for large-scale 3d pdes.  Advances in Neural Information Processing Systems, 36, 2024.
  * [49] B. List, L.-W. Chen, K. Bali, and N. Thuerey.  How temporal unrolling supports neural physics simulators.  arXiv preprint arXiv:2402.12971, 2024.
  * [50] B. Liu, N. Kovachki, Z. Li, K. Azizzadenesheli, A. Anandkumar, A. M. Stuart, and K. Bhattacharya.  A learning-based multiscale method and its application to inelastic impact problems.  Journal of the Mechanics and Physics of Solids, 158:104668, 2022.
  * [51] Q. Liu, W. Zhu, F. Ma, X. Jia, Y. Gao, and J. Wen.  Graph attention network-based fluid simulation model.  AIP Advances, 12(9), 2022.
  * [52] S. Liu, W. Du, Y. Li, Z. Li, V. Bhethanabotla, N. Rampal, O. Yaghi, C. Borgs, A. Anandkumar, H. Guo, et al.  A multi-grained symmetric differential equation model for learning protein-ligand binding dynamics.  arXiv preprint arXiv:2401.15122, 2024.
  * [53] S. Liu, Y. Li, Z. Li, Z. Zheng, C. Duan, Z.-M. Ma, O. Yaghi, A. Anandkumar, C. Borgs, J. Chayes, et al.  Symmetry-informed geometric representation for molecules, proteins, and crystalline materials.  Advances in Neural Information Processing Systems, 36, 2024.
  * [54] Y. Liu, A. Sholokhov, H. Mansour, and S. Nabi.  Physics-informed koopman network for time-series prediction of dynamical systems.  In ICLR 2024 Workshop on AI4DifferentialEquations In Science, 2024.
  * [55] M. Liu-Schiaffini, J. Berner, B. Bonev, T. Kurth, K. Azizzadenesheli, and A. Anandkumar.  Neural operators with localized integral and differential kernels.  arXiv preprint arXiv:2402.16845, 2024.
  * [56] L. Lu, P. Jin, G. Pang, Z. Zhang, and G. E. Karniadakis.  Learning nonlinear operators via deeponet based on the universal approximation theorem of operators.  Nature machine intelligence, 3(3):218–229, 2021.
  * [57] D. J. Lucia, P. S. Beran, and W. A. Silva.  Reduced-order modeling: new approaches for computational physics.  Progress in aerospace sciences, 40(1-2):51–117, 2004.
  * [58] P. Ma, P. Y. Chen, B. Deng, J. B. Tenenbaum, T. Du, C. Gan, and W. Matusik.  Learning neural constitutive laws from motion observations for generalizable pde dynamics.  In International Conference on Machine Learning, pages 23279–23300. PMLR, 2023.
  * [59] Q. Ma, H. Wu, L. Xing, J. Wang, and M. Long.  Eulagnet: Eulerian fluid prediction with lagrangian dynamics.  arXiv preprint arXiv:2402.02425, 2024.
  * [60] H. A. Majid and F. Tudisco.  Mixture of neural operators: Incorporating historical information for longer rollouts.  In ICLR 2024 Workshop on AI4DifferentialEquations In Science, 2024.
  * [61] Z. Mao, A. D. Jagtap, and G. E. Karniadakis.  Physics-informed neural networks for high-speed flows.  Computer Methods in Applied Mechanics and Engineering, 360:112789, 2020.
  * [62] R. Maulik, B. Lusch, and P. Balaprakash.  Reduced-order modeling of advection-dominated systems with recurrent neural networks and convolutional autoencoders.  Physics of Fluids, 33(3):037106, 03 2021.
  * [63] J. J. Monaghan.  Smoothed particle hydrodynamics.  In: Annual review of astronomy and astrophysics. Vol. 30 (A93-25826 09-90), p. 543-574., 30:543–574, 1992.
  * [64] H. C. Nam, J. Berner, and A. Anandkumar.  Solving Poisson equations using neural walk-on-spheres.  arXiv preprint arXiv:2406.03494, 2024.
  * [65] R. Ni and A. H. Qureshi.  Ntfields: Neural time fields for physics-informed robot motion planning.  arXiv preprint arXiv:2210.00120, 2022.
  * [66] S. Pan, S. L. Brunton, and J. N. Kutz.  Neural implicit flow: a mesh-agnostic dimensionality reduction paradigm of spatio-temporal data.  Journal of Machine Learning Research, 24(41):1–60, 2023.
  * [67] J. Pathak, S. Subramanian, P. Harrington, S. Raja, A. Chattopadhyay, M. Mardani, T. Kurth, D. Hall, Z. Li, K. Azizzadenesheli, et al.  Fourcastnet: A global data-driven high-resolution weather model using adaptive fourier neural operators.  arXiv preprint arXiv:2202.11214, 2022.
  * [68] B. Peherstorfer.  Breaking the kolmogorov barrier with nonlinear model reduction.  Notices of the American Mathematical Society, 69(5):725–733, 2022.
  * [69] J. Peng, H. Viswanath, K. Tiwari, and A. Bera.  Graph-based decentralized task allocation for multi-robot target localization.  arXiv preprint arXiv:2309.08896, 2023.
  * [70] A. Peyvan, V. Oommen, A. D. Jagtap, and G. E. Karniadakis.  Riemannonets: Interpretable neural operators for riemann problems.  arXiv preprint arXiv:2401.08886, 2024.
  * [71] T. Pfaff, M. Fortunato, A. Sanchez-Gonzalez, and P. W. Battaglia.  Learning mesh-based simulation with graph networks.  arXiv preprint arXiv:2010.03409, 2020.
  * [72] F. Pichi, B. Moya, and J. S. Hesthaven.  A graph convolutional autoencoder approach to model order reduction for parametrized pdes.  Journal of Computational Physics, 501:112762, 2024.
  * [73] M. A. Rahman, M. A. Florez, A. Anandkumar, Z. E. Ross, and K. Azizzadenesheli.  Generative adversarial neural operators.  arXiv preprint arXiv:2205.03017, 2022.
  * [74] M. A. Rahman, R. J. George, M. Elleithy, D. Leibovici, Z. Li, B. Bonev, C. White, J. Berner, R. A. Yeh, J. Kossaifi, et al.  Pretraining codomain attention neural operators for solving multiphysics pdes.  arXiv preprint arXiv:2403.12553, 2024.
  * [75] M. A. Rahman, J. Ghosh, H. Viswanath, K. Azizzadenesheli, and A. Bera.  Pacmo: Partner dependent human motion generation in dyadic human activity using neural operators.  arXiv preprint arXiv:2211.16210, 2022.
  * [76] M. A. Rahman and R. A. Yeh.  Truly scale-equivariant deep nets with fourier layers.  Advances in Neural Information Processing Systems, 36, 2024.
  * [77] M. Raissi, P. Perdikaris, and G. E. Karniadakis.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  Journal of Computational Physics, 378:686–707, 2019.
  * [78] L. Richter and J. Berner.  Robust sde-based variational formulations for solving linear pdes via deep learning.  In International Conference on Machine Learning, pages 18649–18666. PMLR, 2022.
  * [79] 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, pages 8459–8468. PMLR, 2020.
  * [80] R. Sarkar, K. S. S. Aripirala, V. S. Jadhav, S. S. Sakhinana, and V. Runkana.  Pointsage: Mesh-independent superresolution approach to fluid flow predictions.  In ICLR 2024 Workshop on AI4DifferentialEquations In Science, 2024.
  * [81] B. Shih, A. Peyvan, Z. Zhang, and G. E. Karniadakis.  Transformers as neural operators for solutions of differential equations with finite regularity, 2024.
  * [82] K. Shukla, V. Oommen, A. Peyvan, M. Penwarden, N. Plewacki, L. Bravo, A. Ghoshal, R. M. Kirby, and G. E. Karniadakis.  Deep neural operators as accurate surrogates for shape optimization.  Engineering Applications of Artificial Intelligence, 129:107615, 2024.
  * [83] J. Sirignano and K. Spiliopoulos.  DGM: A deep learning algorithm for solving partial differential equations.  Journal of computational physics, 375:1339–1364, 2018.
  * [84] M. Stoffel, F. Bamer, and B. Markert.  Deep convolutional neural networks in structural dynamics under consideration of viscoplastic material behaviour.  Mechanics Research Communications, 108:103565, 2020.
  * [85] B. Ummenhofer, L. Prantl, N. Thuerey, and V. Koltun.  Lagrangian fluid simulation with continuous convolutions.  In International Conference on Learning Representations, 2020.
  * [86] R. Vinuesa and S. L. Brunton.  The potential of machine learning to enhance computational fluid dynamics.  arXiv preprint arXiv:2110.02085, pages 1–13, 2021.
  * [87] H. Viswanath, M. A. Rahman, R. Bhaskara, and A. Bera.  Adafnio: Adaptive fourier neural interpolation operator for video frame interpolation.  arXiv preprint arXiv:2211.10791, 2022.
  * [88] H. Viswanath, M. A. Rahman, A. Vyas, A. Shor, B. Medeiros, S. Hernandez, S. E. Prameela, and A. Bera.  Neural operator: Is data all you need to model the world? an insight into the impact of physics informed machine learning.  arXiv preprint arXiv:2301.13331, 2023.
  * [89] S. Wang, J. H. Seidman, S. Sankaran, H. Wang, G. J. Pappas, and P. Perdikaris.  Bridging operator learning and conditioned neural fields: A unifying perspective.  arXiv preprint arXiv:2405.13998, 2024.
  * [90] T. Wen, K. Lee, and Y. Choi.  Reduced-order modeling for parameterized pdes via implicit neural representations.  arXiv preprint arXiv:2311.16410, 2023.
  * [91] C. White, J. Berner, J. Kossaifi, M. Elleithy, D. Pitt, D. Leibovici, Z. Li, K. Azizzadenesheli, and A. Anandkumar.  Physics-informed neural operators with exact differentiation on arbitrary geometries.  In The Symbiosis of Deep Learning and Differential Equations III, 2023.
  * [92] A. Wikner, J. Harvey, M. Girvan, B. R. Hunt, A. Pomerance, T. Antonsen, and E. Ott.  Stabilizing machine learning prediction of dynamics: Novel noise-inspired regularization tested with reservoir computing.  Neural Networks, 170:94–110, 2024.
  * [93] H. Wu, F. Xu, Y. Duan, Z. Niu, W. Wang, G. Lu, K. Wang, Y. Liang, and Y. Wang.  Spatio-temporal fluid dynamics modeling via physical-awareness and parameter diffusion guidance.  arXiv preprint arXiv:2403.13850, 2024.
  * [94] M. Xu, J. Han, A. Lou, J. Kossaifi, A. Ramanathan, K. Azizzadenesheli, J. Leskovec, S. Ermon, and A. Anandkumar.  Equivariant graph neural operator for modeling 3d dynamics.  arXiv preprint arXiv:2401.11037, 2024.
  * [95] Y. Yin, M. Kirchmeyer, J.-Y. Franceschi, A. Rakotomamonjy, and patrick gallinari.  Continuous PDE dynamics forecasting with implicit neural representations.  In The Eleventh International Conference on Learning Representations, 2023.
  * [96] X. Zhang, L. Wang, J. Helwig, Y. Luo, C. Fu, Y. Xie, M. Liu, Y. Lin, Z. Xu, K. Yan, et al.  Artificial intelligence for science in quantum, atomistic, and continuum systems.  arXiv preprint arXiv:2307.08423, 2023.
  * [97] Q. Zhou and A. Jacobson.  Thingi10k: A dataset of 10,000 3d-printing models, 2016.

##  Appendix A Additional related works

#### Time series dynamical systems

Simulating temporal dynamics in an auto-regressive manner is a particularly
challenging task due to error accumulations during long rollout [92, 49].
There have been many works that learn temporal PDEs and CFD, including [60,
54, 80, 93, 34, 36, 59, 33]. Some works have proposed neural network-based
approaches to model 3D Lagrangian dynamics, such as [85], who propose a
convolutional neural network-based approach to model the behavior of Newtonian
fluids in 3D systems. [79] propose a more general graph-based framework, but
the network suffers from high computation time on very dense graphs and is
restricted to learning physics in the full-order setting.

##  Appendix B Operator Learning

###  B.1 Background

Here, we summarize the important ingredients of neural operators. For more
details, please refer to [45]. Operator learning is a machine learning
paradigm where a neural network is trained to map between infinite-dimensional
function spaces. Let $\mathcal{G}:\mathcal{V}\rightarrow\mathcal{A}$ be a
nonlinear map between the two function spaces $\mathcal{V}$ and $\mathcal{A}$.
A neural operator is an operator parameterized by a neural network given by

$\mathcal{G_{\theta}}\colon\mathcal{V}\rightarrow\mathcal{A},\quad\theta\in\mathbb{R}^{P},$
(15)

that approximates this function mapping in the finite-dimensional space. The
learning problem can be formulated as

$\min_{\theta\in\mathbb{R}^{P}}\,\mathbb{E}_{v\sim
D}\left[\|\mathcal{G}_{\theta}(v)-\mathcal{G}(v)\|_{\mathcal{V}}^{2}\right],$
(16)

where $\|\cdot\|_{\mathcal{V}}$ is a norm on $\mathcal{V}$ and $D$ is a
probability distribution on $\mathcal{V}$. In practice, the above optimization
is posed as an empirical risk-minimization problem, defined as

$\min_{\theta\in\mathbb{R}^{P}}\frac{1}{N}\sum_{i=1}^{N}\|\mathcal{G}_{\theta}(v^{(i)})-a^{(i)}\|^{2}_{\mathcal{V}}.$
(17)

A neural operator $\mathcal{G}_{\theta}$ learns the mapping between two
functions through a sequence of point-wise and integral operators, defined as

$\mathcal{G}_{\theta}=\mathcal{L}\circ\mathcal{J}_{1}\circ...\circ\mathcal{J}_{L}\circ\mathcal{P}$
(18)

The lifting and projection layers $\mathcal{L}$ and $\mathcal{P}$ are
learnable pointwise operators that output a function with a higher and lower-
dimensional co-domain, respectively. The intermediate layers
$\mathcal{J}_{\ell}$ perform kernel integration operations with a learnable
kernel function as in (4).

##  Appendix C Interaction network

The interaction network proposed in [5] learns a relation-centric function
$\mathit{f}$ that encodes spatial interactions between the interacting nodes
within a system as a function of their interaction attributes $\mathit{r}$.
This can be represented as

$e_{t+1}=\mathit{f}_{R}(x_{1,t},x_{2,t},r)$ (19)

A node-centered function predicts the temporal dynamics of the node as a
function of the spatial interactions as follows

$x_{1,t+1}=f_{o}(e_{t+1},x_{1,t})$ (20)

In a system of $m$ nodes, the spatial interactions are represented as a graph,
where the neighborhood is defined by a ball of radius $\mathrm{r}$. This graph
is represented as $G(O,R)$, where $O$ is the collection of objects and $R$ is
the relationships between them. The interaction between them is defined as

$\mathcal{I}(G)=f_{o}(a(G,X,f_{R}(\langle x_{i},x_{j},r_{ij}\rangle)))$ (21)

Where $a$ is an aggregation function that combines all the interactions, $X$
is the set of external effects, not part of the system, such as gravitational
acceleration, etc.

###  C.1 Input and output representation

We define $\mathbf{X}_{t_{n}}\in\mathbb{R}^{Q\times d}$ to be the pointwise
positions of $Q$ particles within a $d$-dimensional system at time $n$. A
sequence of $N$ time steps is denoted as
$\mathbf{X}_{t_{0:N}}=(\mathbf{X}_{t_{0}},\dots,\mathbf{X}_{t_{N}})$. In
particular,
$\\{\mathbf{X}^{0}_{t_{n}},\dots,\mathbf{X}^{Q}_{t_{n}}\\}\in\mathbf{X}_{t_{n}}$
are the individual particles within the system. We define velocity at time $n$
as $\mathbf{V}_{t_{n}}\in\mathbb{R}^{Q\times d}$ as
$\mathbf{X}_{t_{n}}-\mathbf{X}_{t_{n-1}}$. Similarly, acceleration at time $n$
is defined as $\mathbf{A}_{t_{n}}=\mathbf{V}_{t_{n}}-\mathbf{V}_{t_{n-1}}$ or
$\mathbf{A}_{t_{n}}=\mathbf{X}_{t_{n+1}}-2\mathbf{X}_{t_{n}}+\mathbf{X}_{t_{n-1}}$.
In all these cases, $\Delta t$ is set to one for simplicity. The model learns
the instantaneous acceleration, denoted as

$\mathbf{A}_{t_{n}}=\mathcal{G}_{\theta}(\tilde{\mathbf{V}}_{t_{n-w:n}})$ (22)

where, $w$ is the window used for past time step instances and
$\tilde{\mathbf{V}}_{t_{n-w:n}}$ denotes the velocity sequence with added
noise. The particle types (water, sand, plasticine, etc.) are represented as
embeddings of dimension 16.

###  C.2 Hyperparameters

The models were implemented using Pytorch library and trained on CUDA. The
graphs were built using Pytorch Geometric module. All models were trained on
NVIDIA RTX 3060 GPUs for 1e6 steps.

The input to the model is a state vector matrix corresponding to $w$ = 6
previous time steps of each particle, along with features that represent the
material of each particle. A radius graph is constructed for the set of
particles within the input space, such that edges are added between particles
that are within the radius $r$. The nodes of the graph are the velocity
sequences for all the particles within the sparse graph.

The graph is constructed using radius_graph defined in Pytorch Geometric. The
node features and the edge features, which include the distance from the
boundary points, are encoded into latent vectors of size 128 using 2 MLPs. The
model uses one layer of interaction operator, followed by one layer of graph
neural operator. The latents are then processed by two layers of FNO. The
decoder layers are symmetric to the encoder layers. The decoder GNO uses a
projection layer with 16 channels that lifts the output to 128 channels, which
is then projected back to the physical dimensions of the input graph (2D or
3D). All MLPs within the GNO and FNO framework use gelu activation function.

#### Training noise

To prevent noise accumulation during long rollouts, the velocity sequence is
corrupted with random walk noise during training. The noise is sampled from a
normal distribution $\mathcal{N}(0,\sigma^{2})$.

#### Normalization

All velocity sequences are standardized to zero mean and unit standard
deviation. The dataset statistics are computed during training. Global mean
and variance values from the training dataset are used to compute statistics.

#### Optimizers

Optimization is done with Adam optimizer, with an initial learning rate of
1e-4 and a batch size of 4. The learning rate was decayed exponentially from
$10^{-4}$ to $10^{-6}$ using a scheduler, with a gamma of $0.1^{1/5e6}$

##  Appendix D Pseudocode

Algorithm 1 Predicting Lagrangian dynamics with GIOROM

Input: Reduced-order velocities
$\mathbf{V}_{t_{n-w:n}}=\\{\mathbf{V}^{k}_{t_{n-w:n}}\\}^{Q}_{k=1}$, full-
order points $\mathbf{\bar{X}}=\\{\mathbf{X}^{j}\\}^{P}_{j=1}$
Output:  Full-order deformation
$\mathbf{\bar{X}}_{t_{n+1}}=\\{\mathbf{X}_{t_{n+1}}^{j}\\}^{P}_{j=1}=\\{{\boldsymbol{\phi}}_{t_{n+1}}(\mathbf{X}^{j})\\}^{P}_{j=1}$

1:$\mathbf{A}_{t_{n}}$ $\leftarrow$ $\texttt{G}_{\theta}$
($\mathbf{V}_{t_{n-w:n}}$) $\triangleright$ See Section 4

2:$\mathbf{V}_{t_{n+1}}$ $\leftarrow$ $\mathbf{V}_{t_{n}}+\Delta
t\,\mathbf{A}_{t_{n}}$ $\triangleright$ See Equation (3)

3:$\mathbf{X}_{t_{n+1}}$ $\leftarrow$ $\mathbf{X}_{t_{n}}+\Delta
t\,\mathbf{V}_{t_{n+1}}$ $\triangleright$ See Equation (2)

4:$\mathbf{\bar{X}}_{t_{n+1}}\leftarrow\texttt{NeuralField}(\mathbf{X}_{t_{n+1}},\mathbf{\bar{X}})$
$\triangleright$ See Section 4.3

Algorithm 2 Training of the neural operator

Input: Reduced-order position sequence $\mathbf{X}_{t_{n-w:n}}$, ground truth
acceleration $\mathbf{\hat{A}}_{t_{n}}$
Output:  Reduced-order acceleration $\mathbf{A}_{t_{n}}$

1:$\mathbf{V}_{t_{n-w:n}}\leftarrow(\mathbf{X}_{t_{n-w:n}}-\mathbf{X}_{t_{n-w:n-1}})/\Delta
t$

2:$\mathbf{\tilde{V}}_{t_{n-w-1:n}}\leftarrow\mathbf{V}_{t_{n-w:n}}+\mathcal{N}(0,\sigma^{2})$
$\triangleright$ As explained in Section C.2

3:edges $\leftarrow$ radius$\\_$graph($\mathbf{X}_{t_{n}}$, radius)

4:edge_feats $\leftarrow$ MLP($\mathbf{X}_{t_{n}},$ edge)

5:node_feats $\leftarrow$ MLP($\mathbf{\tilde{V}}_{t_{n-w:n}})$

6:node_feats, edge_feats $\leftarrow$ IO(node_feats, edge_feats)
$\triangleright$ As defined in equation (9)

7:lgrid $\leftarrow$ linspace ([min, max])

8:latents $\leftarrow$ GNO($\mathbf{X}_{t_{n}}$, lgrid, node_feats)
$\triangleright$ See Equation (5)

9:acc $\leftarrow$ FNO (latents) $\triangleright$ See Equation (8)

10:acc$\\_$spatial $\leftarrow$ GNO (lgrid, $\mathbf{X}_{t_{n}}$, acc)

11:$\mathbf{A}_{t_{n}}$ $\leftarrow$ IO (acc_spatial, edge_feats)

12:loss $\leftarrow$ MSE($\mathbf{A}_{t_{n}}$, $\mathbf{\hat{A}}_{t_{n}}$)

##  Appendix E Additional dataset details

We model six classes of materials - elastic, plasticine, granular, Newtonian
fluids, non-Newtonian fluids, and multi-material simulations. We model them in
both 2D and 3D.

#### Newtonian Fluids

For Newtonian fluids, we use WaterDrop dataset created by [71], which is
generated using the material point method (MPM). We prepared a second dataset
to show geometry invariance in simulating the trajectories of Newtonian fluids
using the NCLAW framework. We generated 60 trajectories of 2000 time steps
with random initial velocities and 4 different geometries - Stanford bunny,
Stanford armadillo, blub (goldfish), and spot (cow). Lastly, we trained our
model on Water-3D dataset published by [71], using smooth particle
hydrodynamics (SPH).

#### Plasticine (von Mises Yield)

Using the NCLAW simulator, we generated 60 trajectories of 2000 time steps
with random initial velocities and 4 different geometries - Stanford bunny,
Stanford armadillo, blub (goldfish), and spot (cow). These trajectories were
generated using MPM.

#### Granular material (Drucker Prager sand flows)

We trained the model on 3 datasets to simulate granular media. We generated
trajectories using NCLAW simulator and on the datasets released by [71]

#### Elasticity

To simulate elasticity, we generated simulations using meshes from Thingi10k
dataset [97]. We generated 24 trajectories for 6 geometries to train the
model.

#### Non-Newtonian fluids

To train the model on non-Newtonian fluids, we used the Goop and Goop-3D
datasets. We also generated Jelly trajectories using NCLAW simulator.

#### Multimaterial

We simulated multi-material trajectories in 2D using the dataset published by
[71].

##  Appendix F Training details

Ground truth acceleration is computed from position sequences before adding
noise to the input. This is adjusted by removing the velocity noise
accumulated at the end of the random walk. This ensures that the model
corrects the noise present in the velocity.

It is to be noted that this loss is defined as a 1-step loss function over a
pair of consecutive time steps $k$ and $k+1$, imposing a strong inductive bias
towards a Markovian system. Optimizing the model for rollout over $K$ steps
would overlook the effects of instantaneous physical states (influence of
gravity, etc.), thus resulting in greater one-step errors, which would
eventually accumulate and result in larger errors during rollout.

The model was validated by full rollouts on 5 held-out validation sets per
material simulation, with performance measured by the MSE between predicted
particle positions and ground-truth particle positions.

We test our model on multiple materials, ranging from Newtonian fluids to
elastic solids, in both 2D and 3D settings. We empirically show that our model
is at least 5x faster than other neural network-based solvers with comparable
parameter counts on the same simulation trajectory. Furthermore, we show that
this speed-up doesn’t compromise the accuracy of rollout predictions. We also
highlight the generalization capability of our model to unseen initial
trajectories, graph densities, and unknown geometries. Lastly, we show the
capability of our model to perform interactive physics simulations.

The effects of the interaction operator can be seen in Figure 6. In the
absence of the interaction operator, the model fails to converge during
training.

Figure 6: Loss curves for different neural solvers. The plot highlights the
importance of parameterization of GNO with interaction operator (IO). In the
absence of this, the model fails to converge.

###  F.1 Ensuring discretization invariance

We employ a sampling strategy to ensure the model processes different
discretizations of the point cloud at each step in the training process. At
each training step, we choose an arbitrary percentage of points from the point
cloud randomly, (1%-3%) for elastic materials, (10% - 15%) for plasticine,
(20%-25%) for water, and (35%-45%) for sand. A similar strategy is employed
during inference, such that the model infers on a previously unseen sampling
percentage.

##  Appendix G Additional results

Table 6 highlights the performance of two additional baseline models. To train
the FNO on the irregular grid, the input is processed through a set of MLP
layers, which generate the latents on which the FNO layers are applied.

Table 6: This table provides additional baseline comparisons Model | Material | Noise | 1-step Loss | Rollout-Loss
---|---|---|---|---
FNO | Plasticity | 1e-4 | 2.2e-7 | 0.46
FNO | Elasticity | 3e-4 | 1.92e-8 | 1.03
FNO | Water | 3e-4 | 1.9e-7 | 0.14
FNO | Sand | 3e-4 | 2.4e-7 | 0.31
FNO | Jelly | 3e-4 | 6.73e-8 | 1.10
FNO | Multi-material | 3e-4 | 5.89e-8 | 0.975
MeshGraphNet | Plasticity | 1e-4 | 2.86e-7 | 0.096
MeshGraphNet | Elasticity | 3e-4 | 2.42e-7 | 0.083
MeshGraphNet | Water | 3e-4 | 2.35e-7 | 0.086
MeshGraphNet | Sand | 3e-4 | 7.5e-8 | 0.045
MeshGraphNet | Jelly | 3e-4 | 8.3e-8 | 0.059
MeshGraphNet | Multi-material | 3e-4 | 4.78e-8 | 0.033
Table 7: This table showcases the performance of GIOROM on several datasets,
categorized into six families of material - Newtonian fluids, elasticity,
plasticine, granular material, multi-material, and non-Newtonian fluids.

Dataset | Connectivity | Full graph | Connectivity | Sparse graph
---|---|---|---|---
|  | 1-step loss | Rollout loss |  | 1-step loss | Rollout Loss
Newtonian Fluids
WaterDrop 2D (MPM) | 0.015 | 4.33e-9 | 0.0088 | 0.045 | 5.43e-9 | 0.0098
Water 3D (MPM) | 0.035 | 4.21e-8 | 0.011 | 0.095 | 5.65e-8 | 0.016
Water (NCLAW) | 0.060 | 7.05e-9 | 0.0074 | 0.105 | 3.1e-8 | 0.010
Plasticine
Plasticine (NCLAW) | 0.060 | 1.2e-8 | 0.00944 | 0.100 | 6.4e-8 | 0.0089
Elasticity
Elasticity (FEM) | 0.100 | 1.07e-8 | 0.00349 | 0.080 | 3.6e-8 | 0.0075
Owl (FEM) | 0.090 | 2.05e-7 | 0.0069 | 0.047 | 6.8e-8 | 0.0087
Granular Material
Sand 2D (MPM) | 0.015 | 7.88e-9 | 0.011 | 0.050 | 8e-9 | 0.0095
Sand 3D (MPM) | 0.060 | 6.2e-9 | 0.00077 | 0.095 | 1.38e-8 | 0.0018
Sand (NCLAW) | 0.060 | 7.03e-9 | 0.0063 | 0.105 | 9.3e-9 | 0.0089
Non-Newtonian Fluids
Goop 2D (MPM) | 0.015 | 9.45e-9 | 0.025 | 0.030 | 7.5e-9 | 0.0078
Goop 3D (MPM) | 0.060 | 7.9e-9 | 0.00089 | 0.095 | 1.5e-8 | 0.0019
Multi-Material 2D (MPM) | 0.015 | 5.5e-9 | 0.050 | 0.025 | 7.0e-9 | 0.0056

Figure 7: The above figures depict full-order inference by GIOROM on 2D point clouds. (a) depicts granular flow, (b) represents the trajectory of jelly-like substance under gravity. (c) shows the effects of external force on a highly elastic object. (d) depicts the interaction of granular media and Newtonian fluid. Table 8: The table denotes the various training and testing geometries. Shape Params | Armadillo | Bunny | Spot | Blub
---|---|---|---|---
Mean Curvature | 2.9e-3 | 1.1e-2 | 1.3e-2 | 5.9e-3
Dirichlet Energy | 2.3e-4 | 2.5e-3 | 2.4e-3 | 8.1e-4
Figure 8: This figure highlights the performance of GIOROM on additional 3D
trajectories for water, sand, and elasticity on previously unseen geometries.
Sand and water simulations were inferred on “spot” and elasticity on
“armadillo” Figure 9: This figure showcases the resolution agnosticism through
the ability of GIOROM to infer the physics on sparse graphs while retaining
the volume of the point cloud, which is important in order to upscale it to
its original size. The top row depicts the sparse graph, while the bottom row
depicts the full-order ground truth. The discretization is randomly chosen
through a sampling threshold between 20-40% of the volume. Due to randomness
in sampling, these sets of points are very unlikely to be seen during training
Figure 10: The above graphs highlight the difference in computation times for
three graph-based models. (a) contrasts the increase in computation time with
the graph radius while (b) shows how the computation times change with an
increase in graph density.

#### Evaluation

The evaluation metrics used to evaluate the models are particle-wise one-step
MSE and rollout MSE on the held-out evaluation sets. The rollout velocity and
positions are computed using semi-implicit Euler integration as

$\mathbf{V}_{t_{k+1}}=\frac{\Delta\mathbf{X}_{t_{k}}}{\Delta t}+\Delta
t\cdot\mathcal{G}_{\theta}(\mathbf{V}_{t_{k-C:k}})$ (23)
$\mathbf{X}_{t_{k+1}}=\mathbf{X}_{t_{k}}+\Delta t\cdot\mathbf{V}_{t_{k+1}}$
(24)

In our calculations, we assume $\Delta t$ to be 1.

Table 9: The table highlights the model’s invariance to upsampling. Upsampling Rate | Loss
---|---
2$\times$ | 0.009
4$\times$ | 0.007
10$\times$ | 0.0089
100$\times$ | 0.0094

##  Appendix H Additional ablations

#### Number of message-passing layers

We show that the key bottleneck in terms of speed is the message-passing
operation within the interaction network. This operation scales with the
number of edges as $E$ = $O(K^{2})$, where $K$ is the number of nodes.

Table 10 shows the effect of message-passing layers on the computation time
per forward pass for a rollout of 200 time steps and the resulting accuracy of
the simulation. Based on the results shown in the table, a single layer of
message-passing significantly improves the speed of convergence of the
architecture, as opposed to using only integral transform layers of the graph
neural operator, which fails to converge in all of the settings.

Table 10: This table shows that the number of message-passing layers results
in a negligible improvement in rollout Loss.

Number of Message-Passing layers | Connectivity | Graph Size | Inference Time per step | Loss
---|---|---|---|---
2 | 0.077 | 2247 | 3.6 | 0.0087
4 | 0.077 | 2247 | 3.8 | 0.0091
6 | 0.077 | 2247 | 4.2 | 0.0085
8 | 0.077 | 2247 | 4.3 | 0.0083

#### Graph radius and viscosity

The following table highlights the changes in viscosity as the graph radius
varies during inference.

Table 11: This table highlights the increase in viscosity, measured by the average minimum inter-particle distance over 200 time steps and the average maximum velocity over 200 time steps. Higher values of minimum inter-particle distance denote a more rigid graph where the particles don’t collide with each other as often, and a lower average particle velocity indicates a more constrained flow Radius | Average min inter-Particle Distance | Average Max Particle Velocity
---|---|---
0.015 | 6.06e-5 | 1.7e-2
0.020 | 6.32e-5 | 8.1e-3
0.025 | 6.4e-5 | 7.6e-3
0.030 | 6.42e-5 | 6.6e-3
0.035 | 6.51e-5 | 5.7e-3
0.040 | 6.58e-5 | 5.0e-3
0.045 | 6.71e-5 | 4.6e-3
0.050 | 7.04e-5 | 4.2e-3

#### Discretization invariance w.r.t the latent grid

We show that the model is agnostic to the resolution of the latent grid during
inference. If the latent grid is too small, the performance degrades due to
data loss. However, with larger latent grid sizes, there is no significant
improvement in performance. The model was tested on latent grid dimensions of
8, 16, 32, 64, and 128. The results are shown in Figure 11a.

#### Effects of latent grid

The latent grid allows the Fourier neural operator to learn the temporal
dynamics on a regular grid of fixed size. This allows it to learn the dynamics
of non-uniform and complex geometries. Table 12 shows the performance of the
model in the absence of the latent grid.

Table 12: The table showcases invariance to grid sizes greater than 8. At sizes less than 16, the model fails to perform as well due to data loss Grid Size | Rollout Loss
---|---
128 | 0.0072
64 | 0.0074
32 | 0.0081
16 | 0.0075
8 | 0.0110

###  H.1 Effects of noise on rollout accuracy

Temporal auto-regressive models suffer from corruption of simulations due to
noise accumulation. It is important to choose the right noise scale to ensure
that the model is robust to this noise accumulation. We experimented with
different noise standard deviations and found that values between 1e-4 and
3e-4 resulted in the most stable rollouts. Table 13 highlights the effects of
noise on one-step and rollout losses. It can be observed that lower noise
scales result in higher one-step accuracy but lower rollout accuracy because
the training and validation distributions become more similar in the absence
of noise.

Table 13: The table presents the one-step and rollout-loss for different noise standard deviations on the WaterDrop2D dataset. Similar trends can be observed on other datasets. Noise | 1-step Loss | Rollout Loss
---|---|---
9e-4 | 9.15e-8 | 0.0236
6e-4 | 6.7e-8 | 0.0104
3e-4 | 3.5e-8 | 0.0075
1e-4 | 1.7e-8 | 0.0096
9e-5 | 1.53e-8 | 0.0120
6e-5 | 1.24e-8 | 0.0110
3e-5 | 1.09e-8 | 0.0260
9e-6 | 7.86e-9 | 0.0180
6e-6 | 9.31e-9 | 0.0132
3e-6 | 6.9e-9 | 0.017
Figure 11: The above graphs denote the effect of various hyperparameters on
rollout loss. (a) denotes the effect of latent grid, (b) and (d) noise on
rollout and one-step losses respectively and (c) Number of training
trajectories on the WaterDrop2D dataset

#### Effect of training dataset size on generalizability

We performed experiments to see if the performance would improve significantly
with the addition of new trajectories in the WaterDrop2D dataset. We observed
that the rollout loss steadily decreases with the addition of new
trajectories. However, this decrease is less apparent after 200 trajectories.
The loss is much higher when the number of trajectories is less than 100.

Table 14: This table shows the trends in rollout loss with the number of training trajectories. The model generalizes fairly well when trained on 150 trajectories, after which there’s a gradual improvement in performance Training Data size (# Trajectories) | Rollout Error
---|---
50 | 0.010
150 | 0.0098
200 | 0.0095
400 | 0.0092
1000 | 0.0089

###  H.2 Design decisions with minimal impact

We performed hyper-parameter tuning on WaterDrop dataset and found that the
following parameters have the least impact on the overall model performance.

#### Time window for input velocity sequence

The window used for input velocity sequence didn’t affect the accuracy of the
output by a significant amount. We experimented with window sizes of [2, 3, 5,
6, 7]. A window size of 2 allows for the network to be a Markov process, with
the model predicting the acceleration at a time step from the corresponding
acceleration at the previous time step. This can be leveraged for interactive
manipulation of the material within the simulation.

#### Number of MLP layers in FNO, GNO, and IO

We experimented with [1, 2, 4, 8, 10] as the number of layers within the
architecture. There was no significant change in accuracy across the different
configurations.

#### Graph reduction and discretization invariance

We approximately account for the integration weights of the GNO in (5) by
computing the mean of the kernel values in each neighborhood. Note that we
would not obtain a neural operator when using a sum as a reduction method
since the values would diverge in the limit of finer discretizations.

#### Number of Fourier modes

We tested the model on three configurations of Fourier modes: (16, 16), (20,
20), and (24, 24). There was no significant change in performance. However,
fewer Fourier modes resulted in faster runtimes.

#### Hidden, lifting, and projection channels in FNO

The hidden channels for FNO were chosen as [16, 32, 64, 128, 256]. The models
had the best performance when the channels were set to 32, with no significant
improvement in increase but a slight degradation in performance when reduced.

#### GNO hidden layer size

We tried three configurations of GNO hidden layers for both the non-linear
kernel and the linear kernel, i.e, [32, 64], [512, 256], and [64, 512, 1024,
256]. In the latter cases, the models became significantly bulkier without a
noticeable change in performance.

#### GNO radius

We experimented with the following radii - [0.0004, 0.0015, 0.015, 0.045,
0.100]. In each of these cases, there wasn’t a noticeable change in the
performance or the inference times.

#### Latent bounds for projected 3D grids

To optimize the model parameter count, we projected the 3D grids to 2D space
and unfolded the 3D manifolds into 2D for the FNO. We used tanh to bound the
2D space between [-1, 1]. We experimented other bounds, such as [grid.min(),
grid.max()] and [0, 1] (using sigmoid). All three of the choices had similar
performance metrics.

#### Projection layers

We experimented with three configurations for projection layers, i.e., [1, 2,
5], and observed only minor variations in the performance. To optimize for the
parameter count, we chose a single projection layer in the final model.

Figure 12: The above figure depicts the effects of several hyperparameters on
the rollout-loss.