# Implicit Neural Spatial Representations for Time-dependent PDEs

Honglin Chen  Rundi Wu  Eitan Grinspun  Changxi Zheng  Peter Yichen Chen

###### Abstract

Implicit Neural Spatial Representation (INSR) has emerged as an effective
representation of spatially-dependent vector fields. This work explores
solving time-dependent PDEs with INSR. Classical PDE solvers introduce both
temporal and spatial discretizations. Common spatial discretizations include
meshes and meshless point clouds, where each degree-of-freedom corresponds to
a location in space. While these _explicit_ spatial correspondences are
intuitive to model and understand, these representations are not necessarily
optimal for accuracy, memory usage, or adaptivity. Keeping the classical
temporal discretization unchanged (e.g., explicit/implicit Euler), we explore
INSR as an alternative spatial discretization, where spatial information is
_implicitly_ stored in the neural network weights. The network weights then
evolve over time via time integration. Our approach does _not_ require any
training data generated by existing solvers because our approach is the solver
itself. We validate our approach on various PDEs with examples involving large
elastic deformations, turbulent fluids, and multi-scale phenomena. While
slower to compute than traditional representations, our approach exhibits
higher accuracy and lower memory consumption. Whereas classical solvers can
dynamically adapt their spatial representation only by resorting to complex
remeshing algorithms, our INSR approach is intrinsically adaptive. By tapping
into the rich literature of classic time integrators, e.g., operator-splitting
schemes, our method enables challenging simulations in contact mechanics and
turbulent flows where previous neural-physics approaches struggle. Videos and
codes are available on the project page.111Project webpage:
http://www.cs.columbia.edu/cg/INSR-PDE/

Machine Learning, ICML

##  1 Introduction

Implicit neural spatial representation (INSR) [50, 78] parameterizes a
spatially-dependent vector field with a neural network. It has been proven to
be instrumental in computer graphics and vision applications, including
volumetric rendering [43], 3D reconstruction [42], signal processing [19], and
geometry processing [79]. While existing works have mostly focused on static
representations, this work aims to explore these neural representations for
dynamic simulations where the vector fields evolve over time. In particular,
we explore a fundamental class of physics simulation tasks governed by partial
differential equations (PDEs) with both _spatial_ and _temporal_ dependence,

$\displaystyle\begin{split}{\boldsymbol{\mathcal{F}}}({{\boldsymbol{f}}},{\boldsymbol{\nabla}}{{\boldsymbol{f}}},{\boldsymbol{\nabla}}^{2}{{\boldsymbol{f}}},\ldots,\dot{{{\boldsymbol{f}}}},\ddot{{{\boldsymbol{f}}}},\ldots)={\boldsymbol{0}},\\\
{{\boldsymbol{f}}}({\boldsymbol{x}},t):\Omega\times\mathcal{T}\to\mathbb{R}^{d}\
,\end{split}$ (1)

where $\Omega\in\mathbb{R}^{m}$ and $\mathcal{T}\in\mathbb{R}$ are the spatial
and temporal domains, respectively. Examples include the Navier-Stokes
equations for fluid dynamics and the elastodynamic equation for solid
mechanics.

To computationally simulate these problems, classical methods introduce both
_spatial_ and _temporal_ discretizations. On the one hand, _temporal_
discretization breaks down the entire temporal range into a finite number of
time steps $\\{t_{n}\\}_{n=0}^{{T}}$, where ${T}$ is the number of temporal
discretization samples, and $\Delta t=t_{n+1}-t_{n}$ is the time step size.
The solution to Equation 1 then becomes a list of spatially-dependent vector
fields: $\\{{{\boldsymbol{f}}}^{n}({\boldsymbol{x}})\\}_{n=0}^{{T}}$.
Classical methods then sequentially integrate from one time step ($n$) to the
next ($n+1$), using a wide range of time integrators, such as
explicit/implicit Euler [2]. On the other hand, _spatial_ discretization
represents these spatially-dependent vector fields
${{\boldsymbol{f}}}^{n}({\boldsymbol{x}})$ using grids, meshes, or point
clouds (meshless particles). For example, the grid-based linear finite element
method (FEM) [25] defines a shape function ${{N}^{i}}$ on each grid node and
represents the spatially dependent vector field as
${{\boldsymbol{f}}}^{n}({\boldsymbol{x}})=\sum_{i=1}^{{P}}{{\boldsymbol{f}}}^{n}_{i}{{N}^{i}}$,
where ${P}$ is the number of spatial samples.

While widely adopted in scientific computing applications, these classic
_spatial_ discretizations are not without drawbacks:

  1. 1.

Spatial discretization errors abound in fluid simulations as artificial
numerical diffusion [32], dissipation [21], and viscosity [56]. These errors
also appear in solid simulations as inaccurate collision resolution [45] and
numerical fractures [57].

  2. 2.

Memory usage spikes with the number of spatial samples ${P}$ [47].

  3. 3.

Adaptive meshing [48] and data structures [61] can reduce memory footprints
but are often computationally expensive and challenging to implement.

In this work, we ask: _what are the (dis)advantages of replacing a classical
numerical method’s spatial discretization with INSR while keeping intact the
time integrator?_ Unlike traditional representations that _explicitly_
discretize the spatial vector via spatial primitives (e.g., points), INSRs
_implicitly_ encode the field through neural network weights. In other words,
the field is parameterized by a neural network (typically multilayer
perceptrons), i.e.,
${{\boldsymbol{f}}}^{n}({\boldsymbol{x}})={{\boldsymbol{f}}_{\theta^{n}}}({\boldsymbol{x}})$
with $\theta^{n}$ being the network weights. As such, the memory usage for
storing the spatial field is independent of the number of spatial samples, but
rather it is determined by the number of neural network weights. We show that
under the same memory constraint, INSRs indeed achieve higher accuracy than
traditional discrete representations. Furthermore, INSRs are adaptive by
construction [78], allocating the network weights to resolve field details at
_any_ spatial location without changing the network architecture.

We emphasize that our contribution is orthogonal to the choice of _temporal_
discretization. Notably, our INSR approach works with various classic time
integrators, including explicit/implicit/midpoint Euler, variational time
integrators [28], and even operator splitting schemes [15]. Indeed, our focus
contrasts with other neural approaches, e.g., physics-informed neural network
(PINN) [54]. Whereas prior work has focused on overcoming low data
availability, efficiently solving inverse problems, or addressing high-
dimensional problems, our unique focus is on exploring incorporating various
existing classical time integrators. This ability to leverage classic time
integrators is particularly effective in highly nonlinear problems, e.g.,
turbulence, where previous neural-PDE approaches struggle, e.g., PINN.

In summary, we make the following contributions:

  * •

We present INSR as an alternative spatial representation for various time-
dependent physics simulation problems.

  * •

Compared to the classic grid, mesh, point cloud (meshless) spatial
representations, our INSR approach trades wall-clock runtime in favor of three
benefits: lower discretization error, lower memory usage, and built-in
adaptivity.

  * •

Utilizing various classic time integrators, including variational time
integrators and operator splitting schemes, INSR-based simulations capture
challenging cases in contact mechanics and turbulent flows where previous
neural-physics approaches fail.

Figure 1: Evolve neural field over time. We represent the field of interest
using a neural network ${{\boldsymbol{f}}_{\theta^{n}}}$, whose weights
$\theta^{n}$ are updated at each time step via optimization-based time
integration (Equation 3). In this case, the spatial domain $\Omega$ is the
volume encompassed by the undeformed object, and the represented field
${{\boldsymbol{f}}}$ is the deformation map. The governing PDE is the
elastodynamic equation.

##  2 Related Works

Implicit Neural Spatial Representation (INSR) uses neural networks to
parameterize spatially-dependent functions [13, 42, 43, 20, 10, 49, 11], e.g.,
a signed-distance-field, where the input is an arbitrary spatial location and
the output is its distance to the surface [50]. The nonlinear neural network’s
enormous expressivity makes INSR more accurate than its classic mesh-based and
meshless counterparts under the same memory constraint. For example, with the
same number of neural network weights as the number of mesh vertices (or
meshless particles), INSR-SDF captures more geometric details than a triangle
mesh [69]. Indeed, memory consumption of traditional representations scales
poorly with spatial resolutions. Adaptive discretizations can reduce memory,
but their generations are expensive. By contrast, neural representations are
adaptive by construction and can use their representation capacities at
arbitrary locations of interest without memory increases or data structure
alterations [78]. Because of the above-mentioned advantages, researchers have
used INSR for many other applications such as image processing [12, 62], 3D
reconstruction [72, 80], generative modeling [59, 76, 9] and geometry
processing [79, 63]. Besides, time-dependent implicit representations have
also been explored for capturing scene dynamics [51, 52].

Regarding PDE applications, INSR has been successively deployed in strictly
spatially dependent PDEs, including elastostatics [81], elliptic PDEs [14],
and geometry processing [79]. For time-dependent PDEs, Mehta et al. [41]
propose a framework for evolving INSR weights over time. However, their
approach specializes in level sets. Du & Zaki [18], Bruna et al. [6] also
evolve INSR’s network weights over time, with the goal of removing PINN’s
limited time range as well as solving high-dimensional problems where meshing
is impossible. By contrast, our work focuses on low-dimensional problems that
heavily rely on classical FEM methods and explore the INSR solver as a more
accurate, memory-efficient, and adaptive alternative.

#### Machine Learning (ML) for PDEs

One line of ML-PDE works train on data from classical solvers (or real
experiments), e.g., graph neural network (GNN) approaches [58, 53], neural
operator approaches [35, 36], and DeepONet [37]. After training, these methods
solve a new problem faster than the solver on which it was trained. However,
these methods typically do not generalize to arbitrary initial/boundary
conditions, material parameters, or geometries that are drastically different
from the ones presented in the training data [73]. Another line of ML-PDE
works does not require training data from classical solvers at all since these
methods are _solvers_ themselves, i.e., given the PDE and initial/boundary
conditions, these methods can directly solve the PDE just like the classical
solvers. These methods usually employ a physics-informed loss term (e.g.,
$\|{\boldsymbol{\mathcal{F}}}\|^{2}$) that incorporates the governing-PDE [54,
71, 74]. Our work also belongs to this “solver-type”.

Previous physics-informed approaches generally treat the temporal domain as a
continuous variable, and the PDE loss term would incorporate the entire
spatiotemporal-dependent PDE (${\boldsymbol{\mathcal{F}}}$). Krishnapriyan et
al. [31] break down the temporal domain into several subdomains and obtain
better long-term temporal integration. However, the time variable is still
treated as a continuous variable within each subdomain. By contrast, our
approach discretizes the temporal domain just like the classical solvers.
While this explicit temporal discretization can address the long-term
prediction limitation of PINN, it is not the primary goal of this work.
Instead, we show that by tapping into the rich literature of classical time
integration schemes, we can model challenging problems in contact mechanics
and turbulent flows where previous neural-PDE approaches struggle. Raissi et
al. [54], Wessels et al. [75] also explore a time-discrete approach, but their
methods specialize in Runge-Kutta schemes, while our general formulation
supports a wide range of classical time integrators, including operator-
splitting schemes.

##  3 Method: Time Integration on Neural Spatial Representations

Our goal is to solve time-dependent PDEs on neural-network-based spatial
representations. In Section 3.1, we first discuss representing spatial vector
fields with neural networks. Afterward, we will describe how to time-step the
network weights with classic time integrators.

###  3.1 Neural Networks as Spatial Representations

We parameterize each time-discretized spatial vector field with a neural
network: ${{\boldsymbol{f}}}^{n}={{\boldsymbol{f}}_{\theta^{n}}}$, where
$\theta^{n}$ are the neural network weights at time $t_{n}$. Specifically, the
field quantity at an arbitrary spatial location ${\boldsymbol{x}}\in\Omega$
can be queried via network inference
${{\boldsymbol{f}}_{\theta^{n}}}({\boldsymbol{x}})$.

Traditional representations _explicitly_ discretize the spatial vector field
using primitives such as meshes. These primitives _explicitly_ correspond to
spatial locations due to their compactly supported basis functions [25]. By
contrast, INSRs _implicitly_ encode the vector field via neural network
weights, and each weight affects the vector field globally. Such global
support is also an attribute of spectral methods [7, 8]. Compared to spectral
methods, our approach does not need to know the required complexity ahead of
time in order to determine the ideal basis functions [78]. INSR automatically
optimizes its parameters to where field detail is present.

Whereas memory consumption of _explicit_ representations scales poorly with
the number of spatial samples, memory consumption for INSR is independent of
the number of spatial samples. Instead, memory usage (for storing the vector
field) is determined by the number of network weights.

#### Network Architecture

Following the INSR literature, we adopt SIREN [64] (MLP with sinusoidal
activation) as our network architecture for its accuracy and quick convergence
speed advantages. Each MLP has a total of $\alpha$ hidden layers, each layer
of width $\beta$. The specific choice of these hyper-parameters is described
in Section 4.

#### Spatial Gradients

Classic spatial representations compute spatial gradients via basis functions.
Higher-order gradients require higher-order basis functions. By contrast, INSR
is $C^{\infty}$ by construction. We evaluate the gradients via computation-
graph-based auto-differentiation with respect to the input (not the weights).

###  3.2 Time integration

Given previous-time spatial vector fields
$\\{{{\boldsymbol{f}}}^{n}({\boldsymbol{x}})\\}_{k=0}^{n}$, we can compute the
next time-step ($t_{n+1}$) by solving an optimization problem:

$\displaystyle\begin{split}{{\boldsymbol{f}}}^{n+1}=\operatorname*{argmin}_{{{\boldsymbol{f}}}^{n+1}}\sum_{{\boldsymbol{x}}\in\mathcal{M}\subset\Omega}{\boldsymbol{\mathcal{I}}}(&\Delta
t,\\{{{\boldsymbol{f}}}^{k}({\boldsymbol{x}})\\}_{k=0}^{n+1},\\{{\nabla}{{\boldsymbol{f}}}^{k}({\boldsymbol{x}})\\}_{k=0}^{n+1},\\\
&\\{{\nabla}^{2}{{\boldsymbol{f}}}^{k}({\boldsymbol{x}})\\}_{k=0}^{n+1},\ldots)\
.\end{split}$ (2)

For example, the classic explicit/implicit Euler methods can be formulated in
this form [29] as well as variational time integrators derived from Hamilton’s
principle [28]. Operator-splitting style integrators [15] also weave
seamlessly into this formulation by solving multiple optimization problems.
Note that the particular choice of the objective function
${\boldsymbol{\mathcal{I}}}$ depends on the PDE of interest and the time
integrator choice.

This optimization formulation applies to _any_ spatial representation. It has
been explored thoroughly for classic spatial discretizations [3, 4, 22], which
is defined over a finite number of the spatial integration samples
$\mathcal{M}:=\\{{{\boldsymbol{x}}^{j}}\in\Omega\ |\ 1\leq
j\leq|\mathcal{M}|\\}$, e.g., grids.

Applying this formulation to a neural spatial representation, we optimize for

$\displaystyle\begin{split}\theta^{n+1}=\operatorname*{argmin}_{\theta^{n+1}}\sum_{{\boldsymbol{x}}\in\mathcal{M}\subset\Omega}{\boldsymbol{\mathcal{I}}}(&\Delta
t,\\{{{\boldsymbol{f}}_{\theta^{k}}}({\boldsymbol{x}})\\}_{k=0}^{n+1},\\{{\nabla}{{\boldsymbol{f}}_{\theta^{k}}}({\boldsymbol{x}})\\}_{k=0}^{n+1},\\\
&\\{{\nabla}^{2}{{\boldsymbol{f}}_{\theta^{k}}}({\boldsymbol{x}})\\}_{k=0}^{n+1},\ldots)\
\end{split}$ (3)

where $\\{\theta^{k}\\}_{k=0}^{n}$ are the (fixed, not variable) neural
network weights from previous time steps. Figure 1 illustrates our time
integration process, and Algorithm 1 provides the corresponding pseudocode. In
all the examples presented in this work, we solve this time-integration
optimization problem via Adam [30], a first-order stochastic gradient descent
method.

#### Spatial Sampling

_Explicit_ spatial representations (e.g., tetrahedra mesh) are often tied to a
particular spatial sampling; remeshing is sometimes possible but can also have
drawbacks, especially in higher dimensions [1, 48]. By contrast, _implicit_
spatial representations allow for arbitrary spatial sampling by construction
(Equation 3). Following Sitzmann et al. [64], we dynamically sample
$\mathcal{M}$ during optimization. For every gradient descent iteration in
every time step, we use a stochastic sample set $\mathcal{M}$ from the spatial
domain $\Omega$; $\mathcal{M}$ corresponds to the “mini-batch” in stochastic
gradient descent, with batch size $|\mathcal{M}|$.

#### Boundary Condition

PDEs are typically accompanied by spatial (e.g., Dirichlet or Neumann)
boundary conditions, which we formulate as additional penalty terms in the
objective Equation 3,

$\displaystyle\begin{split}\theta^{n+1}=\operatorname*{argmin}_{\theta^{n+1}}&\sum_{{\boldsymbol{x}}\in\mathcal{M}\subset\Omega}{\boldsymbol{\mathcal{I}}}(\Delta
t,\\{{{\boldsymbol{f}}_{\theta^{k}}}({\boldsymbol{x}})\\}_{k=0}^{n+1},\\{{\nabla}{{\boldsymbol{f}}_{\theta^{k}}}({\boldsymbol{x}})\\}_{k=0}^{n+1},\\\
\ldots)+\lambda&\sum_{{\boldsymbol{x}}^{b}\in\mathcal{M}^{b}\subset\partial\Omega}{\boldsymbol{\mathcal{C}}}({{\boldsymbol{f}}_{\theta^{n+1}}}({\boldsymbol{x}}^{b}),{\nabla}{{\boldsymbol{f}}_{\theta^{n+1}}}({\boldsymbol{x}}^{b}),\ldots)\
,\end{split}$ (4)

where $\lambda$ is the weighting factor and $\partial\Omega$ is the boundary
of the spatial domain. The particular choice of the boundary constraint
function ${\boldsymbol{\mathcal{C}}}$ depends on the problem of interest.

Algorithm 1 Time integration

0: $\text{initial network weights }\theta^{0},\text{timestep size }\Delta
t,$$\text{number of timesteps }N,\text{time integrator
}{\boldsymbol{\mathcal{I}}},\text{spatial domain }\Omega$

1: $n\leftarrow 0$

2: while $n<N$ do

3: $\theta^{n+1}\leftarrow\theta^{n}$

4: while not converged do

5: randomly sample $\mathcal{M}\subset\Omega$

6:
$L_{\theta^{n+1}}=\sum\limits_{{\boldsymbol{x}}\in\mathcal{M}}{\boldsymbol{\mathcal{I}}}(\Delta
t,\\{{{\boldsymbol{f}}_{\theta^{k}}}({\boldsymbol{x}})\\}_{k=0}^{n+1},\\{{\nabla}{{\boldsymbol{f}}_{\theta^{k}}}({\boldsymbol{x}})\\}_{k=0}^{n+1},\ldots)$

7: $\theta^{n+1}\leftarrow\theta^{n+1}-\alpha\nabla L_{\theta^{n+1}}$

8: $n\leftarrow n+1$

9: end while

10: end while

#### Initial Condition

The neural network is initialized using the given initial condition, i.e., the
field value at time $t=0$, by optimizing

$\displaystyle\theta^{0}=\operatorname*{argmin}_{\theta^{0}}\sum_{{\boldsymbol{x}}\in\mathcal{M}\subset\Omega}||{{\boldsymbol{f}}}_{\theta^{0}}({\boldsymbol{x}})-\hat{{{\boldsymbol{f}}}^{0}}({\boldsymbol{x}})||_{2}^{2}\
,$ (5)

where $\hat{{{\boldsymbol{f}}}^{0}}$ is the given initial condition. Similar
to Equation 3, we solve this optimization problem using Adam [30] and
stochastically sample $\mathcal{M}$ at each gradient descent iteration.

##  4 Experiments

In this section, we evaluate our method on three classic time-dependent PDEs:
the advection equation, the incompressible Euler equations, and the
elastodynamic equation.

Figure 2: 1D advection example: A Gaussian-shaped wave initially centered at $x=-1.5$ moves rightward with a constant velocity of $0.25$. From left to right, we show the mean absolute error plot over time and solutions at $t=3s$ and $t=12s$. Error is computed using $500$ uniform spatial points. Using an energy-preserving midpoint time integrator, our solution (blue) well approximates the ground truth (grey) over time, while the grid-based finite difference method (green) tends to diffuse over time.  Table 1: Quantitative results for the 1D advection example (Figure 2). Error: mean absolute error over a total of $240$ time steps, compared to the ground truth analytical solution. Error is evaluated over $500$ uniform spatial samples. Time: runtime for a total of $240$ time steps. Memory: memory usage for storing the spatial representations. Methods | Error | Time | Memory
---|---|---|---
Ours | 0.0030 | 5.33h | 3.520KB
Grid (same memory) | 0.0146 | 1.13s | 3.520KB
Grid (same error) | 0.0029 | 1.80s | 27.35KB

#### Baselines

From classical solvers, we compare with three baselines: (1) the grid-based
finite difference method [21], (2) the tetrahedral-mesh-based finite element
method [25], (3) the meshless-particle-based material point method [26]. From
neural-network-based, physics-informed approaches, we compare with another
three baselines: (4) the original PINN [54], (5) PINN with temporal sub-
domains [31], (6) physics-informed DeepONet [74, 73]. We focus on comparing
these physics-informed neural approaches because, like our method, they do not
require any training data from the classic solvers.

#### Comparison

To strike an apple-to-apple comparison, we use the same time integrator and
contact model for our approach and all the classical baselines. The only
difference is the spatial representation. Notably, one can also use more
advanced time integrators and contact models than the ones used in this work.
Nevertheless, since the baselines and our approach adopt the same time
integrators, our advantages on spatial discretization remain. For both classic
and neural baselines, we ensure that they use the same amount of memory for
storing the spatial representation as our method, e.g., grid size and the
number of network layers.

We refer readers to Appendices A and B for other implementation details (e.g.,
initial / boundary conditions, baseline setups) and additional results. The
temporal evolutions of the PDEs are best illustrated by the supplementary
video.

###  4.1 Advection Equation

Consider the classic 1D advection equation,

$\displaystyle\dfrac{\partial u}{\partial t}+(a\cdot{\nabla})u=0\ ,$ (6)

where $a$ is the advection velocity, and the vector field of interest is the
advected quantity ${{\boldsymbol{f}}}=u$.

#### Time Integration

We adopt the same time integration scheme in both the discrete grid
representation and ours. Choosing the energy-preserving midpoint method [44]
yields the time integration operator,

$\displaystyle{\boldsymbol{\mathcal{I}}}$
$\displaystyle=\|\frac{{u}^{n+1}({\boldsymbol{x}})-{u}^{n}({\boldsymbol{x}})}{\Delta
t}+(a\cdot{\nabla})(\frac{{u}^{n+1}({\boldsymbol{x}})+{u}^{n}({\boldsymbol{x}})}{2})\|^{2}_{2}\
.$ (7)

#### Results

Figure 2 shows an example where a Gaussian-shaped wave moves with constant
velocity $a=0.25$. Under the same memory usage for storing the spatial
representations, our approach uses $\alpha=2$ hidden layers of width
$\beta=20$, and the finite difference grid resolution is $901$.

Using midpoint time integrator, the grid-based method (grid-midpoint) diffuses
over time due to its spatial discretization, which is a well-known numerical
issue [17, 60]. On the contrary, our result (ours-midpoint) does not suffer
from numerical dissipation and agrees well with the ground truth at all
frames. We also tried the implicit Euler time integrator (see ours-implicit)
and found it inherits its property of energy dissipation. Choosing the
midpoint time integrator helps us preserve energy and obtain high-accuracy
results. In Table 1, we report the quantitative evaluation result for our 1D
advection example in Figure 2.

Figure 3: 2D Taylor-Green vortex simulation. Left: the mean squared error of the velocity field for $100$ timesteps. Error is computed using $48^{2}$ uniform spatial points. Right: velocity magnitude of solutions from the ground truth, ours, and the grid-based method (grid size $48$) at timestep $n=100$. Under the same memory usage (for storing the spatial representation), our solution has a significantly smaller error than the grid-based method.  Table 2: Quantitative results for the 2D Taylor-Green fluid example (Figure 3). Error: mean squared error of velocity field over total $100$ time steps, compared to the ground truth analytical solution. Time: runtime for a total of $100$ time steps. Memory: memory usage for storing the spatial representations. Methods | Error | Time | Memory
---|---|---|---
Ours | 3.35e-4 | 14.02h | 25.887KB
Grid (same memory) | 4.83e-3 | 2.91s | 27.00KB
Grid (same error) | 3.24e-4 | 189.4s | 12.00MB

###  4.2 Incompressible Euler Equations

Figure 4: Two vortices of different scales. We show the advected density field after $2.5$ seconds from the reference (top-left), our method (top-right), the grid-based method of resolution $48$ (bottom-left), and PINN (bottom-right). The reference is obtained using the high-resolution grid-based method (we use resolution $1024$) and serves as a good approximation of the ground truth. Our MLP ($\alpha=3,\beta=32$) has the same memory footprint as grids of resolution $48$. PINN uses the same MLP network as ours. Under the same memory constraint, our approach suffers less dissipation, captures more vorticity, and best resembles the reference solution, whose grids take $\sim 450\times$ memory compared to our network. See Figure 10 for the initial condition of this example.  Table 3: Quantitative results for the two-vortices fluid example in Figure 4. Error: average absolute error of kinetic energy over a total of $50$ timesteps, compared to the reference solution. Kinetic energy is computed using $1024^{2}$ uniform samples. Time: runtime for a total of $50$ timesteps. Memory: memory usage for storing spatial representations. Ours and Grid-based [67] use the operator splitting scheme; Ours-residual, PINN [54], PINN-sub [31] and piDeepONet [74] use the residual of the Euler equation as the objectives (see Equation 23 and Equation 24).  Methods | Error | Time | Memory
---|---|---|---
Ours | 2.24e2 | 10.81h | 25.887KB
Ours-residual | 2.97e4 | 10.07h | 25.887KB
Grid-based | 1.07e4 | 1.78s | 27.00KB
PINN | 2.25e4 | 7.88h | 26.137KB
PINN-sub | 3.21e4 | 20.83h | 26.137KB
piDeepONet | 2.11e4 | 9.42h | 65.855KB

In the incompressible Euler Equations

$\displaystyle\begin{split}\rho_{f}(\dfrac{\partial{\boldsymbol{u}}}{\partial
t}+{\boldsymbol{u}}\cdot{\nabla}{\boldsymbol{u}})&=-{\nabla}p+\rho_{f}{\boldsymbol{g}},\\\
{\nabla}\cdot{\boldsymbol{u}}&=0,\end{split}$ (8)

the vector field of interest is the fluid velocity field
${{\boldsymbol{f}}}={\boldsymbol{u}}$; $p$ is the pressure, ${\boldsymbol{g}}$
is the external force, and $\rho_{f}$ is the fluid density. In our
experiments, we consider $\rho_{f}=1$ and ${\boldsymbol{g}}=0$. The pressure
field $p$ is represented with another MLP network.

#### Time Integration

We apply the Chorin-style operator splitting scheme [15, 67] to both the
neural spatial and finite-difference grid representations. This scheme
converts the highly nonlinear PDE into three linear PDEs, which significantly
eases the challenge of solving. The entire scheme breaks down into three
sequential steps: advection (adv), pressure projection (pro), and velocity
correction (cor).

Advection uses a semi-Lagrangian method, encoded by the operator [68]

$\displaystyle{\boldsymbol{\mathcal{I}}}_{adv}$
$\displaystyle=\|{{\boldsymbol{u}}_{adv}}^{n+1}({\boldsymbol{x}})-{{\boldsymbol{u}}}^{n}({\boldsymbol{x}}_{backtrack})\|^{2}_{2}\
,$ (9)

whose optimization yields the advected velocity
${{\boldsymbol{u}}_{adv}}^{n+1}$. The backtracked location is given by
${\boldsymbol{x}}_{backtrack}={\boldsymbol{x}}-\Delta
t{{\boldsymbol{u}}}^{n}({\boldsymbol{x}})$. While traditional discrete
representations compute the backtracked velocity using interpolation (e.g.,
linear basis function), our approach _requires no interpolation_ , only direct
evaluation via network inference at ${\boldsymbol{x}}_{backtrack}$.

Pressure projection is encapsulated by the operator

$\displaystyle{\boldsymbol{\mathcal{I}}}_{pro}$
$\displaystyle=\|{\nabla}^{2}{p}^{n+1}({\boldsymbol{x}})-{\nabla}\cdot{{\boldsymbol{u}}_{adv}}^{n+1}({\boldsymbol{x}})\|^{2}_{2}.$
(10)

Plugging ${\boldsymbol{\mathcal{I}}}_{pro}$ into the optimization solver, we
obtain the pressure ${p}^{n+1}$ that enforces incompressibility. Note that the
MLP that represents the velocity field ${\boldsymbol{u}}_{adv}$ is kept fixed
in this step.

Velocity correction is formulated by the operator

$\displaystyle{\boldsymbol{\mathcal{I}}}_{cor}$
$\displaystyle=\|{{\boldsymbol{u}}}^{n+1}-({{\boldsymbol{u}}_{adv}}^{n+1}({\boldsymbol{x}})-{\nabla}{p}^{n+1}({\boldsymbol{x}}))\|^{2}_{2}\
,$ (11)

which adds the pressure gradient to the advected velocity yielding the
_incompressible_ velocity ${{\boldsymbol{u}}}^{n+1}$.

#### Results

We first test our method on the 2D Taylor-Green vortex with zero viscosity
[70, 5]. The closed-form analytical solution is given by:
${\boldsymbol{u}}({\boldsymbol{x}},t)=(\sin x\cos y,-\cos x\sin y)$ for
${\boldsymbol{x}}\in[0,2\pi]\times[0,2\pi]$. To compare under the same memory
usage (for storing the velocity field), we use $\alpha=3$ hidden layers of
width $\beta=32$ for our MLP and set the grid size to $48$ for the grid-based
projection method [67]. We set $\Delta t=0.05$ and execute both methods for
$100$ timesteps. In Figure 3, we show the mean squared error of the solved
velocity field over time. In Table 2, we report the quantitative evaluation
result for our 2D Taylor-Green example in Figure 3. This example demonstrates
that our method excellently preserves a stationary solution and obtains a much
smaller error than the grid-based method.

For discrete grid representation, efficiently capturing multi-scale details
usually requires difficult-to-implement adaptive data structures [61].
Instead, INSRs are adaptive by construction [78] and enable us to capture more
details under the same memory storage. We set up an example where the initial
velocity field is composed of two Taylor-Green vortices of different scales
(see Figure 10 for illustration).

Figure 5: Elastic tension test. We use $\alpha=3$ hidden layers of width
$\beta=68$ for our MLP, which takes the same memory as the FEM mesh (0.8K
vertices, 1.5K faces) and MPM point cloud (1.7K points).  Figure 6: Error of
the elastic tension test. We visualize the L2 distance error of our result and
the FEM result (0.8K vertices). Errors are calculated against the reference
result obtained by high-resolution FEM. Our result obtains smaller errors.

In Figure 4 and Table 3, we show our results on this example and compare with
the grid-based projection method [67], PINN [54], PINN with temporal sub-
domains (PINN-sub) [31] and piDeepONet [74]. All methods are compared under
the same memory storage for spatial representations. For PINN and PINN-sub, we
use the same MLP structure as ours. Detailed setups for these baselines can be
found in Section A.3. We execute our approach, the grid-based method, and
PINN-sub for $50$ timesteps with $\Delta t=0.05$, and train PINN and
piDeepONet with the same temporal range of $2.5$ seconds. Our approach can
capture the fine details of the smaller vortex, best approximate the reference
solution, and has the most negligible energy dissipation.

Table 4: Quantitative results for the 2D elasticity tension example in Figure 6. Error: infinity norm of L2 distance w.r.t. the high-resolution ground truth. Time: total runtime until convergence. Memory: memory usage for storing spatial representations.  Methods | Error | Time | Memory
---|---|---|---
Ours | 8.82e-2 | 38.33m | 56.32KB
Mesh-based | 1.99e-1 | 22.04s | 54.00KB

Moreover, the superior accuracy of our approach also comes from the usage of
the operator-splitting time integration scheme. While PINN, PINN-sub, and
piDeepONet also use implicit neural representations like ours, they treat time
as a continuous variable (i.e., part of the network inputs). Therefore, they
cannot use the operator-splitting scheme but only employ the residual of the
Euler equations as the training objective, which is more challenging to
optimize. The time-discrete PINN proposed by Raissi et al. [54] specializes in
Runge-Kutta schemes and does not support operator splitting in its current
form either. We verify these observations by changing the objectives of our
method to a similar training objective used by them, i.e., the residual of the
Euler equation (Equation 23) with a fully-implicit-time-discretization
(Equation 24). After such change, we obtain a significantly worse result (see
ours-residual in Table 3), which confirms the necessity of using the operator-
splitting scheme. This experiment demonstrates that just replacing
PINN/SIREN’s time-dependent formulation is insufficient, but the particular
choice of temporal discretization matters.

Figure 7: An elastic square collides with a circle. Under the same memory
footprint (for storing the spatial representation), the result of mesh-based
FEM (top) conforms poorly at the collision interface. In contrast, our result
(bottom) fits the boundary more gracefully.

###  4.3 Elastodynamic Equation

Lastly, we study the Elastodynamic equation

$\displaystyle\rho_{0}\ddot{\boldsymbol{\phi}}={\nabla}\cdot{\boldsymbol{P}}({\boldsymbol{F}})+\rho_{0}{\boldsymbol{b}}$
(12)

that describe the motions of deformable solids [23]. The vector field of
interest is the deformation map ${{\boldsymbol{f}}}=\boldsymbol{\phi}$. Here
$\rho_{0}$ is the density in the reference space, ${\boldsymbol{P}}$ is the
first Piola-Kirchhoff stress, ${\boldsymbol{F}}={\nabla}\boldsymbol{\phi}$ is
the deformation gradient, $\dot{\boldsymbol{\phi}}$ and
$\ddot{\boldsymbol{\phi}}$ are the velocity and acceleration, and
${\boldsymbol{b}}$ is the body force.

We assume a hyper-elasticity constitutive law, i.e.,
${\boldsymbol{P}}=\tfrac{\partial\Psi}{\partial{\boldsymbol{F}}}$, where
$\Psi$ is the energy density function. In particular, we assume a variant of
the stable Neo-Hookean energy [65]

$\displaystyle\Psi=\frac{\lambda}{2}\operatorname{\sf
tr}^{2}(\Sigma-{\boldsymbol{I}})+\mu(\operatorname{det}({\boldsymbol{F}})-1)^{2},$
(13)

where $\lambda$ and $\mu$ are the first and second lame parameters, $\Sigma$
are the singular values of the deformation gradient ${\boldsymbol{F}}$, and
$\operatorname{det}({\boldsymbol{F}})$ is the determinant of the deformation
gradient ${\boldsymbol{F}}$. When $\mu=0$, the elastic energy recovers the As-
Rigid-As-Possible energy [66].

#### Time Integration

We apply the variational time integration scheme [22, 27] to the (1)
tetrahedral finite element method, (2) the material point method, and (3) our
neural representation,

$\displaystyle\begin{split}{\boldsymbol{\mathcal{I}}}=&\underbrace{\frac{1}{2}\rho_{0}({\dot{\boldsymbol{\phi}}}^{n+1}-{\dot{\boldsymbol{\phi}}}^{n})^{T}({\dot{\boldsymbol{\phi}}}^{n+1}-{\dot{\boldsymbol{\phi}}}^{n})}_{\text{kinetic
energy}}+\\\
&\underbrace{\vphantom{\frac{1}{2}}\Psi({\boldsymbol{\phi}}^{n+1})}_{\begin{subarray}{c}\text{elastic}\\\
\text{energy}\end{subarray}}\underbrace{\vphantom{\frac{1}{2}}-\rho_{0}{\boldsymbol{b}}^{T}{\boldsymbol{\phi}}^{n+1}}_{\begin{subarray}{c}\text{external
force}\\\ \text{potential}\end{subarray}}\ ,\end{split}$ (14)

where
${\dot{\boldsymbol{\phi}}}^{n+1}=({\boldsymbol{\phi}}^{n+1}-{\boldsymbol{\phi}}^{n})/\Delta
t$, $\rho_{0}$ is the density, ${\boldsymbol{b}}$ is the external force. We
can also incorporate boundary conditions, e.g., positional and contact
constraints, by introducing additional energy terms [4, 34] (see Section A.4).
These energy terms allow us to simulate challenging contact problems where the
material impacts a collision surface at high speed.

#### Results

Table 5: Quantitative results for the collision example in Figure 7. Error:
maximum overlapping distance between the square and the circle. Time: runtime
for a total of $10$ timesteps. Memory: memory usage for storing spatial
representations.

| Error | Time | Memory
---|---|---|---
Methods | 0.2s | 0.6s | 0.8s
Ours | 1.62e-2 | 1.04e-2 | 1.06e-2 | 23.0m | 56.32KB
Mesh-based | 3.45e-2 | 5.47e-2 | 4.23e-2 | 98.2s | 54.00KB

We first evaluate our method on a typical 2D example for elastic tension test,
and compare with the traditional finite element method (FEM, [25, 55]) using
tetrahedral mesh representation. We use $\alpha=3$ hidden layers of width
$\beta=68$ for our MLP, which takes the same memory as the meshes used by FEM
(0.8K vertices, 1.5K faces). The geometry is rendered as point clouds for our
method in all our examples. As shown in Figure 6 and Table 4, our method
obtains a more minor error than FEM. Classic mesh-less particle technique
(material point method, MPM) suffers incorrect numerical fracture in this
example (see Figure 5).

Figure 8: A bunny collides with the ground in 3D.  Using INSR, our method
(green, left) captures more intricate geometry details and complex dynamics
compared to the traditional mesh-based FEM (blue, right) under the same memory
usage.

Figure 9:  PINN’s result on the square-circle collision example in Figure 7.
Our approach successfully captures the correct contact behaviors, but PINN
fails to do so under the same memory usage and network structure. To
demonstrate this, we present PINN’s failure results trained with both 2K and
10K ADAM iterations (top). Despite making progress in minimizing the loss
(bottom) under various hyper-parameters, such as learning rate and collision
ratio, PINN was unable to properly capture the highly discontinuous contact
elasticity. Moreover, even with additional training iterations (i.e., more
than 10K), the training loss did not improve.

By using INSR, our method can capture more intricate details than the
traditional discrete representations under the same memory usage. In Figure 7,
we show that our method allows the deformed square to gracefully fit the
boundary of the sphere during the non-trivial collision. In contrast, the
mesh-based FEM struggles to produce smooth results due to its insufficient
mesh resolution. To alleviate such artifacts, the mesh-based FEM either needs
to increase resolutions, thus inducing higher memory cost, or conducts complex
remeshing [48]. As shown in Figure 8 and Figure 13, our method allows for more
complex dynamics and fine geometry details compared to the mesh-based FEM
under the same memory footprint.

Note that we adopt the same collision detection and handling strategy for both
the neural representation and the mesh-based representation (FEM).
Specifically, we use a spring-like penalty force and the corresponding energy
to move the collided point out of its collision surface, similar to [40, 77].
Since our approach and FEM share the same time integration scheme and the same
collision handling method, the difference strictly stems from the underlying
spatial representations.

Finally, in Figure 9, we compare our time-independent formulation to PINN’s
time-dependent approach on non-trivial collision cases. PINN struggles to
capture the correct contact behavior when extreme nonlinearities and
discontinuities involve. The most likely reason is that PINN models the time
continuously, but collision is highly discontinuous in time. Therefore, the
loss function of PINN is prone to abrupt change when the collision happens,
and collision forces are involved. Our approach is able to adopt the
variational time integrator and use the incremental potential as our loss
function, which is known for its stability.

##  5 Discussion and Conclusion

This work explores INSR for numerically modeling time-dependent PDEs. Combined
with a wide range of classical time integrators, the INSR solver captures
various advection, elasticity, and fluid phenomena and outperforms previous
physics-informed neural network approaches on multiscale, turbulent flows and
contact mechanics problems. Compared to classic _explicit_ representations
(e.g., grids and meshes), our approach offers improved accuracy, reduced
memory, and automatic adaptivity.

While offering important benefits, INSR-based PDE time-stepping requires
longer wall-clock computation time than existing methods. (For reference, PINN
also takes longer to solve forward problems than classic FEMs. See also Table
1 by Zehnder et al. [81] and Section 7 by Yang et al. [79].) Optimizing
_globally-supported_ neural network weights takes longer than optimizing
_locally-supported_ grid values, even if there are fewer neural network
weights than the number of grid nodes. For instance, for the bunny example
(Figure 8), our neural network optimization takes around 30 minutes per
timestep while the corresponding FEM simulation takes less than 1 minute.

Facing this wall clock vs. memory/accuracy/adaptivity trade-off, we believe an
exciting future direction is hybrid mesh-neural spatial representations that
aim for the best of both worlds. Indeed, our work does not advocate INSR as
the “perfect” spatial representation for solving PDEs. Instead, we view our
work as a stepping stone toward future hybrid mesh-neural PDE solvers. For
this matter, recent hybrid representation works [46, 69, 39] have offered
promising results in reducing training time from hours to seconds while
keeping the expressiveness of INSR. We hope our work will serve as a benchmark
for future representations.

Our work demonstrates the effectiveness of INSR in solving time-dependent PDEs
and observes empirical convergence under refinement (see Figure 12). Future
work may consider a theoretical analysis of convergence and stability. More
challenging physical phenomena, such as turbulence and intricate contacts, are
also important future directions. Currently, our work enforces “soft” boundary
conditions. Enforcing “hard” boundary conditions on a neural network is
another exciting direction [38].

## Acknowledgements

This research was partially supported by National Science Foundation (1910839)
and a Natural Sciences and Engineering Research Council of Canada (NSERC)
Discovery Grant [RGPIN-2021-03733]. We thank Henrique Teles Maia for
proofreading and all the anonymous reviewers for their helpful comments and
suggestions. Lastly, we want to thank all our math professors (James Scott,
Joseph Teran, among many others) who instilled in us a love for PDEs.

## References

  * Alliez et al. [2002] Alliez, P., Meyer, M., and Desbrun, M.  Interactive geometry remeshing.  _ACM Transactions on Graphics (TOG)_ , 21(3):347–354, 2002.
  * Ascher & Petzold [1998] Ascher, U. M. and Petzold, L. R.  _Computer methods for ordinary differential equations and differential-algebraic equations_ , volume 61.  Siam, 1998.
  * Batty et al. [2007] Batty, C., Bertails, F., and Bridson, R.  A fast variational framework for accurate solid-fluid coupling.  _ACM Transactions on Graphics (TOG)_ , 26(3):100–es, 2007.
  * Bouaziz et al. [2014] Bouaziz, S., Martin, S., Liu, T., Kavan, L., and Pauly, M.  Projective dynamics: Fusing constraint projections for fast simulation.  _ACM transactions on graphics (TOG)_ , 33(4):1–11, 2014.
  * Brachet et al. [1983] Brachet, M. E., Meiron, D. I., Orszag, S. A., Nickel, B., Morf, R. H., and Frisch, U.  Small-scale structure of the taylor–green vortex.  _Journal of Fluid Mechanics_ , 130:411–452, 1983.
  * Bruna et al. [2022] Bruna, J., Peherstorfer, B., and Vanden-Eijnden, E.  Neural galerkin scheme with active learning for high-dimensional evolution equations.  _arXiv preprint arXiv:2203.01360_ , 2022.
  * Canuto et al. [2007a] Canuto, C., Hussaini, M. Y., Quarteroni, A., and Zang, T. A.  _Spectral methods: fundamentals in single domains_.  Springer Science & Business Media, 2007a.
  * Canuto et al. [2007b] Canuto, C., Hussaini, M. Y., Quarteroni, A., and Zang, T. A.  _Spectral methods: evolution to complex geometries and applications to fluid dynamics_.  Springer Science & Business Media, 2007b.
  * Chan et al. [2022] Chan, E. R., Lin, C. Z., Chan, M. A., Nagano, K., Pan, B., De Mello, S., Gallo, O., Guibas, L. J., Tremblay, J., Khamis, S., et al.  Efficient geometry-aware 3d generative adversarial networks.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pp. 16123–16133, 2022.
  * Chen et al. [2021a] Chen, P. Y., Chiaramonte, M., Grinspun, E., and Carlberg, K.  Model reduction for the material point method via an implicit neural representation of the deformation map.  _arXiv preprint arXiv:2109.12390_ , 2021a.
  * Chen et al. [2022] Chen, P. Y., Xiang, J., Cho, D. H., Pershing, G., Maia, H. T., Chiaramonte, M., Carlberg, K., and Grinspun, E.  CROM: Continuous reduced-order modeling of PDEs using implicit neural representations.  _arXiv preprint arXiv:2206.02607_ , 2022.
  * Chen et al. [2021b] Chen, Y., Liu, S., and Wang, X.  Learning continuous image representation with local implicit image function.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , pp. 8628–8638, 2021b.
  * Chen & Zhang [2019] Chen, Z. and Zhang, H.  Learning implicit fields for generative shape modeling.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pp. 5939–5948, 2019.
  * Chiaramonte et al. [2013] Chiaramonte, M., Kiener, M., et al.  Solving differential equations using neural networks.  _Machine Learning Project_ , 1, 2013.
  * Chorin [1968] Chorin, A. J.  Numerical solution of the navier-stokes equations.  _Mathematics of computation_ , 22(104):745–762, 1968.
  * Chuang & Barba [2022] Chuang, P.-Y. and Barba, L. A.  Experience report of physics-informed neural networks in fluid simulations: pitfalls and frustration.  _arXiv preprint arXiv:2205.14249_ , 2022.
  * Courant et al. [1952] Courant, R., Isaacson, E., and Rees, M.  On the solution of nonlinear hyperbolic differential equations by finite differences.  _Communications on pure and applied mathematics_ , 5(3):243–255, 1952.
  * Du & Zaki [2021] Du, Y. and Zaki, T. A.  Evolutional deep neural network.  _Physical Review E_ , 104(4):045303, 2021.
  * Du et al. [2021] Du, Y., Collins, K., Tenenbaum, J., and Sitzmann, V.  Learning signal-agnostic manifolds of neural fields.  _Advances in Neural Information Processing Systems_ , 34:8320–8331, 2021.
  * Dupont et al. [2022] Dupont, E., Kim, H., Eslami, S., Rezende, D., and Rosenbaum, D.  From data to functa: Your data point is a function and you should treat it like one.  _arXiv preprint arXiv:2201.12204_ , 2022.
  * Fedkiw et al. [2001] Fedkiw, R., Stam, J., and Jensen, H. W.  Visual simulation of smoke.  In _Proceedings of the 28th annual conference on Computer graphics and interactive techniques_ , pp. 15–22, 2001.
  * Gast et al. [2015] Gast, T. F., Schroeder, C., Stomakhin, A., Jiang, C., and Teran, J. M.  Optimization integrator for large time steps.  _IEEE transactions on visualization and computer graphics_ , 21(10):1103–1115, 2015.
  * Gonzalez & Stuart [2008] Gonzalez, O. and Stuart, A. M.  _A first course in continuum mechanics_ , volume 42.  Cambridge University Press, 2008.
  * Hu et al. [2019] Hu, Y., Li, T.-M., Anderson, L., Ragan-Kelley, J., and Durand, F.  Taichi: a language for high-performance computation on spatially sparse data structures.  _ACM Transactions on Graphics (TOG)_ , 38(6):1–16, 2019.
  * Hughes [2012] Hughes, T. J.  _The finite element method: linear static and dynamic finite element analysis_.  Courier Corporation, 2012.
  * Jiang et al. [2016] Jiang, C., Schroeder, C., Teran, J., Stomakhin, A., and Selle, A.  The material point method for simulating continuum materials.  In _ACM SIGGRAPH 2016 Courses_ , SIGGRAPH ’16. Association for Computing Machinery, 2016.
  * Kane et al. [2000a] Kane, C., Marsden, J. E., Ortiz, M., and West, M.  Variational integrators and the newmark algorithm for conservative and dissipative mechanical systems.  _International Journal for Numerical Methods in Engineering_ , 49(10):1295–1325, 2000a.
  * Kane et al. [2000b] Kane, C., Marsden, J. E., Ortiz, M., and West, M.  Variational integrators and the newmark algorithm for conservative and dissipative mechanical systems.  _International Journal for numerical methods in engineering_ , 49(10):1295–1325, 2000b.
  * Kharevych et al. [2006] Kharevych, L., Wei, W., Tong, Y., Kanso, E., Marsden, J. E., Schröder, P., and Desbrun, M.  _Geometric, variational integrators for computer animation_.  Eurographics Association, 2006.
  * Kingma & Ba [2014] Kingma, D. P. and Ba, J.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Krishnapriyan et al. [2021] Krishnapriyan, A., Gholami, A., Zhe, S., Kirby, R., and Mahoney, M. W.  Characterizing possible failure modes in physics-informed neural networks.  _Advances in Neural Information Processing Systems_ , 34:26548–26560, 2021.
  * Lantz [1971] Lantz, R.  Quantitative evaluation of numerical diffusion (truncation error).  _Society of Petroleum Engineers Journal_ , 11(03):315–320, 1971.
  * Levin [2020] Levin, D. I.  Bartels: A lightweight collection of routines for physics simulation, 2020.  https://github.com/dilevin/Bartels.
  * Li et al. [2020a] Li, M., Ferguson, Z., Schneider, T., Langlois, T. R., Zorin, D., Panozzo, D., Jiang, C., and Kaufman, D. M.  Incremental potential contact: intersection-and inversion-free, large-deformation dynamics.  _ACM Trans. Graph._ , 39(4):49, 2020a.
  * Li et al. [2020b] Li, Z., Kovachki, N., Azizzadenesheli, K., Liu, B., Bhattacharya, K., Stuart, A., and Anandkumar, A.  Fourier neural operator for parametric partial differential equations.  _arXiv preprint arXiv:2010.08895_ , 2020b.
  * Li et al. [2020c] Li, Z., Kovachki, N., Azizzadenesheli, K., Liu, B., Bhattacharya, K., Stuart, A., and Anandkumar, A.  Neural operator: Graph kernel network for partial differential equations.  _arXiv preprint arXiv:2003.03485_ , 2020c.
  * Lu et al. [2019] Lu, L., Jin, P., and Karniadakis, G. E.  Deeponet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators.  _arXiv preprint arXiv:1910.03193_ , 2019.
  * Lu et al. [2021] Lu, L., Pestourie, R., Yao, W., Wang, Z., Verdugo, F., and Johnson, S. G.  Physics-informed neural networks with hard constraints for inverse design.  _SIAM Journal on Scientific Computing_ , 43(6):B1105–B1132, 2021.
  * Martel et al. [2021] Martel, J. N., Lindell, D. B., Lin, C. Z., Chan, E. R., Monteiro, M., and Wetzstein, G.  Acorn: Adaptive coordinate networks for neural scene representation.  _arXiv preprint arXiv:2105.02788_ , 2021.
  * McAdams et al. [2011] McAdams, A., Zhu, Y., Selle, A., Empey, M., Tamstorf, R., Teran, J., and Sifakis, E.  Efficient elasticity for character skinning with contact and collisions.  In _ACM SIGGRAPH 2011 papers_ , pp. 1–12. 2011.
  * Mehta et al. [2022] Mehta, I., Chandraker, M., and Ramamoorthi, R.  A level set theory for neural implicit evolution under explicit flows.  In _Computer Vision–ECCV 2022: 17th European Conference, Tel Aviv, Israel, October 23–27, 2022, Proceedings, Part II_ , pp. 711–729. Springer, 2022.
  * Mescheder et al. [2019] Mescheder, L., Oechsle, M., Niemeyer, M., Nowozin, S., and Geiger, A.  Occupancy networks: Learning 3d reconstruction in function space.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pp. 4460–4470, 2019.
  * Mildenhall et al. [2020] Mildenhall, B., Srinivasan, P. P., Tancik, M., Barron, J. T., Ramamoorthi, R., and Ng, R.  Nerf: Representing scenes as neural radiance fields for view synthesis.  In _European conference on computer vision_ , pp. 405–421. Springer, 2020.
  * Mullen et al. [2009] Mullen, P., Crane, K., Pavlov, D., Tong, Y., and Desbrun, M.  Energy-preserving integrators for fluid animation.  _ACM Transactions on Graphics (TOG)_ , 28(3):1–8, 2009.
  * Müller et al. [2015] Müller, M., Chentanez, N., Kim, T.-Y., and Macklin, M.  Air meshes for robust collision handling.  _ACM Transactions on Graphics (TOG)_ , 34(4):1–9, 2015.
  * Müller et al. [2022] Müller, T., Evans, A., Schied, C., and Keller, A.  Instant neural graphics primitives with a multiresolution hash encoding.  _arXiv preprint arXiv:2201.05989_ , 2022.
  * Museth [2013] Museth, K.  Vdb: High-resolution sparse volumes with dynamic topology.  _ACM transactions on graphics (TOG)_ , 32(3):1–22, 2013.
  * Narain et al. [2012] Narain, R., Samii, A., and O’brien, J. F.  Adaptive anisotropic remeshing for cloth simulation.  _ACM transactions on graphics (TOG)_ , 31(6):1–10, 2012.
  * Pan et al. [2022] Pan, S., Brunton, S. L., and Kutz, J. N.  Neural implicit flow: a mesh-agnostic dimensionality reduction paradigm of spatio-temporal data.  _arXiv preprint arXiv:2204.03216_ , 2022.
  * Park et al. [2019] Park, J. J., Florence, P., Straub, J., Newcombe, R., and Lovegrove, S.  Deepsdf: Learning continuous signed distance functions for shape representation.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pp. 165–174, 2019.
  * Park et al. [2021a] Park, K., Sinha, U., Barron, J. T., Bouaziz, S., Goldman, D. B., Seitz, S. M., and Martin-Brualla, R.  Nerfies: Deformable neural radiance fields.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pp. 5865–5874, 2021a.
  * Park et al. [2021b] Park, K., Sinha, U., Hedman, P., Barron, J. T., Bouaziz, S., Goldman, D. B., Martin-Brualla, R., and Seitz, S. M.  Hypernerf: a higher-dimensional representation for topologically varying neural radiance fields.  _ACM Transactions on Graphics (TOG)_ , 40(6):1–12, 2021b.
  * Pfaff et al. [2020] Pfaff, T., Fortunato, M., Sanchez-Gonzalez, A., and Battaglia, P. W.  Learning mesh-based simulation with graph networks.  _arXiv preprint arXiv:2010.03409_ , 2020.
  * Raissi et al. [2019] Raissi, M., Perdikaris, P., and Karniadakis, G. E.  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\.
  * Reddy [2019] Reddy, J. N.  _Introduction to the Finite Element Method_.  McGraw-Hill Education, New York, 4th edition edition, 2019.
  * Roache [1998] Roache, P. J.  _Fundamentals of computational fluid dynamics_.  Hermosa Publishers, 1998.
  * Sadeghirad et al. [2011] Sadeghirad, A., Brannon, R. M., and Burghardt, J.  A convected particle domain interpolation technique to extend applicability of the material point method for problems involving massive deformations.  _International Journal for numerical methods in Engineering_ , 86(12):1435–1456, 2011.
  * Sanchez-Gonzalez et al. [2020] Sanchez-Gonzalez, A., Godwin, J., Pfaff, T., Ying, R., Leskovec, J., and Battaglia, P.  Learning to simulate complex physics with graph networks.  In _International Conference on Machine Learning_ , pp. 8459–8468. PMLR, 2020.
  * Schwarz et al. [2020] Schwarz, K., Liao, Y., Niemeyer, M., and Geiger, A.  Graf: Generative radiance fields for 3d-aware image synthesis.  _Advances in Neural Information Processing Systems_ , 33:20154–20166, 2020.
  * Selle et al. [2008] Selle, A., Fedkiw, R., Kim, B., Liu, Y., and Rossignac, J.  An unconditionally stable maccormack method.  _Journal of Scientific Computing_ , 35(2):350–371, 2008.
  * Setaluri et al. [2014] Setaluri, R., Aanjaneya, M., Bauer, S., and Sifakis, E.  Spgrid: A sparse paged grid structure applied to adaptive smoke simulation.  _ACM Transactions on Graphics (TOG)_ , 33(6):1–12, 2014.
  * Shaham et al. [2021] Shaham, T. R., Gharbi, M., Zhang, R., Shechtman, E., and Michaeli, T.  Spatially-adaptive pixelwise networks for fast image translation.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pp. 14882–14891, 2021.
  * Sharp & Jacobson [2022] Sharp, N. and Jacobson, A.  Spelunking the deep: guaranteed queries on general neural implicit surfaces via range analysis.  _ACM Transactions on Graphics (TOG)_ , 41(4):1–16, 2022.
  * Sitzmann et al. [2020] Sitzmann, V., Martel, J., Bergman, A., Lindell, D., and Wetzstein, G.  Implicit neural representations with periodic activation functions.  _Advances in Neural Information Processing Systems_ , 33, 2020.
  * Smith et al. [2018] Smith, B., Goes, F. D., and Kim, T.  Stable neo-hookean flesh simulation.  _ACM Trans. Graph._ , mar 2018.
  * Sorkine & Alexa [2007] Sorkine, O. and Alexa, M.  As-rigid-as-possible surface modeling.  In _Symposium on Geometry processing_ , volume 4, pp. 109–116, 2007\.
  * Stam [1999] Stam, J.  Stable fluids.  In _Proceedings of the 26th annual conference on Computer graphics and interactive techniques_ , pp. 121–128, 1999.
  * Staniforth & Côté [1991] Staniforth, A. and Côté, J.  Semi-lagrangian integration schemes for atmospheric models—a review.  _Monthly weather review_ , 119(9):2206–2223, 1991\.
  * Takikawa et al. [2021] Takikawa, T., Litalien, J., Yin, K., Kreis, K., Loop, C., Nowrouzezahrai, D., Jacobson, A., McGuire, M., and Fidler, S.  Neural geometric level of detail: Real-time rendering with implicit 3d shapes.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pp. 11358–11367, 2021.
  * Taylor & Green [1937] Taylor, G. I. and Green, A. E.  Mechanism of the production of small eddies from large ones.  _Proceedings of the Royal Society of London. Series A-Mathematical and Physical Sciences_ , 158(895):499–521, 1937\.
  * Wandel et al. [2020] Wandel, N., Weinmann, M., and Klein, R.  Learning incompressible fluid dynamics from scratch–towards fast, differentiable fluid models that generalize.  _arXiv preprint arXiv:2006.08762_ , 2020.
  * Wang et al. [2021a] Wang, P., Liu, L., Liu, Y., Theobalt, C., Komura, T., and Wang, W.  Neus: Learning neural implicit surfaces by volume rendering for multi-view reconstruction.  _Advances in Neural Information Processing Systems_ , 34:27171–27183, 2021a.
  * Wang & Perdikaris [2021] Wang, S. and Perdikaris, P.  Long-time integration of parametric evolution equations with physics-informed deeponets.  _arXiv preprint arXiv:2106.05384_ , 2021.
  * Wang et al. [2021b] Wang, S., Wang, H., and Perdikaris, P.  Learning the solution operator of parametric partial differential equations with physics-informed deeponets.  _Science advances_ , 7(40):eabi8605, 2021b.
  * Wessels et al. [2020] Wessels, H., Weißenfels, C., and Wriggers, P.  The neural particle method–an updated lagrangian physics informed neural network for computational fluid dynamics.  _Computer Methods in Applied Mechanics and Engineering_ , 368:113127, 2020.
  * Wu & Zheng [2022] Wu, R. and Zheng, C.  Learning to generate 3d shapes from a single example.  _ACM Transactions on Graphics (TOG)_ , 41(6):1–19, 2022.
  * Xian et al. [2019] Xian, Z., Tong, X., and Liu, T.  A scalable galerkin multigrid method for real-time simulation of deformable objects.  _ACM Trans. Graph._ , 38(6), nov 2019.
  * Xie et al. [2021] Xie, Y., Takikawa, T., Saito, S., Litany, O., Yan, S., Khan, N., Tombari, F., Tompkin, J., Sitzmann, V., and Sridhar, S.  Neural fields in visual computing and beyond.  _arXiv preprint arXiv:2111.11426_ , 2021.
  * Yang et al. [2021] Yang, G., Belongie, S., Hariharan, B., and Koltun, V.  Geometry processing with neural fields.  _Advances in Neural Information Processing Systems_ , 34, 2021.
  * Yariv et al. [2020] Yariv, L., Kasten, Y., Moran, D., Galun, M., Atzmon, M., Ronen, B., and Lipman, Y.  Multiview neural surface reconstruction by disentangling geometry and appearance.  _Advances in Neural Information Processing Systems_ , 33:2492–2502, 2020.
  * Zehnder et al. [2021] Zehnder, J., Li, Y., Coros, S., and Thomaszewski, B.  Ntopo: Mesh-free topology optimization using implicit neural representations.  _Advances in Neural Information Processing Systems_ , 34:10368–10381, 2021.

##  Appendix A Implementation Details

###  A.1 Optimization

We solve our time-integration optimization problem (Equation 3) with the Adam
optimizer [30]. For all examples in our experiments, we set an initial
learning rate $\textbf{lr}_{0}$ and reduce it by a factor of $0.1$ if the loss
value does not decrease for $\textbf{iter}_{\text{p}}$ iterations. We stop the
optimization process when the learning rate is lower than
$\textbf{lr}_{\text{min}}$ or until it reaches a maximum of
$\textbf{iter}_{\text{max}}$ iterations. Specific values of these hyper-
parameters are described for each example below. We implemented our method
using the PyTorch library and performed experiments on an NVIDIA GeForce RTX
3090 GPU.

###  A.2 Advection Equation

For our advection example in Figure 2, the 1D spatial domain is
$\Omega=[-2,2]$. We consider the Dirichlet boundary condition, i.e., the
advected quantity at boundaries equals zero. Hence we set the boundary
constraint term in Equation 4 as

$\displaystyle{\boldsymbol{\mathcal{C}}}=||{u}^{n+1}({\boldsymbol{x}})||_{2}^{2},$
(15)

with the weighting factor $\lambda=1$. The initial condition for this example
is

$\displaystyle\hat{u^{0}}({\boldsymbol{x}})=e{-\frac{({\boldsymbol{x}}-\mu)^{2}}{2\sigma^{2}}},$
(16)

with $\mu=-1.5$ and $\sigma=0.1$. We set the optimization hyper-parameters
$\textbf{lr}_{0}=1e{-4}$, $\textbf{lr}_{\text{min}}=1e{-8}$,
$\textbf{iter}_{\text{p}}=500$ and $\textbf{iter}_{\text{max}}=20000$. For
each gradient descent iteration, we randomly sample $|\mathcal{M}|=5000$
points within the spatial domain $[-2,2]$. For this example, our method takes
$\sim 80s$ to compute per timestep, while the grid-based method (using the
same memory) takes $\sim 4e{-3}s$.

###  A.3 Incompressible Euler Equations

For our 2D fluid examples, the spatial domain is $\Omega=[-1,1]\times[-1,1]$.
We consider solid boundary conditions, i.e., the fluid cannot go through the
boundaries. Recall that we adopt the operator splitting scheme. Therefore, the
boundary constraint terms for the three sequential steps are

$\displaystyle\begin{split}{\boldsymbol{\mathcal{C}}}_{adv}=||{\boldsymbol{u}}_{adv\perp}^{n+1}({\boldsymbol{x}})||_{2}^{2}\\\
{\boldsymbol{\mathcal{C}}}_{pro}=||{\nabla}_{\perp}p^{n+1}({\boldsymbol{x}})||_{2}^{2}\\\
{\boldsymbol{\mathcal{C}}}_{cor}=||{\boldsymbol{u}}_{\perp}^{n+1}({\boldsymbol{x}})||_{2}^{2}\\\
\end{split}$ (17)

where $\perp$ indicates the perpendicular direction against the boundary. The
weighting factor $\lambda=1$.

Figure 10: Initial condition for the example in Figure 4. Left: velocity field
(Equation 19). Right: density field (Equation 22).

#### 2D Taylor-Green vortex

Standard 2D Taylor-Green is originally defined in domain
$[0,2\pi]\times[0,2\pi]$. We translate and scale the domain to
$[-1,1]\times[-1,1]$ such that the input range fits our MLP with the SIREN
activation [64]. Therefore, the initial condition for the velocity field
becomes

$\displaystyle\begin{split}\hat{{\boldsymbol{u}}^{0}}({\boldsymbol{x}})=(&\frac{1}{\pi}\sin[\pi(x+1)]\cos[\pi(y+1)],\\\
&-\frac{1}{\pi}\cos[\pi(x+1)]\sin[\pi[y+1]]).\end{split}$ (18)

After the simulation, we convert it back to domain $[0,2\pi]\times[0,2\pi]$
for evaluation and comparison. We set the optimization hyper-parameters
$\textbf{lr}_{0}=1e{-5}$, $\textbf{lr}_{\text{min}}=1e{-8}$,
$\textbf{iter}_{\text{p}}=500$ and $\textbf{iter}_{\text{max}}=20000$. The
size of the sample set is $|\mathcal{M}|=256^{2}$. For this example, our
method takes $\sim 10\text{min}$ to compute per timestep, while the grid-based
method (using the same memory) takes $\sim 0.03s$.

Table 6: Experiment setup for the elasticity examples. $\rho_{0}$ is the density. $\lambda$ and $\mu$ are the first and second lame parameters. $\alpha$ and $\beta$ are the number of hidden layers and the dimension of the hidden features. $\textbf{lr}_{0}$ and $\textbf{iter}_{\text{max}}$ are the initial learning rate and the maximum number of iterations. $\textbf{t}_{\text{avg}}$ is the average training time per time step. Note that the density $\rho_{0}$ and timestep size dt are reported as N/A for the quasistatic example Stretch (2D) (Figure 5).  Example | Dim | $|\mathcal{M}|$ | dt | $\rho_{0}$ | $\lambda$ | $\mu$ | $\alpha$ | $\beta$ | $\textbf{lr}_{0}$ | $\textbf{iter}_{\text{max}}$ |  $\textbf{t}_{\text{avg}}$(s)
---|---|---|---|---|---|---|---|---|---|---|---
Collision (2D) (Figure 7) | 2 | $100^{2}$ | 0.1 | $1e1$ | $2e1$ | $1e3$ | 3 | 68 | $1e{-5}$ | $1e{4}$ | $1.38e{2}$
Stretch (2D) (Figure 5) | 2 | $100^{2}$ | N/A | N/A | $1e0$ | $1e{3}$ | 3 | 68 | $1e{-4}$ | $5e{4}$ | $2.30e{3}$
Bunny (Figure 8) | 3 | $20^{3}$ | 0.1 | $1e0$ | $1e2$ | $1e3$ | 3 | 66 | $1e{-5}$ | $2e{4}$ | $1.70e{3}$
Spot (Figure 1) | 3 | $20^{3}$ | 0.1 | $1e0$ | $1e2$ | $1e3$ | 3 | 66 | $1e{-3}$ | $5e{3}$ | $1.74e{3}$
Lucy (Figure 13) | 3 | $20^{3}$ | 0.1 | $1e0$ | $1e3$ | $1e3$ | 3 | 128 | $1e{-4}$ | $2e{4}$ | $1.16e{3}$

#### Two vortices of different scale

For the example shown in Figure 4, the initial condition for the velocity
field is

$\displaystyle\hat{{\boldsymbol{u}}^{0}}({\boldsymbol{x}})=\begin{cases}\hat{{\boldsymbol{u}}^{1}}({\boldsymbol{x}}),&{\boldsymbol{x}}\in[-1,0]^{2}\\\
\hat{{\boldsymbol{u}}^{2}}({\boldsymbol{x}}),&{\boldsymbol{x}}\in[\frac{7}{4},1]^{2}\\\
(0,0),&\text{otherwise.}\end{cases}$ (19)

where

$\displaystyle\begin{split}\hat{{\boldsymbol{u}}^{1}}({\boldsymbol{x}})=&(\sin[2\pi(x+1)]\cos[2\pi(y+1)],\\\
&-\cos[2\pi(x+1)]\sin[2\pi(y+1)]),\end{split}$ (20)
$\displaystyle\begin{split}\hat{{\boldsymbol{u}}^{2}}({\boldsymbol{x}})=&(\sin[8\pi(x-\frac{7}{4})]\cos[8\pi(y-\frac{7}{4})],\\\
&-\cos[8\pi(x-\frac{7}{4})]\sin[8\pi(y-\frac{7}{4})]).\end{split}$ (21)

The density field that we advect is initialized as

$\displaystyle\hat{d^{0}}({\boldsymbol{x}})=\begin{cases}1&||2{\boldsymbol{x}}+1||\leq
0.5\text{ or }||8{\boldsymbol{x}}+7||\leq 0.5\\\
0&\text{otherwise.}\end{cases}$ (22)

Figure 10 visually illustrates the above initial conditions. We set the
optimization hyper-parameters $\textbf{lr}_{0}=1e{-5}$,
$\textbf{lr}_{\text{min}}=1e{-8}$, $\textbf{iter}_{\text{p}}=500$ and
$\textbf{iter}_{\text{max}}=20000$. The size of the sample set is
$|\mathcal{M}|=128^{2}$. For this example, our method takes $\sim
10\text{min}$ to compute per timestep.

#### Baseline setups

For PINN [54] and PINN-sub [31], we use the same MLP structure as ours
($\alpha=3$ hidden layers of width $\beta=32$, with SIREN activation). For
piDeepONet [74], we use $\alpha=3$ hidden layers of width $\beta=32$ with Tanh
activation for both the branch net and trunk net. Since they do not support
the operator-splitting scheme, we follow the previous literature [16] and use
the residual of incompressible Euler equation as the physics-informed training
objective,

$\displaystyle\begin{split}\mathcal{L}_{\theta_{\boldsymbol{u}},\theta_{p}}&=\underbrace{\sum_{\begin{subarray}{c}{\boldsymbol{x}}\in\mathcal{M}\subset\Omega\\\
t\in\mathcal{T}\end{subarray}}||\dfrac{\partial{\boldsymbol{u}}}{\partial
t}+{\boldsymbol{u}}\cdot{\nabla}{\boldsymbol{u}}+{\nabla}p||_{2}^{2}+||{\nabla}\cdot{\boldsymbol{u}}||_{2}^{2}}_{\text{PDE
residual}}\\\
&+\underbrace{\sum_{\begin{subarray}{c}{\boldsymbol{x}}\in\Omega\\\
t=0\end{subarray}}||{\boldsymbol{u}}-\hat{{\boldsymbol{u}}^{0}}||_{2}^{2}}_{\text{initial
condition}}+\underbrace{\sum_{\begin{subarray}{c}{\boldsymbol{x}}\in\partial\Omega\\\
t\in\mathcal{T}\end{subarray}}||{\boldsymbol{u}}_{\perp}({\boldsymbol{x}})||_{2}^{2}}_{\text{boundary
condition}}\end{split}$ (23)

where
${\boldsymbol{u}}={\boldsymbol{u}}_{\theta_{\boldsymbol{u}}}({\boldsymbol{x}},t)$
and $p={\boldsymbol{u}}_{\theta_{p}}({\boldsymbol{x}},t)$ are parameterized by
MLPs. The number of spatial samples used for each training iteration is the
same as ours. We train their models until convergence.

Another baseline (Ours-residual in Table 3) uses an implicit-time-discretized
version of this objective function for time integration, i.e.,

$\displaystyle\begin{split}{\boldsymbol{\mathcal{I}}}&=\sum_{\begin{subarray}{c}{\boldsymbol{x}}\in\mathcal{M}\subset\Omega\\\
t\in\mathcal{T}\end{subarray}}||\frac{{{\boldsymbol{u}}}^{n+1}-{{\boldsymbol{u}}}^{n}}{t}+{{\boldsymbol{u}}}^{n+1}\cdot{\nabla}{{\boldsymbol{u}}}^{n+1}+{\nabla}{p}^{n+1}||_{2}^{2}\\\
&+\sum_{\begin{subarray}{c}{\boldsymbol{x}}\in\mathcal{M}\subset\Omega\\\
t\in\mathcal{T}\end{subarray}}||{\nabla}\cdot{{\boldsymbol{u}}}^{n+1}||_{2}^{2}+\sum_{\begin{subarray}{c}{\boldsymbol{x}}\in\partial\Omega\\\
t\in\mathcal{T}\end{subarray}}||{{\boldsymbol{u}}}^{n+1}_{\perp}({\boldsymbol{x}})||_{2}^{2}.\end{split}$
(24)

Figure 11: Twisting test. Quasistatic simulation in 3D. The right end is
twisted 45 degrees.

Figure 12: Sampling convergence test. We use a different number of spatial
samples $|\mathcal{M}|$ and optimize for the same number of gradient descent
iterations. On the right, we further report the error with respect to the
reference result (using $|\mathcal{M}|=50^{3}$). As we increase the number of
spatial samples, the simulation result converges.

Figure 13: The statue collides with the ground and deforms. Our implicit
neural representation (green, right) is capable of capturing more fine
geometry details compared to the traditional tetrahedral mesh representation
(blue, left) under the same memory footprint.

###  A.4 Elastodynamic Equation

#### Initial and Boundary Conditions

For our 2D elasticity examples in Figure 5 and Figure 7, the 2D spatial domain
is $\Omega=[-1,1]\times[-1,1]$. For our 3D elasticity examples in Figure 12,
the 3D spatial domain is $\Omega=[-1,1]\times[-1,1]\times[-1,1]$. For our 2D
and 3D examples involving nonregular geometry (Figure 8, Figure 1 and Figure
13), the spatial domain is the interior of the shape, including the boundary.
The initial condition for all the elasticity examples is

$\displaystyle\begin{split}\hat{\boldsymbol{\phi}^{0}}({\boldsymbol{x}})=(0,0)&\text{
(2D)},\\\ \hat{\boldsymbol{\phi}^{0}}({\boldsymbol{x}})=(0,0,0)&\text{
(3D)}\end{split}$ (25)

The boundary constraint for elasticity examples involves positional
constraints or collision constraints. Positional constraints, or Dirichlet
boundary conditions, can be realized by defining the position of the
constraint set $\partial\Omega$ as the desired goal positions
$\overline{\boldsymbol{\phi}}_{\partial\Omega}$:

$\displaystyle{\boldsymbol{\mathcal{I}}}_{\text{pos}}$
$\displaystyle=\|{\boldsymbol{\phi}}^{n+1}_{\partial\Omega}-\overline{\boldsymbol{\phi}}_{\partial\Omega}\|^{2}_{2}.$
(26)

Collision constraints can be handled by adding unilateral constraints
dynamically and viewing the collision penalty force as an external force.
Specifically, for a colliding point ${\boldsymbol{q}}_{c}$, we first find the
closest surface point ${\boldsymbol{b}}_{c}$ with normal
${\boldsymbol{n}}_{c}$, and define our spring-like collision penalty force as:

$\displaystyle{\boldsymbol{f}}_{\text{col}}=k_{\text{col}}(({\boldsymbol{b}}_{c}-{\boldsymbol{q}}_{c})^{\top}{\boldsymbol{n}}_{c}){\boldsymbol{n}}_{c}.$
(27)

where $k_{\text{col}}$ is the ratio for the collision penalty force.

The corresponding collision energy can be defined as the work exerted by the
collision force:

$\displaystyle{\boldsymbol{\mathcal{I}}}_{\text{col}}=\rho_{0}{\boldsymbol{f}}_{\text{col}}^{T}{\boldsymbol{\phi}}^{n+1}.$
(28)

#### Experiment Setup

For all the 2D comparisons under the same memory usage, we use $\alpha=3$
hidden layers of width $\beta=68$ with SIREN activation function [64] for our
MLP, which takes the same memory (57 KB) as the FEM mesh (0.8K vertices, 1.5K
faces) and MPM point cloud (1.7K points) in use. We initialize the 2D
deformation field of the network to be zero using $|\mathcal{M}|=1000^{2}$
uniform and random samples. Then we train the network using
$|\mathcal{M}|=100^{2}$ uniform and random samples at each training iteration.
We use Bartels [33] and Taichi [24] to perform the FEM and MPM simulation,
respectively. We run our FEM and MPM comparison on CPU using a MacBook Pro
with the Apple M2 processor and 24GB of RAM.

For the 3D comparison under the same memory usage, for the bunny example
(Figure 8), we use $\alpha=3$ hidden layers of width $\beta=66$ with SIREN
activation function for our MLP, which takes the same memory (53 KB) as the
FEM mesh (0.5K vertices, 1.5K tetrahedra) in use. For the statue example
(Figure 13), we use $\alpha=3$ hidden layers of width $\beta=128$ with SIREN
activation function for our MLP, which takes the same memory (197 KB) as the
FEM mesh (2.0K vertices, 7.0K tetrahedra) in use. We initialize the 3D
deformation field of the network to be zero using $|\mathcal{M}|=100^{3}$
uniform and random samples. Then we train the network using
$|\mathcal{M}|=20^{3}$ uniform and random samples at each training iteration.
Here for simplicity, we use the mesh vertices as the uniform samples. We
further report all the parameters and experiment setup in Table 6. In
addition, we set the hyper-parameters $\textbf{iter}_{\text{p}}=800$ and
$\textbf{lr}_{\text{min}}=1e{-8}$ for all elasticity examples and assume a
constant density in the reference space.

The geometry (i.e., the undeformed shape) can be any representation (e.g.,
analytical, mesh, or level set), as long as it allows sampling within the
volume. For our examples in Figure 5 and Figure 7, the geometry (a square) is
represented analytically. For examples in Figure 8 and Figure 13, the geometry
is represented using the original high-resolution mesh.

For sampling of the shapes involving nonregular geometry, for simplicity, we
choose to use a triangle or tetrahedral mesh and perform sampling within the
volume _during the training_. An ideal alternative would be adopting the
implicit representation of the surface and performing rejection sampling based
on it.

For rendering, we sample a sufficient number of points from the undeformed
shape and evaluate the trained model at time $t$ on the sample positions to
predict their deformation. Thus the deformed shape $x^{t}$ at each time step
$t$ can be obtained by applying the deformation field $\boldsymbol{\phi}^{t}$
on the sample points of undeformed shape $x^{0}$, i.e.,
$x^{t}=x^{0}+\boldsymbol{\phi}^{t}(x^{0})$. For visualization, we only sample
the surface of the shape in 3D cases. Then we render the shape as a dense
point cloud.

##  Appendix B Additional Results

###  B.1 Elastodynamic Equation

In Figure 12, we demonstrate that our method exhibits volume-preserving
property on a 3D twisting example. In Figure 13, we provide another example
involving complex contact-induced deformations.

The L2 distance errors in Figure 6 and Table 4 are computed using the vertices
of the reference FEM mesh. Specifically, the L2 distance is defined between
the deformed positions of those vertices in the reference solution and in
our/compared FEM solution. In our solution, the deformed positions of those
vertices are queried via direct network inference. In the compared low-
resolution FEM solution, the deformed positions of those vertices are
calculated using the barycentric interpolation.

Finally, we demonstrate our method’s qualitative and quantitative convergence
when increasing the number of spatial samples for training. In Figure 12, we
compare the quasistatic stretching results when using a different number of
spatial samples (recall Section 3.2 Spatial Sampling), and report the error
with respect to the reference result (using $|\mathcal{M}|=50^{3}$).


# ARCOQ: Arabic Closest Opposite Questions Dataset

Sandra Rizkallah, Amir F. Atiya, and Samir Shaheen
Computer Engineering Department
Faculty of Engineering, Cairo University
Giza, Egypt

###### Abstract

This paper presents a dataset for closest opposite questions in Arabic
language. The dataset is the first of its kind for the Arabic language. It is
beneficial for the assessment of systems on the aspect of antonymy detection.
The structure is similar to that of the Graduate Record Examination (GRE)
closest opposite questions dataset for the English language. The introduced
dataset consists of 500 questions, each contains a query word for which the
closest opposite needs to be determined from among a set of candidate words.
Each question is also associated with the correct answer. We publish the
dataset publicly in addition to providing standard splits of the dataset into
development and test sets. Moreover, the paper provides a benchmark for the
performance of different Arabic word embedding models on the introduced
dataset.

Keywords: Arabic dataset; Closest Opposite; Antonymy Detection; Polarity

utf8

##  1 Introduction

Antonymy detection is vital for widespread applications including sentiment
analysis, contradiction detection, question answering, machine translation and
information retrieval [1, 2, 3, 4]. In sentiment analysis, one needs to
differentiate between positive and negative attitudes. Some methods were
developed for this purpose making use of antoyms and antonymy detection. In
order to assess them, it would be a good idea to have a benchmark dataset for
antonyms. This would be helpful in comparing different antonymy detection
methods. It can also be beneficial in assessing the capability of the methods
in selecting the most accurate antonym. Educational applications are another
example. As a step to automatic answer generation, we need to be able to have
an approach that can accurately select the correct answer from a set of
choices in multiple choice questions.

Assessing the capability of a system to detect antonyms is addressed in the
English language, for example in [5, 6] a dataset is gathered from Graduate
Record Examination (GRE) closest opposite questions. The mentioned dataset is
used by a number of works addressing the English language. In [7], Multi-
Relational Latent Semantic Analysis (MRLSA) is presented showing the
performance of the model on identifying antonyms relations. In [8], a
Bayesian-based model is proposed where neural word embeddings are generated to
differentiate synonyms and antonyms. In [9], word embeddings are trained in a
way such that antonyms can be captured. In [10], algorithms are evaluated for
detecting antonymy. In [11], vector space is fune-tuned to include antonymy
relations. In [12, 13], new spherical word embeddings are learned to capture
antonymy. In [14], an opinion mining approach is presented utilizing antonymy
relations. In [15], contextual relations are addressed with a focus on
antonymy.

On the other hand, few works address antonymy in the Arabic language. In [16],
the study investigates antonymy use in the Modern Standard Arabic (MSA) texts.
In [17], ontological lexicon enrichment is performed through semi-automated
extraction of lexical relations focusing on antonyms. In [18], a pattern-based
approach is used with the aim of extending the Arabic WordNet with more
relations specially, antonyms. In [19], Arabic word vectors are embedded in a
unit sphere to enable antonymy detection where antonyms vectors lie on
opposite poles of the sphere. In all the aforementioned works, no unified
benchmark is used to measure the performance of the proposed systems.

To our knowledge, datasets that allow for measuring the performance of
antonymy detection do not exist for the Arabic language. To fill this gap, we
present a new dataset for the Arabic language that can be used for assessing
the performance of systems on the aspect of antonymy detection. The dataset
contains 500 questions in the form of multiple choice questions. The target is
to choose the closest opposite word (antonym) to a query word given in the
question from among the given candidate choices. Each question is structured
as follows: “query_word: candidate_answer1 candidate_answer2 candidate_answer3
:: correct_answer”. Moreover, we present a benchmark that records the
performance of different Arabic word embedding models in detecting antonyms on
the introduced dataset. We considered word embedding methods because these are
the dominant approach in antonymy detection. The dataset is published publicly
111https://github.com/SandraRizkallah/ARCOQ-dataset in order to urge the
research community to build upon the introduced benchmark and create new
benchmarks for the Arabic language.

The rest of the paper is structured as follows: Section 2 discusses the data
collection and processing procedures. In Section 3, the benchmark is
presented, explaining the methods used and discussing the obtained results on
the introduced dataset. Finally, Section 4 shows the conclusion and future
work.

##  2 Data Collection and Processing

We collected the dataset from the following sources:

  1. 1.

Educational exercises
Choosing the closest opposite for a given Arabic word from among a set of
candidate words is a question that is commonly included in the educational
exercises and exams for the Arabic language subject. We have explored examples
of the Arabic subject exams and exercises for the different Egyptian
educational stages: primary, preparatory and secondary. Such exercises are
provided by the Egyptian Ministry of Education and Technical Education through
Egyptian Knowledge Bank (EKB) platform [20].

  2. 2.

Translating GRE closest opposite questions
In [5, 6] a GRE closest opposite questions English dataset is introduced. We
have automatically translated all the questions in this dataset from English
to Arabic using “google-trans-new” python package [21]. After automatic
translation, we manually removed the questions that we considered to be not
perfectly translated.

  3. 3.

Native Arabic speakers
We asked native Arabic speakers (specifically 4 speakers) to generate
questions of the required type and form i.e. one query Arabic word with three
candidate words from which the closest opposite of the query word should be
chosen.

After collecting the questions as explained above, we perform some processing
and checks:

  * •

Remove any diacritics (tashkeel).

  * •

Remove any elongations in the letters (tatweel).

  * •

Remove any duplicated questions.

  * •

Make sure that the correct answer is among the given choices in the question
(this is a sanity check).

Finally, the dataset consists of 500 questions, each in the following form:

“query_word: candidate_answer1 candidate_answer2 candidate_answer3 ::
correct_answer”

These questions are further sorted alphabetically based on the query word.
Moreover, we provide standard splits for the dataset into development and test
sets. We randomly split the dataset into 20% development set and 80% test set:

  * •

Dev 100: contains 100 questions

  * •

Test 400: contains 400 questions

We followed the lead of the GRE dataset [5, 6] in having a development set to
test set ratio of 1:4. The reason is that there are not many parameters to
tune, for example in most methods pre-trianed embedding vectors are available.

Figure 1 shows example questions from the ARCOQ dataset we established.

Figure 1: Example Questions from ARCOQ dataset

##  3 Benchmark

We established a benchmark using our ARCOQ dataset. The performance of
different Arabic word embedding models in detecting antonymy is measured. The
benchmark includes results for the development set, the test set and the whole
dataset. The metrics used for the performance evaluation are those introduced
in [5, 6] as follows:

$Precision=\frac{\mbox{number of correctly answered questions}}{\mbox{number
of answered questions}}$ (1) $Recall=\frac{\mbox{number of correctly answered
questions}}{\mbox{total number of questions}}$ (2)
$Fscore=\frac{(2*Precision*Recall)}{(Precision+Recall)}$ (3)

A question is answered by the model only if the model has vectors for all the
words in the question (query and all candidate answers).

The benchmark includes the major Arabic word embedding models:

  1. 1.

Polyglot (2013) [22]

  2. 2.

Altowayan (2016) [23]

  3. 3.

AraVec: full grams CBOW 100 twitter (2017) [24]

  4. 4.

AraVec: full grams SG 100 twitter (2017) [24]

  5. 5.

fastText (2018) [25]

  6. 6.

ArWordVec: SG_300_3 (2020) [26]

  7. 7.

ArWordVec: SG_300_5 (2020) [26]

  8. 8.

AraBERT (2020) [27]

  9. 9.

ArSphere (2021) [19]

The models from 1–8 are distributional models that are pre-trained on large
Arabic text corpora. In such models, antonyms vectors can be placed far away
from each other (small similarity between vectors) or close to each other
(high similarity between vectors). The latter is because antonyms often occur
in the same context. On the other hand, Arsphere is a semi-supervised model
that embeds antonyms vectors at opposite poles of the sphere (negative
similarity between vectors).

For answering a given question in the dataset, the following procedure is
followed:

  1. 1.

Get the vector of the query word from the designated model.

  2. 2.

Get the vector of each candidate answer from the designated model.

  3. 3.

Get the similarity scores between the query word and each candidate answer:

     * •

Normalize all the obtained vectors.

     * •

Compute the dot product between the normalized vector of the query word and
the normalized vector of each candidate answer.

     * •

Having 3 computed similarity scores, choose the answer of the question as the
word with highest/lowest similarity score with the query word.

     * •

Highest or lowest similarity score?, we performed two experiments:

       * –

Experiment 1: Lowest similarity score is chosen for all the models. Such
choice is made since in embedding models the similarity between vectors
reflect the underlying semantics so synonyms will have high similarity scores
while antonyms or unrelated words conversely have low similarity scores.

       * –

Experiment 2: Highest similarity score is chosen for models 1–8 while also the
lowest similarity score is chosen for ArSphere. This choice is made because in
distributional models antonyms often occur in the same context. Therefore,
vectors of synonyms as well as antonyms will mostly have high similarity
scores.

The choice for ArSphere model is always the lowest similarity score. This is
because this model is specifically designed to incorporate antonym relations
where vectors of antonyms are designed by having negative similarity scores
(numbers below zero).

  4. 4.

Finally, the performance metrics are computed over all the questions.

###  3.1 Results

Table 1 and Table 2 include the results of Experiments 1 and 2 respectively.

Table 1: Benchmark Experiment1

Approach | Dev 100 | Test 400 | All 500
---|---|---|---
| Precision | Recall | Fscore | Precision | Recall | Fscore | Precision | Recall | Fscore
Polyglot (2013) | 0.23 | 0.14 | 0.18 | 0.22 | 0.14 | 0.17 | 0.22 | 0.14 | 0.17
Altowayan (2016) | 0.24 | 0.19 | 0.21 | 0.25 | 0.21 | 0.23 | 0.25 | 0.2 | 0.22
AraVec_cbow (2017) | 0.17 | 0.17 | 0.17 | 0.22 | 0.22 | 0.22 | 0.19 | 0.19 | 0.19
AraVec_sg (2017) | 0.18 | 0.18 | 0.18 | 0.19 | 0.19 | 0.19 | 0.19 | 0.19 | 0.19
fastText (2018) | 0.1 | 0.1 | 0.1 | 0.15 | 0.15 | 0.15 | 0.14 | 0.14 | 0.14
ArWordVec_SG3 (2020) | 0.17 | 0.07 | 0.1 | 0.16 | 0.07 | 0.1 | 0.16 | 0.07 | 0.1
ArWordVec_SG5 (2020) | 0.12 | 0.05 | 0.07 | 0.15 | 0.07 | 0.09 | 0.14 | 0.06 | 0.09
AraBERT (2020) | 0.27 | 0.27 | 0.27 | 0.25 | 0.25 | 0.25 | 0.25 | 0.25 | 0.25
ArSphere (2021) | 0.93 | 0.93 | 0.93 | 0.76 | 0.76 | 0.76 | 0.8 | 0.8 | 0.8

Table 2: Benchmark Experiment2

Approach | Dev 100 | Test 400 | All 500
---|---|---|---
| Precision | Recall | Fscore | Precision | Recall | Fscore | Precision | Recall | Fscore
Polyglot (2013) | 0.57 | 0.34 | 0.43 | 0.51 | 0.33 | 0.4 | 0.52 | 0.33 | 0.41
Altowayan (2016) | 0.47 | 0.37 | 0.41 | 0.4 | 0.33 | 0.36 | 0.42 | 0.33 | 0.37
AraVec_cbow (2017) | 0.52 | 0.51 | 0.52 | 0.51 | 0.5 | 0.5 | 0.51 | 0.5 | 0.5
AraVec_sg (2017) | 0.59 | 0.58 | 0.59 | 0.56 | 0.56 | 0.56 | 0.57 | 0.56 | 0.56
fastText (2018) | 0.63 | 0.62 | 0.62 | 0.59 | 0.58 | 0.59 | 0.6 | 0.59 | 0.59
ArWordVec_SG3 (2020) | 0.64 | 0.27 | 0.38 | 0.6 | 0.27 | 0.37 | 0.61 | 0.27 | 0.38
ArWordVec_SG5 (2020) | 0.57 | 0.24 | 0.34 | 0.61 | 0.28 | 0.38 | 0.6 | 0.27 | 0.37
AraBERT (2020) | 0.49 | 0.49 | 0.49 | 0.45 | 0.45 | 0.45 | 0.46 | 0.46 | 0.46
ArSphere (2021) | 0.93 | 0.93 | 0.93 | 0.76 | 0.76 | 0.76 | 0.8 | 0.8 | 0.8

###  3.2 Comments on the Results

All the models performed better in Experiment 2. This is because
distributional models tend to place vectors of words that occur in the same
context in close proximity (having high similarity scores) which is mostly the
case for antonyms. It can be shown that the best performing model is ArSphere
which is specifically trained to detect antonymy. The other models are not
specifically trained to detect antonymy and they perform well in other
relations such as detecting synonymy. From the results, it is also shown that
the second best performing model is fastText [25] followed by AraVec SG model
[24].

##  4 Conclusion and Future Work

We have established the first closest opposite questions dataset for the
Arabic language: “ARCOQ”. The dataset is beneficial for antonymy detection
specifically and question answering in general. Moreover, a benchmark is
developed comparing the major existing Arabic embedding models on ARCOQ. The
benchmark shows that there still exists room for improving Arabic word
embeddings on the aspect of antonymy. In the future, the presented benchmark
can be built upon including other models either than embedding models.

## References

  * [1] Silke Scheible, Sabine Schulte Im Walde, and Sylvia Springorum.  Uncovering distributional differences between synonyms and antonyms in a word space model.  In Proceedings of the Sixth International Joint Conference on Natural Language Processing, pages 489–497, 2013.
  * [2] Natalia Silveira.  Learning to identify antonyms.  URL: http://cs229. stanford. edu/proj2013/Silveira-LearningToIdentifyAntonymy. pdf, 2013.
  * [3] Luyang Li, Bing Qin, and Ting Liu.  Contradiction detection with contradiction-specific word embedding.  Algorithms, 10(2):59, 2017.
  * [4] Van-Tan Bui, Khac-Quy Dinht, and Phuong-Thai Nguyen.  Vietnamese antonyms detection based on specialized word embeddings using semantic knowledge and distributional information.  In 2020 12th International Conference on Knowledge and Systems Engineering (KSE), pages 159–164. IEEE, 2020.
  * [5] Saif Mohammad, Bonnie Dorr, and Graeme Hirst.  Computing word-pair antonymy.  In Proceedings of the 2008 Conference on Empirical Methods in Natural Language Processing, pages 982–991, 2008.
  * [6] Saif M Mohammad, Bonnie J Dorr, Graeme Hirst, and Peter D Turney.  Computing lexical contrast.  Computational Linguistics, 39(3):555–590, 2013.
  * [7] Kai-Wei Chang, Wen-tau Yih, and Christopher Meek.  Multi-relational latent semantic analysis.  In Proceedings of the 2013 Conference on Empirical Methods in Natural Language Processing, pages 1602–1612, 2013.
  * [8] Jingwei Zhang, Jeremy Salwen, Michael Glass, and Alfio Gliozzo.  Word semantic representations using bayesian probabilistic tensor factorization.  In Proceedings of the 2014 Conference on Empirical Methods in Natural Language Processing (EMNLP), pages 1522–1531, 2014.
  * [9] Masataka Ono, Makoto Miwa, and Yutaka Sasaki.  Word embedding-based antonym detection using thesauri and distributional information.  In Proceedings of the 2015 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, pages 984–989, 2015.
  * [10] Dean J Jones and Gunjan Mansingh.  Not just dissimilar, but opposite: An algorithm for measuring similarity and oppositeness between words.  In 2016 International Conference on Data Science and Engineering (ICDSE), pages 1–6. IEEE, 2016.
  * [11] Ivan Vulić.  Injecting lexical contrast into word vectors by guiding vector space specialisation.  In Proceedings of The Third Workshop on Representation Learning for NLP, pages 137–143, 2018.
  * [12] Sandra Rizkallah, Amir F Atiya, and Samir Shaheen.  A polarity capturing sphere for word to vector representation.  Applied Sciences, 10(12):4386, 2020.
  * [13] Sandra Rizkallah, Amir F Atiya, and Samir Shaheen.  New vector-space embeddings for recommender systems.  Applied Sciences, 11(14):6477, 2021.
  * [14] Sandra Rizkallah, Amir F Atiya, and Samir Shaheen.  Learning spherical word vectors for opinion mining and applying on hotel reviews.  In International Conference on Intelligent Systems Design and Applications, pages 200–211. Springer, 2020.
  * [15] Wenqiang Lei, Yisong Miao, Runpeng Xie, Bonnie Webber, Meichun Liu, Tat-Seng Chua, and Nancy F Chen.  Have we solved the hard problem? it’s not easy! contextual lexical contrast as a means to probe neural coherence.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 35, pages 13208–13216, 2021.
  * [16] Rukayah AlHedayani.  Antonymy in modern standard Arabic.  PhD thesis, University of Sussex, 2016.
  * [17] Maha Al-Yahya, Sawsan Al-Malak, and Luluh Aldhubayi.  Ontological lexicon enrichment: The badea system for semi-automated extraction of antonymy relations from arabic language corpora.  Malaysian Journal of Computer Science, 29(1):56–73, 2016.
  * [18] Mohamed Ali Batita and Mounir Zrigui.  The enrichment of arabic wordnet antonym relations.  In International Conference on Computational Linguistics and Intelligent Text Processing, pages 342–353. Springer, 2017.
  * [19] Sandra Rizkallah, Amir Atiya, Samir Shaheen, and Hossam ElDin Mahgoub.  Arsphere: Arabic word vectors embedded in a polar sphere.  Working Paper, 2021.
  * [20] Egyptian Knowledge Bank (EKB).  https://lms.ekb.eg/repository/discovery?sort=recommended\&strict=0.  [Online; accessed 25-July-2021].
  * [21] google-trans-new.  https://pypi.org/project/google-trans-new/.  [Online; accessed 12-July-2021].
  * [22] Rami Al-Rfou, Bryan Perozzi, and Steven Skiena.  Polyglot: Distributed word representations for multilingual nlp.  arXiv preprint arXiv:1307.1662, 2013.
  * [23] A Aziz Altowayan and Lixin Tao.  Word embeddings for arabic sentiment analysis.  In 2016 IEEE International Conference on Big Data (Big Data), pages 3820–3825. IEEE, 2016.
  * [24] Abu Bakr Soliman, Kareem Eissa, and Samhaa R El-Beltagy.  Aravec: A set of arabic word embedding models for use in arabic nlp.  Procedia Computer Science, 117:256–265, 2017.
  * [25] Edouard Grave, Piotr Bojanowski, Prakhar Gupta, Armand Joulin, and Tomas Mikolov.  Learning word vectors for 157 languages.  In Proceedings of the International Conference on Language Resources and Evaluation (LREC 2018), 2018.
  * [26] Mohammed M Fouad, Ahmed Mahany, Naif Aljohani, Rabeeh Ayaz Abbasi, and Saeed-Ul Hassan.  Arwordvec: efficient word embedding models for arabic tweets.  Soft Computing, 24(11):8061–8068, 2020.
  * [27] Wissam Antoun, Fady Baly, and Hazem Hajj.  Arabert: Transformer-based model for arabic language understanding.  In LREC 2020 Workshop Language Resources and Evaluation Conference 11–16 May 2020, page 9.