# Portable, Massively Parallel Implementation of a Material Point Method for
Compressible Flows

Paolo Joseph Baioni(1), Tommaso Benacchio(2),
Luigi Capone(3), Carlo de Falco(1)

###### Abstract

The recent evolution of software and hardware technologies is leading to a
renewed computational interest in Particle-In-Cell (PIC) methods such as the
Material Point Method (MPM). Indeed, provided some critical aspects are
properly handled, PIC methods can be cast in formulations suitable to the
requirements of data locality and fine-grained parallelism of modern hardware
accelerators as Graphics Processing Units (GPUs). Such a rapid and continuous
technological development increases also the importance of generic and
portable implementations. While continuum mechanics simulations have already
shown the capabilities of MPM on a wide range of phenomena, the use of the
method in compressible fluid dynamics is less frequent, especially in the
supersonic regime.

In this paper we present a portable, highly parallel, GPU based MPM solver for
compressible gas dynamics. The implementation aims to reach a good compromise
between portability and efficiency and to give a first assessment of the
potential of this approach in reproducing high speed gas flows, also taking
into account solid obstacles. The proposed model constitutes a new step
towards the realization of a monolithic MPM solver for Fluid-Structure
Interaction (FSI) problems at all Mach numbers up to the supersonic regime.

(1) MOX – Modelling and Scientific Computing
Dipartimento di Matematica, Politecnico di Milano
Piazza Leonardo da Vinci 32, 20133 Milano, Italy
<EMAIL_ADDRESS>

(2) Danish Meteorological Institute
Sankt Kjelds Plads 11
2100 Copenhagen, Denmark
<EMAIL_ADDRESS>

(3) Leonardo Labs, Leonardo S.p.A.
Torre Fiumara, Via Raffaele Pieragostini, 80, 16149 Genova, Italy
<EMAIL_ADDRESS>

Keywords: Material Point Method; GPUs; Supersonic Gas Dynamics; Performance
Portability; High-Performance Computing.

This work has not yet been peer-reviewed and is provided by the contributing
authors as a means to ensure timely dissemination of scholarly and technical
work on a noncommercial basis. Copyright and all rights therein are maintained
by the authors or by other copyright owners. It is understood that all persons
copying this information will adhere to the terms and constraints invoked by
each author’s copyright. This work may not be reposted without explicit
permission of the copyright owner. This work has been submitted for
publication. Copyright in this work may be transferred without further notice.

##  1 Introduction and Motivation

Hardware and software technology evolution has led to a renewed interest in
particle-based numerical techniques for differential problems, such as the
Material Point Method (MPM, [45, 56]), that was developed stemming from the
Particle In Cell method Fluid Implicit Particle Method [26, 7] to treat
history-dependent materials, thanks to their suitability to modern hardware
accelerators requirements [25, 20]. Indeed, parallel computing is a key
capability for effective advancement of industrial-grade, large-scale
computations, particularly in the fields of high-speed compressible
computational fluid dynamics and Fluid-Structure Interaction (FSI, [32])
simulations. On the other hand, rapid hardware and software evolution makes
researchers lean towards generic and portable implementations, which enables
developers to adapt fast to sudden technology innovations.
In view of exponentially increasing High Performance Computing (HPC) resources
and the method’s accurate handling of a wide range of phenomena, MPM is
becoming increasingly relevant for high-speed CFD and FSI simulations. In MPM,
continuum media are discretized by means of Lagrangian particles, the material
points, each one storing every state variable evaluated at its position,
enabling the description of a wide range of materials [23, 43, 29] as well as
large deformations [1, 52, 16] and fracture phenomena [53]. Notwithstanding
its Lagrangian character, MPM also employs a background Cartesian grid to
compute differential quantities and solve the motion equation, thus mediating
the particle-particle interactions, and taking advantage of both Eulerian and
Lagrangian approaches. These features make it an interesting candidate for the
development of a highly parallel, Graphics Processing Units (GPUs)-based FSI
solver. Indeed, in order to reach a high degree of efficiency, GPU-based
solvers require a large degree of data locality and fine-level parallelism.
MPM implementations can naturally achieve these, since:

  * •

every GPU thread can manage up to very few grid cells or particles;

  * •

the link between cells and particles can be easily built thanks to the
structured character of the grid;

  * •

accuracy in the discretisation of continuum bodies is obtained by means of the
Lagrangian particles, so that the grid does not have to match the solid
boundary.

Thus, there is no need of computationally expensive, unstructured, moving or
deforming meshes and a structured cartesian grid can be used instead,
simplifying also techniques such as Adaptive Mesh Refinement and Domain
Decomposition in a parallel context [47, 54].
Since MPM first introduction, various enhancements of the method were
proposed, solving some of its drawbacks such as grid-crossing error [49] and
volumetric locking [15], and improving its accuracy and conservation
properties. To address instabilities which arise when particles cross the
grid, variants of the method were proposed, such as GIMP [4], CPDI [39], iMPM
[46] and IGA-MPM [34], and numerical analysis studies on the use of B-Splines
as basis functions in MPM were conducted [42]. An affine and a polynomial
mapping was later added to the B-Splines MPM [30, 21] resulting in an angular
momentum-conserving model. Other studies explored least squares techniques
[55] enabling the recovery of the affine/polynomial splines MPM [28], the
effect of spatial and temporal discretization errors, and the application of
symplectic integrators [5]. In light of these developments, MPM represents an
effective method for the simulation of a wide range of materials and
phenomena.
However, in spite of recent implementations for computer graphics applications
[22, 54], MPM implementations on GPUs have not yet been applied to
compressible fluid dynamics. In fact, even when addressing gas dynamics and
solids mechanics, MPM has been used only for the solids, while the
compressible fluid dynamic problem has been solved through a second order WENO
scheme [11]; moreover, such implementation is run on CPUs. Existing studies on
MPM for compressible flows focused on numerical convergence and error analysis
[51] and on modelling and expanding the range of applications for the method
[58, 44]. Nonetheless, to the best of the authors’ knowledge, no studies exist
that are dedicated to the GPU based, High Performance Computing aspects of MPM
implementations for strongly compressible fluid dynamics problems.
In previous work [2] we ported an existing MPM algorithm [58, 57, 12] to
CUDA-C [35], adapting it to GPUs architecture and obtaining promising speed-
ups with respect to a C++ CPU version. The CUDA-C programming model was chosen
because of its known performance capabilities on NVIDIA A100 GPUs [36, 8].
However, the model had shortcomings in its low-level character and specificity
to the underlying architecture. In addition, while it can be made to work on,
e.g., AMD GPUs via HIP with minimal changes, it cannot be ported to different
architectures such as CPUs without substantial code modifications.
One of NVIDIA’s solutions for performance portability involves its own C++
compiler nvc++ [19, 37], using which one rewrites algorithm steps relying on
C++ Standard Template Library (STL), specifying a parallel execution policy.
To give an example, in this framework it is possible to parallelize the
following for loop:

⬇

1 for (index_t i = 0; i < n; ++i)

2 {do some operations}

into a call to a for_each with a proper lambda or functor, e.g.:

⬇

1 std::for_each (std::execution::par, first_counting_iterator,
last_counting_iterator, [](index_t i){do some operations});

Although the approach does not prevent the developer from introducing data-
races or deadlocks, the code can be written only once, and then, depending on
the option specified during compilation, run in parallel on multicore CPUs or
GPUs.
However, nvc++ high-level features come at a price, as the developer loses the
control on memory handling and, in particular, on when Host-to-Device and
Device-to-Host memory copies are performed. The code runs on host, any
invocation to an algorithm with a parallel execution policy expresses the
developer’s preference to run that algorithm in parallel. If compiled for
GPUs, memory is copied to the device, computations are performed and then data
is copied back to host. In physics and in engineering dynamics simulations,
where each time step depends on the previous one, this means that the maximum
reachable level of parallelization corresponds to a single call to a parallel
algorithm inside a traditional for loop for time advancing. Since memory
fetches can be significantly more costly than computations, this is
suboptimal.
The present contribution expands the scope of [2] and overcomes the
performance issues with nvc++ by designing a portable MPM code for transonic
and supersonic compressible gas dynamics based on Thrust [38], a C++ template
library for CUDA based on the Standard Template Library (STL). In addition to
nvc++ features, Thrust introduces two types of STL like vectors, one for the
host (tipically, the CPU) and one for the device (tipically, the GPU); this
allows the user to specify when memory should be copied via explicit calls.
For example, it’s possible to perform a single Host to Device copy when
reading the initial data file, and one Device To Host copy whenever saving
partial results on output data files, strongly reducing execution time related
to memory fetches.
Moreover, as additional portability-enhancing features included in the
implementation:

  * •

Thrust syntax is mostly compatible with STL syntax, adding just a few define
preprocessor macros or templates. Thanks to this, almost the same source code
can be compiled for parallel execution on GPUs with NVIDIA CUDA compiler nvcc,
on AMD GPUs with the AMD compiler for ROCmThrust, and on CPUs with g++ and
clang++ compilers (when it is linked with a parallel library such as TBB,
Intel’s Threading Building Blocks).

  * •

Thrust itself allows for modifications to the default device and host system,
so that it can be run serially and in parallel, with TBB or OpenMP on CPUs for
the host, and also on the GPU via CUDA-C++ for the device only. The code can
be compiled and executed on a CPU-only machine as well, provided CUDA is
installed and THRUST_ DEVICE_SYSTEM is defined to be THRUST_DEVICE_SYST
EM_{CPP,OMP,TBB} at compilation stage.

Notwithstanding the portability objectives, we designed the code for
optimization on modern graphics hardware accelerators and made further
modifications to the original algorithm [12, 58]. For example, we:

  * •

added treatment of obstacles in gas dynamics simulations,

  * •

added a particle re-ordering step,

  * •

postponed particles movement to the end of the time loop,

and thus values of basis function evaluated at previous particle positions
don’t need to be stored and can be computed on the fly instead, which is
cheaper on modern hardware accelerators.
Specific development choices are described that enable optimization on modern
GPU-based HPC architectures and adaptability to different hardware, as well as
generality of the approach for the purposes of its portability also towards
different software. The accuracy and performance of the implementation are
evaluated on a number of compressible flow benchmarks and results are compared
with available approaches in the literature.
The paper is structured as follows. Section 2 contains the analytical and
numerical formulation of the method. Section 3 describes the chosen algorithm
and the current and implementation choices. Section 4 illustrates the results
obtained on benchmark test cases, and the final section 5 contains a
discussion of the results, draws the conclusions of the performed work, and
prospects future developments.

##  2 Model Formulation

###  2.1 Model Equations and Weak Formulation

As mentioned, MPM takes advantages of both the Eulerian and the Lagrangian
perspective; specifically, the Lagrangian particles discretisation approach
affords flexibility in modelling various kind of phenomena and materials,
while the Eulerian grid discretisation is beneficial for scalability as the
grid allows the mediation of particle-particle interactions. For the
mathematical modelling we choose the Lagrangian framework, and we model high
speed gases as inviscid fluids in a domain $\Omega$ governed by the
compressible Euler equations [24, 57]:

$\begin{cases}\dfrac{d\rho}{dt}+\rho\nabla\cdot\mathbf{v}=0\\\\[5.69054pt]
\rho\dfrac{d\mathbf{v}}{dt}=\nabla\cdot\bm{\sigma}+\rho\mathbf{b}\\\\[5.69054pt]
\rho\dfrac{de}{dt}=\bm{\sigma}:\bm{\dot{\varepsilon}}\end{cases}$ (1)

where $\rho$, $\mathbf{v}$, $\bm{\sigma}$, $\rho\mathbf{b}$, $p$, $e$ and
$\bm{\dot{\varepsilon}}$ represent respectively the mass density, the
velocity, the stress tensor, the external forces density, the pressure, the
specific internal energy and the strain rate tensor of the fluid.
In this work, we:

  * •

enforce mass conservation by keeping constant the number of Lagrangian
material points (particles) along the simulations, as well as each particle
mass,

  * •

choose to model the fluid as a bi-atomic perfect gas, governed by the state
equation

$p=(\gamma-1)\rho e$ (2)

where $e$ is the specific internal energy and $\gamma$ is the ratio between
specific heat at constant pressure and constant volume,
$\gamma=c_{p}/c_{v}=1.4$ for a bi-atomic gas.

To derive the final discretisation of the equations, a weak form of the force
balance expressed by the motion equation is derived [45, 44]:

$\int_{\Omega}\rho\mathbf{\dot{v}}\cdot\mathbf{w}d\Omega=-\int_{\Omega}\bm{\sigma}:\nabla\mathbf{w}d\Omega+\int_{\partial\Omega}\bm{\tau}\cdot\mathbf{w}d\Gamma+\int_{\Omega}\rho\mathbf{b}\cdot\mathbf{w}d\Omega$
(3)

where, $\mathbf{\dot{v}}$ denotes the total, or Lagrangian, velocity
derivative, $\mathbf{w}$ is a generic test function and $\bm{\tau}$ the stress
on the surface. It is worth noting that it is possible to derive the same
formulation also by adopting the standard Finite Element procedure, provided
material points are considered as moving quadrature points [17]. In fact,
standard MPM can be seen as a modification of traditional FEM, in which the
integration points are not related to the Finite Element mesh and are
Lagrangian, and thus suitable for history dependent models too [12].
This dynamically comprehensive yet relatively simple equation set is chosen as
it enables testing the performance portability of the MPM implementation on
compressible gas dynamics benchmarks [50, 41] and comparison with existing
approaches [58, 44].
In the following, we simplify further the model assuming the fluid to be
inviscid, and thus the stress tensor to be due to pressure only, and we
neglect external forces. Euler equations become:

$\begin{cases}\dfrac{d\rho}{dt}+\rho\nabla\cdot\mathbf{v}=0\\\\[5.69054pt]
\rho\dfrac{d\mathbf{v}}{dt}=-\nabla p\\\\[5.69054pt]
\rho\dfrac{de}{dt}=-p\nabla\cdot\mathbf{v}\end{cases}$ (4)

The continuity and energy conservation equations will be directly discretized
and used to update particle properties in the algorithm (see section 3 and
algorithm 1), while the discretisation of the weak form of the motion equation
is treated in more depth in the next section.

###  2.2 Numerical Method

The link between the Eulerian grid and the Lagrangian particles is given by
the Particle-to-Grid (P2G) and Grid-to-Particle (G2P) mappings, which are the
building blocks of the computational algorithm, and that can be derived from
the discretization of the weak form (3). The latter can be obtained expressing
the mass density via of Dirac delta functions [44]:

$\rho(\mathbf{x},t)=\sum_{k=1}^{N_{p}}M_{p}\delta(\mathbf{x}-\mathbf{x_{p}})$
(5)

Here and in the following, when it is useful to distinguish, grid quantities
are written in lower case Roman letters, while particle ones are written in
either upper case Roman or Greek letters.
Particle-to-Grid and Grid-to-Particle mappings read [45]:

$a=\mathrm{P2G}(A;X,Y):=\left[\sum_{k=0}^{Np-1}A_{k}\
u_{i}(X_{k},Y_{k})\right]$ (6)
$A=\mathrm{G2P}(a;X,Y):=\left[\sum_{i=0}^{(nx+1)(ny+1)-1}a_{i}\
u_{i}(X_{k},Y_{k})\right]$ (7) $\nabla
a=\mathrm{P2GD}(A;X,Y):=\left[\sum_{k=0}^{Np-1}A_{k}\ \nabla
u_{i}(X_{k},Y_{k})\right]$ (8) $\nabla
A=\mathrm{G2PD}(a;X,Y):=\left[\sum_{i=0}^{(nx+1)(ny+1)-1}a_{i}\nabla
u_{i}(X_{k},Y_{K})\right]$ (9)

where $u_{i}(X_{k},Y_{k})$ denotes the basis function of the solution space
and test functions space centred on the $i-th$ grid node and evaluated in the
$k-th$ particle position, $a$ is a generic grid variable, $A$ is its particle
counterpart, $N_{p}$ is the number of particles, $nx,ny$ are the cell numbers
in the $x$ and $y$ direction and $\mathrm{D}$ stands for derivative.
For the numerical discretisation of the Eulerian aspects of the method, we
choose a Cartesian grid of quadrilateral elements. Among the possible choices,
this one lends itself naturally to computationally efficient refinements
through a quad-tree approach [40]. Moreover, we opt for Lagrangian hat
functions as basis for the functional spaces of both the solution and the test
function, setting our treatment in a $\mathbb{Q}_{1}$ framework.
This approach corresponds to bilinear weighting function of standard MPM, and
is known to cause volumetric locking and spurious oscillations in the computed
quantities. The reason can be easily understood by the above equations: when
$u$ is (bi)-linear, $\nabla u$ is constant, but potentially different in
neighbouring cells. Thus, when crossing cells, particles might experience
discontinuous solicitations. This issue can be numerically tackled with many
strategies, among which adding artificial viscosity [56] or choosing higher
order basis functions. For simplicity of implementation, we choose the first,
leaving the second for a future work.
With respect to standard MPM [57, 12], the present algorithm adds a step to
take into account the flow around obstacles. Namely, we subtract the inward
component of both linear momentum and force from the grid nodes whose cell is
marked as an obstacle cell in following equations (stemming from the method
for granular materials [3]):

$(m\mathbf{v}_{i}^{o})^{\prime}=m\mathbf{v}_{i}^{o}-\dfrac{m\mathbf{v}_{i}^{o}\cdot\mathbf{n}_{i}^{o}-|m\mathbf{v}_{i}^{o}\cdot\mathbf{n}_{i}^{o}|}{2}\mathbf{n}_{i}^{o}$
(10)
$(\mathbf{f}_{i}^{o})^{\prime}=\mathbf{f}_{i}^{o}-\dfrac{\mathbf{f}_{i}^{o}\cdot\mathbf{n}_{i}^{o}-|\mathbf{f}_{i}^{o}\cdot\mathbf{n}_{i}^{o}|}{2}\mathbf{n}_{i}^{o}$
(11)

where the subscript denotes the $i-th$ grid node, the superscript $o$ stands
for obstacle and $\mathbf{n}$ is the normal unit vector, positive if directed
outward from the obstacle.
Normal unit vectors to a domain profile $\partial\Omega$ are computed as
gradients of the Signed Distance Function, defined as:

$sd(\mathbf{x})=\begin{cases}-\inf_{y\in\partial\Omega}d(\mathbf{x},\mathbf{y})\quad\mathbf{x}\in\Omega\\\
\inf_{y\in\partial\Omega}d(\mathbf{x},\mathbf{y})\quad\mathbf{x}\notin\Omega\end{cases}$
(12)

where $d(x,y)$ is the Euclidean distance. Taking the gradient yields:

$\nabla sd(\mathbf{x})=\mathbf{n}(\mathbf{x})\quad\mathbf{x}\in\partial\Omega$
(13)

An algorithm for signed distance function computation can be found in [9].
Finally, we use explicit time integration [12], whose stability is subject to
a CFL criterion:

$\Delta t<\alpha\dfrac{\Delta x_{min}}{v_{max}}$ (14)

where $\alpha$ collects constants. The range of admissible time steps is
bounded from above by the ratio between the smallest grid cell dimension
$\Delta x_{min}$ and the highest speed $v_{max}$, considering both particles
velocity and speed of sound $c_{s}=\sqrt{\gamma\tfrac{p}{\rho}}$. The
efficiency cost of enforcing small time steps is balanced by the need of
resolving fast dynamics. In addition, explicit-in-time methods sidestep the
need to solve global linear systems at each time step, thus enjoying better
scalability properties than implicit methods.

##  3 Portable MPM Implementation

This section describes the algorithm that was presented in section 2 and how
it was implemented to ensure portability on different architectures and
efficiency on GPUs. The development of the code benefited in part from
publicly available Open Source libraries:
Thrust111https://github.com/NVIDIA/cccl/tree/main/thrust,
JSON222https://github.com/nlohmann/json, and
quadgrid333https://github.com/carlodefalco/quadgrid. In detail, the algorithm
implementation is composed of the following steps:

  1. 1.

Data generation: the creation of simulations input is performed separately
from code execution. For simple geometries it is handled via ad-hoc C++
scripts. For complex geometries requiring a careful computation of the signed
distance function, an Octave script is used instead, to take advantage of
Octave built-in functions. In both cases, since objects are considered at
rest, normal unit vectors are computed only once, and boundary cells are
marked with a flag corresponding to the boundary condition. Finally, the
initial data file is saved in JSON format.

  2. 2.

MPM class object construction: since the most computationally costly
operations are the memory related ones, such as I/O and copies between host
and device, a class object for the method is constructed only once, reading
grid and particles data from the input file, data is saved on host memory and
copied to device memory. So doing, is possible to construct functor objects
which will be called on the device during the time loop, without the need of
further copies.

  3. 3.

The time advancing loop starts; inner steps, except from output operation, are
implemented as Thrust or STL for _eachs, depending on the architecture chosen
at compile time.

    1. (i)

Grid variables are reset and current grid variables are overwritten with the
ones saved during construction stage, performing a device to device copy. This
step can be performed since all history dependence is maintained by Lagrangian
particles.

    2. (ii)

Equations (6) and (8) are applied to map particles mass, momentum and force to
the grid.

    3. (iii)

Boundary conditions are enforced on grid nodes momentum and force.
Specifically, at each time step and for each boundary node marked as of
obstacle type, inward components of momentum and force are subtracted as in
equations (10) and (11).

    4. (iv)

Grid nodes momentum is advanced in time.

    5. (v)

Velocity and acceleration G2P (7) is applied, and next, velocity is updated.
Acceleration itself is not stored, but computed as the result coming from a
local G2P of the relevant nodes ratio between force and mass. Following
standard MPM [12], we store the G2P velocity in a temporary velocity vector
$\mathbf{\overline{V}}$ and don’t use it to update particles velocity;
however, we postpone its usage.

    6. (vi)

Momentum P2G (6) is performed.

    7. (vii)

Boundary conditions on the updated grid values are enforced.

    8. (viii)

Particles strain rate tensor (velocity gradient) is computed as the G2PD (9)
of the ratio between grid nodes momentum and mass, and then particles
properties are updated.

    9. (ix)

Particles position are updated using $\mathbf{\overline{V}}$. Thus, the time
advancing scheme is the same of [12]. However, particles moving is postponed,
so that, in all the above steps, basis functions and their gradient can be
computed on the fly, without the need to store the value they had at former
particles position.

    10. (x)

Particles-to-grid connectivity map is updated with new positions - this step
is computationally efficient also thanks to the choice of a Cartesian grid.

    11. (xi)

Only if the initial data options ask for $N$ intermediate outputs, a device to
host memory copy and an output from host to disk are performed at $N$ moments
equidistant in time.

    12. (xii)

Optionally, a re-ordering algorithm is applied to keep data locality, i.e. ti
keep contiguous particles in the physical simulation also close in memory. The
algorithm is based on an indexing vector built iterating over the structure
that stores information about which particles are contained in every grid
cell, such that the n-th element of the ordering vector indicates the current
sorting particle index which has to be at n-th position after re-ordering.

    13. (xiii)

A control on the CFL condition (14) is enforced, taking into account particles
velocity (momentum over mass and $\mathbf{\overline{V}}$) and the speed of
sound.

    14. (xiv)

Time step is updated, loop ends.

  4. 4.

A final device-to-host memory copy and output is made.

Algorithm 1 summarise the for loop procedure.

Algorithm 1 Time loop sketch

1:while $t<t_{f}$ do

2: Grid reset

3: Mass, momentum and force P2G  $\qquad m=\mathrm{P2G}(M),\
m\mathbf{v}=\mathrm{P2G}(M\mathbf{V}),\
\mathbf{f}=\mathrm{P2GD}(P\tfrac{M}{\rho})$

4: Boundary conditions enforcement on $m\mathbf{v},\,\mathbf{f}$ & obstacle
treatment as in (10), (11)

5: Momentum equation solution on the grid
$\qquad(m\mathbf{v})^{\prime}=m\mathbf{v}+\Delta t\mathbf{f}$

6: G2P mapping
$\qquad\overline{\mathbf{V}}=\mathrm{G2P}(\tfrac{m\mathbf{v}}{m}),\
{\mathbf{A}}=\mathrm{G2P}(\tfrac{\mathbf{f}}{m})$

7: Velocity update  $\qquad\mathbf{V}^{\prime}=\mathbf{V}+\Delta t*\mathbf{A}$

8: Momentum P2G  $\qquad m\mathbf{v}=\mathrm{P2G}(M\mathbf{V})$

9: Boundary condition enforcement on $m\mathbf{v}$

10: Particles properties update:
$\displaystyle\qquad\nabla\mathbf{V}=\mathrm{G2PD}(\tfrac{m\mathbf{v}}{m})$
$\displaystyle\qquad E^{\prime}=E+\Delta
t\dfrac{-P}{\rho}\nabla\cdot\mathbf{V}$
$\displaystyle\qquad\rho^{\prime}=\dfrac{\rho}{1+\Delta
t\nabla\cdot\mathbf{V}}$ $\displaystyle\qquad
P^{\prime}=(\gamma-1)\rho^{\prime}E^{\prime}$

11: Particles moving $\mathbf{X}^{\prime}=\mathbf{X}+\Delta
t\overline{\mathbf{V}}$

12: Optional particles re-ordering

13: CFL condition (14) check & time advance

14:end while

This algorithm is of special interest for acceleration on GPUs since it can be
written in terms of the four kernels (6), (7), (8), (9), which presents
interesting parallelisation potentialities. Grid-to-Particles operations are
intrinsically Single Instruction Multiple Data (SIMD) parallelizable, since
data is projected from grid nodes to particles, and thus particles are
independent from one another during writing stage. Particles-to-Grid kernel
presents the same potentialities, however, if they are implemented strictly
following their analytical definition, they present a data race: different
particles need access in both read and write mode to the same grid nodes.
In this work, we chose to patch the P2G and G2P kernels written in the
external quadgrid library to adapt them to GPU execution via Thrust rather
than writing ad hoc versions in the code, at the expense of some further
optimisation that could have been added by exploiting the peculiarities of our
test cases, to arrive at more general and portable conclusions. In any case,
preliminary tests have shown that applying a kernel-fusion approach [38] where
possible would not have affected the execution time order of magnitude.
A P2G implementation is then based on the following pillars:

  1. 1.

The definition of a functor for each kernel. For example, the initial P2G on
mass could be implemented as follows:444To compile CUDA or Thrust device code
decorators are necessary; for example, the functor could be modified as:
__device__ void operator () (idx_t k), and the __device__ macro could be
defined as empty during compilation for portability, in case a non-CUDA aware
compiler is used for different architectures.

⬇

1struct functor_t

2{

3 real_t hx, hy;

4 real_iter_t x, y, M, m;

5 idx_iter_t ptcl_to_grid;

6 functor_t (){};

7 void operator () (idx_t k)

8 {

9 real_t N = 0.0, xx = x[k], yy = y[k];

10 idx_t row = grid_to_row (ptcl_to_grid[k]);

11 idx_t col = grid_to_col (ptcl_to_grid[k]);

12 for (idx_t inode = 0; inode < 4; ++inode)

13 {

14 N = shape_function (xx, yy, inode, col, row, hx, hy);

15 idx_t i = global_node(inode, col, row);

16 m[i] += N * M[k];

17 }

18 }

19};

  2. 2.

The definition of the proper counting iterators (or ranges), for example,
following Thrust or Boost Libraries:

⬇

1counting_iterator <idx_t> first_p(0), last_p(num_particles);

  3. 3.

The call to the generic algorithm, e.g.:

⬇

1for_each (execution_policy, first_p, last_p, functor);

Step 3 is the one that requires care when performed with a parallel execution
policy. Indeed, it requires the handling of possible data races. Currently, we
handle them through atomic operations, namely, in the above snippet we
substitute m[i] += N * M[k] with atomicAdd ( &m(i), N*M[k]), having re-defined
m as a real_t*, and:

  * •

relying on CUDA implementation of atomicAdd for execution on the GPU,

  * •

and defining an ad-hoc atomicAdd, based on STL std::atomic_ref, for parallel
execution on the CPU.

Atomics cost is not detrimental for overall performance on the GPU, and
portability to non-CUDA backends is easily obtained, taking care of defining
the ad-hoc atomicAdd for the CPU in such a way to be compliant with CUDA’s
API. However, in Section 5 we propose an alternative algorithm which avoids
the data-race and thus the need of atomics in P2G like kernels.

##  4 Benchmarks

In this section we show the results of 2D numerical simulations of gas
dynamics obtained with the described implementation of MPM, with a focus on
supersonic gas flow over solid obstacles. A first general evaluation of
accuracy is followed by an assessment of the sensitivity of the results to the
choice of numerical parameters. We then test the efficiency and performance of
the implementation in providing parallel scalability on GPU and its
portability to CPUs. The experiments were carried out on NVIDIA A100-40 GPUs
and AMD EPYC Rome 7402 CPUs at Leonardo S.p.A.’s in-house HPC facility
_davinci-1_.

###  4.1 Test cases

We adopt non-dimensional units as in [41], and we model the fluid as a
perfect, bi-atomic gas, with specific heat ratio:

$\gamma=\dfrac{c_{p}}{c_{v}}=1.4$ (15)

In the far field away from the obstacles, we set fluid density and pressure
as:

  * •

$\rho_{\infty}=1.4$;

  * •

$p_{\infty}=1$;

where the subscript $\infty$ denotes unperturbed conditions. Thus, we obtain a
unitary speed of sound:

$c_{s,\infty}=\sqrt{\gamma\ \dfrac{p_{\infty}}{\rho_{\infty}}}=1$ (16)

and an unperturbed Mach number equal to the reference speed:

$M_{\infty}=\tfrac{v_{\infty}}{c_{s,\infty}}=v_{\infty}$ (17)

For time integration, we enforce the time step to be:

$\Delta t\leq\dfrac{1}{2}\dfrac{h_{min}}{v_{max}+c_{s,max}}$ (18)

This requirement is more conservative than necessary, and, in the limits of
the CFL condition, could be partially relaxed. However, we choose it to
guarantee that particles can’t move across more than one grid cell per time
iteration; so doing, in the framework of our algorithm, we avoid the risk that
particles pass through any boundary, and ensure that their motion will respect
boundary conditions.
In all test cases, initial particles position is sampled by a uniform random
distribution in the fluid domain; in the following, we denote with $ppc$ the
average number of fluid particles per fluid cell, with $N_{C}$ the total
number of cells, and with $N_{P}$ the total number of particles.

#### Supersonic flow past solid obstacles

First, we consider two cases of supersonic, Mach 3 flow:

  * (a)

past a cylinder of radius $0.25$ and centre $(0.6,1)$ in a $[0,4]\times[0,2]$
rectangular domain; and

  * (b)

past a rectangular-shaped step $[0.6,3]\times[0,0.2]$ in a $[0,3]\times[0,1]$
rectangular domain.

In both cases, we enforce wind tunnel–like boundary conditions, with inflow
conditions at the left vertical boundary, outflow conditions at the right
horizontal boundary, and slip conditions at the top and bottom boundaries.
Specifically, we set:

  * •

$\mathbf{v}_{inflow}=(3,0)$

  * •

$\mathbf{v}_{top}\cdot\mathbf{n}=\mathbf{v}_{bottom}\cdot\mathbf{n}=0$

where $\mathbf{n}$ denotes the unit normal vector. A particles-recycling
approach is adopted, whereby particles leaving the domain at the outflow
boundary are regenerated in a random cell at the inflow boundary with inflow
conditions, keeping the total particles number constant throughout the
simulation.

Figure 1: Mach 3 flow past a cylinder: pressure at $t=10\,s$, $N_{C}=10^{6}$,
$ppc=16$, $N_{P}=15.6\cdot 10^{6}$ , artificial viscosity as in [56].

In test case (a), as gas particles travel past the cylinder, a shock wave
develops at the front and two oblique shocks in the rear. The first shock is
reflected by walls and interacts with the back shocks; see Figure 1 and its
caption for numerical details.
In the step test (b), a frontal shock develops that is later reflected from
the top wall and successively from the step see Figure 2 and its caption for
numerical details.

Figure 2: Mach 3 flow past a step: pressure at $t=12s$, $N_{C}=326.7\cdot
10^{3}$, $ppc=16$, $N_{P}=4.4\cdot 10^{6}$, artificial viscosity as in [56].

The MPM implementation is able to correctly simulate the main flow features,
in case of both round shapes and sharp edges and corners. Computed shock-waves
angle and flow deflection angle (Figure 1) are found to be in agreement with
experimental references [18] as well as with results with other codes, e.g.,
deal.II, when the same polynomial order is used for basis
functions555https://www.dealii.org/current/doxygen/deal.II/step_69.html#Results.
The lower level of detail in comparison with some results in the literature
([14, 33, 24]) is likely due to the first order space discretisation used
here.

#### Transonic flow past an aerofoil

Next, we test the method in simulating flow past an OAT15a aerofoil at Mach
0.73, to assess the ability of the implementation in dealing with more
complex-shaped obstacles. The signed-distance function approach, in
conjunction with the implemented obstacles treatment (10),(11), naturally fits
the aerofoil profile (Figure 3). Since, rather than using an analytic
expression, the profile and its normal unit vectors are generated reading from
a file containing a set of points, the method proves to be sufficiently
general and not bounded to particular geometric shapes.

Figure 3: Signed distance function for OAT15a aerofoil at
$\alpha=-3.5^{\circ}$

Within the $[0,0.4]\times[0,0.27]$ domain, we enforce inflow conditions at the
left boundary, with

$\mathbf{v}_{\infty}=(0.73,0)$ (19)

and outflow conditions at other boundaries. To test the physical reliability
of the approach, we choose a large attach number to emphasise dynamical
effects and display velocity magnitude and streamlines (Figure 4).

Figure 4: OAT15a Mach 0.73 flow: velocity streamlines at $\alpha=25^{\circ}$,
$t=1s$, $N_{C}=247\cdot 10^{3}$, $ppc=9$, $N_{P}=2.18\cdot 10^{6}$, artificial
viscosity as in [56].

These are found to be in line with what is expected from physical
considerations: we can see airflow speeding up as it approaches the aerofoil,
with an higher velocity magnitude in the top area with respect to the bottom
one.

#### Taylor-Green vortex

In the third test, we address the Taylor-Green Vortex test case, a well
studied, fluid-only benchmark to assess both viscous and inviscid fluid
dynamics codes in the literature, usually at subsonic regime ([10, 6, 48]);
here we simulate the case in supersonic Mach 3 conditions.
In the inviscid case, the dynamics is strongly affected by problem symmetries,
leading to “impermeable box boundary conditions” [6]. These symmetries are:

  * •

within the $[0,\pi]^{2}$ square-shaped computational domain, the fluid never
crosses the lines $x=n\pi$ and $y=n\pi$; and

  * •

close to the boundaries, both the normal velocity $\mathbf{v_{\perp}}$ and the
normal derivative ${\tfrac{\partial\mathbf{v_{//}}}{\partial\mathbf{n}}}$ of
the tangential velocity $\mathbf{v_{//}}$ vanish.

We then set the initial conditions:

  * •

$v_{x}(x,y)=v_{\infty}\sin\left(\tfrac{x}{\pi}\right)\cos\left(\tfrac{y}{\pi}\right)$;

  * •

$v_{y}(x,y)=-v_{\infty}\cos\left(\tfrac{x}{\pi}\right)\sin\left(\tfrac{y}{\pi}\right)$;

  * •

$p(x,y)=p_{\infty}+\tfrac{\rho_{\infty}}{16}\left(\left(\cos\left(\tfrac{2x}{\pi}\right)+\cos\left(\tfrac{2y}{\pi}\right)\right)\cdot
3-2\right)$.

The numerical method captures the flow features with good accuracy as compared
to results in the literature (Figure 5, cf. Figure 1b in [6]), confirming the
general applicability of the MPM algorithm in simulating compressible fluid
flow.

Figure 5: Taylor-Green vortex: computed velocity at $t=2s$, $N_{C}=3.94\cdot
10^{5}$, $ppc=9$, $N_{P}=3.55\cdot 10^{6}$, artificial viscosity as in [56].

###  4.2 Sensitivity to numerical parameters

Having addressed the applicability of the method to high speed, compressible
gas dynamics, we focus on the cylinder test case and test the dependence of
the results on key numerical parameters. Specifically, we vary the number of
grid cells, $N_{C}$, and the average number of particles per cell, $ppc$, and
evaluate the effect on numerical stability and shock-capturing ability.
First, sufficiently high spatial resolution must be used to resolve shocks,
otherwise the scheme becomes too diffusive. For the cylinder test case,
testing evidence suggests that this minimum resolution amounts to around 3500
cells. Besides resolution, supersonic compressible fluid dynamics requires
further care in setting the average number of particles per cell. Indeed,
$ppc$ affects stability in a twofold fashion. First, for low $ppc$ values
(e.g.,1), widespread oscillations are observed all along the simulation, due
to almost randomly distributed empty cells. In a recent study, a $ppc$ value
of 4 was found to be optimal for convergence of gas dynamics tests at lower
speed [44], and, indeed, a $ppc=4$ solves this issue. However, for supersonic
flow past a static obstacle, we observed that using $ppc=4$ might cause the
formation of domain regions of empty or quasi-empty cells, which gets evident
as grid resolution increases, during the initial transient phase from uniform,
undisturbed conditions, as can be seen from Figure 6.

Figure 6: Mach 3 flow past a cylinder test case, various settings: from the
first to the third row, $N_{C}$ changes from $10^{5}$ to $5\cdot 10^{5}$ and
$10^{6}$; from the first to the third column, $ppc$ changes from 4, to 9 and
to 16.

Such emptying causes later formation of nonphysical cavitation-like effects
during the transient that precedes equilibrium, which alter the dynamics, and,
if particularly strong, can cause also spurious shock waves, that interact
with the physical ones. These events can be prevented by increasing the
average particles per cell number: at Mach 3, spurious shock waves become less
relevant with $ppc=9$, and undetectable with $ppc=16$.
These effects, and in particular strong disruptions, are of no interest here,
since the initial transient phase is non-physical, and derives only from the
choice to start the simulations with undisturbed inflow conditions and
particles position sampled from a uniform random distribution instead that
from a previously computed equilibrium condition, e.g. the result of an
already converged simulation. However, it could be a relevant aspect if fast
changing dynamic situations have to be simulated, for example with a moving
cylinder in a FSI framework, since they would alter its dynamics. In such
cases, a sufficient number of particles per cell should be used. By contrast,
a $ppc=4$ value can be employed when simulating equilibrium conditions:
numerical experiments show that equilibrium is reached notwithstanding
spurious instabilities in the initial transient phase (see Figure 7 for a
snapshot).

However, for linear grid basis functions, numerical oscillations are present
for any $ppc$ value. Therefore, additional artificial viscosity is essential
to compute numerical quantities, even though going high order might be
necessary to properly resolve fine structures. For example, in the cylinder
test case the standard deviation on pressure, temporally averaged over enough
different equilibrium instants, can decrease by about one half with the
addition of a proper artificial viscosity (Figure 7 shows also effects of
numerical viscosity on a single snapshot).

Figure 7: Mach 3 flow past a cylinder test case, equilibrium snapshot with
$N_{C}=10^{5}$, $ppc=4$ and no artificial viscosity (above), and
$N_{C}=10^{6}$, $ppc=16$ and artificial viscosity. Despite the very different
setting, which determines a different initial behaviour, the two simulations
foresees according pressure around the cylinder at equilibrium, and only
slight differences in the far rear.

Finally, a sufficient randomness of particles position, especially in
particles recycling, is fundamental to obtain stable and physically meaningful
simulations, in agreement with results in the literature on Particle-In-Cell
methods in plasma physics. On longer runs (i.e. above tens of physical
seconds) a symplectic time integration scheme might be useful [31].

###  4.3 Computational performance

We then focus on the computational performance and scalability of the method,
analysing the computational time required by simulations when varying the grid
resolution and $ppc$ values. It is worth noting that the implemented parallel
kernels can be particle- or grid-based, and, respectively, assign one particle
or one grid node to one GPU thread. Moreover, we implemented both mapping
kernels (P2G and G2P) as particle-based kernels, since this proved to be the
most efficient option.

First, a detailed profiling of each step of algorithm 1 reveals that the G2P
and the P2G kernels take up $85\%$ of the parallel part of the code and more
than $60\%$ of the whole computational time. This is due both to the relative
complexity of those two kernels and to the high number of times they are
called per algorithm iteration. We thus focus the scalability analysis on
these kernels.
The first finding is that both kernels scale linearly with the total number of
particles, in agreement with their implementation (Figures 8 and 9), while
pure grid kernels, such as boundary conditions enforcement, scale with the
number of nodes. For this reason, we applied a global particles reordering
step in the algorithm and studied its effects on performance. Indeed, the
ordering is done so that particles which are in the first grid cell after
moving are assigned the lowest global indices, followed by those in the second
grid cell, and so on moving along the grid. Thus, data which will be accessed
at the next time iteration will be ordered, with contiguous particles having
contiguous indexes, keeping data locality across the simulation, which is of
utmost importance on our target architecture. The reordering turns out to have
a significant impact on the G2P kernels, resulting in a further 4x speed-up
factor, while preserving the linear scaling (Figure 8).

Figure 8: Profiling of main computational kernels: G2P kernels execution time
is compared for the cylinder test case, varying $N_{C}$ and $ppc$, considering
also re-ordering effects. Figure 9: Profiling of main computational kernels:
P2G kernels execution time is compared for the cylinder test case, varying
$N_{C}$ and $ppc$, considering also re-ordering effects.

By contrast, the reordering appears not to have an appreciable effect on the
P2G kernels, with slowdowns and accelerations evident at different particle
values (Figure 9). The reason can be found in the main implementation
difference among the two kernels. P2G introduces a data race, that is
currently dealt with using CUDA’s atomicAdd, and its cost becomes more
apparent with re-ordering. Indeed, even though, to the best of the authors’
knowledge, NVIDIA’s A100 scheduler policy is not accessible, numerical
examples appear to suggest that the scheduler behaves more in a round-robin
fashion rather than in a block partitioning way. So, a particle is assigned to
a GPU thread, until there are some free, then a second particle is assigned,
and so on. In this framework, randomisation of the global order of particles
is somewhat beneficial, since it lowers the probability that particles that
are in the same cells are processed simultaneously. After reordering,
particles close in the domain are also close in memory, enhancing probability
of a memory collision, i.e. of simultaneous access to the same data. Since in
our case such access would be in read-and-write mode, it activates the atomic
operation and degrades performance.

Figure 10: Profiling of P2G kernels: effect of average particles per cell
number on version with re-ordering.

This also explains why reordering appears to not being beneficial for P2G
kernels. Since a higher number of ordered particles per cell increases the
chance of data races, the average number of particles per cell becomes an
important parameter. Indeed, grouping simulations results per $ppc$ value, we
recover linearity, and can also confirm the trend of worsening performance at
larger numbers of ordered particles in a cell. Indeed, comparing the linear
fits in Figure 10, we can see that in the $ppc=4$ and $ppc=9$ cases execution
is faster than in the global fit, and in the $ppc=16$ case it is slower.

#### Portability

Finally, we test the application of the developed implementation also on a
different architecture and compare the performance, addressing the performance
portability of the MPM code. Specifically, we execute the cylinder test case
on $2x24$ CPU cores, running the nine numerical settings described above,
without enforcing particles re-ordering, varying the number of cells and of
particles per cell and studying how this affects the main computational
kernels (Figure 11).

Figure 11: Profiling of main computational kernels: total G2P and P2G kernels
execution time is compared for the cylinder test case, varying $N_{C}$ and
$ppc$, without re-ordering, on $2x24$ CPU cores.

While portability is automatically granted, CPU performance turns out to be
worse than single-GPU performance, even though 2x24 CPU cores are employed in
parallel. This is consistent with expectations that MPM is better suited to
GPU architectures.
As for the computational kernels:

  * •

P2G kernels scale almost linearly with the number of total particles, in
agreement with the beneficial effect of randomisation in reducing the actual
execution of atomic operations;

  * •

G2P kernels maintain quasi-linear scaling only if the average particles-per-
cell number is kept fixed;

  * •

the execution time required by all the other kernels considered together is
considerably smaller than the time required by P2G and G2P ones. This confirm
the computational relevance of P2G and G2P also on different hardware.

The main difference with respect to the GPU test lies in the comparison
between P2G and G2P performance, with the former significantly slower than the
latter on the CPU. This is most probably due to our choice of writing an ad-
hoc atomicAdd function, based on std::atomic_ref, with the same signature of
the more performant CUDA’s atomicAdd, to enhance portability and guarantee
parallelism. An algorithm to directly avoid the data-race present in P2G-like
kernels, and, thus, to provide a natively portable, parallel and atomic
operation free implementation, is proposed in the next section.

##  5 Discussion and conclusion

This paper has showcased the potential of a performance portable
implementation of the Material Point Method for compressible flows. Insight
has been given on algorithm adaptations which will provide further benefits,
moreover, the simulation of standard test cases of supersonic flow has
provided results comparable with the literature and with established software
runs with the same polynomial order (first order - Lagrangian hat functions)
and number of degrees of freedom. These results show that MPM is able to
simulate Mach 3 gas dynamics, giving qualitatively good results with a
prototype developed modifying existing implementations [58, 12], and serve as
further step towards the design of a highly efficient MPM-based supersonic FSI
solver.
As known from the literature [57, 42] and as also noted in a previous
implementation [2], $\mathbb{Q}1$ MPM suffers from strong numerical
oscillation issues, especially in the vicinity of shocks. The issue is tackled
in this paper by adding numerical viscosity in the form of additional terms in
the momentum equation [56]. A more natural and accurate choice could be
obtained by using quadratic B-Spline basis functions [49, 17], a code feature
which is currently under development.
Performance tests on NVIDIA A100 GPUs have confirmed that the main
computational steps of the MPM algorithm are the Particles-to-Grid and Grid-
to-Particles kernels, and that they scale well with the size of the problem.
Moreover, adding a reordering step is beneficial to data locality and,
consequently, performance, while data races in the P2G kernel can be tackled
CUDA’s atomics at acceptable cost, even though alternative exist. Last, this
approach has shown to be easily portable to different hardware architectures,
as CPUs.

Several future developments are planned in connection with the implementation
and modelling choices made in this paper:

  * •

On the modelling side, more complex test cases could be simulated, eventually
leading to a 3D extension and Fluid-Structure Interaction models. In addition,
the effect on parallel performance of higher order B-Splines under development
will need to be carefully investigated.

  * •

For the algorithm implementation, the use of CUDA atomic operation of GPUs
choice was chosen for simplicity, also because it shown an acceptable impact
on computational time and scalability (section 4). To reach the full
parallelisation capabilities proved by the G2P kernels also for the P2G ones,
P2G algorithm adaptations will be necessary. An attractive and and more
portable route would involve avoiding atomics, and writing a grid colouring
algorithm that can be implemented in an architecture-independent way, relying
only on standard libraries functions. Its main steps are:

    1. 1.

Build grid-to-particle connectivity

      1. (i)

Sort particle IDs by cell index

      2. (ii)

Count particles per cell

      3. (iii)

Perform a cumulative sum of particle counts

    2. 2.

Assign a colour to each cell, in such a way that neighbouring cells have
different colours; then, for each colour:

      1. (a)

Perform a for_each on cells

      2. (b)

Call kernels (6) or (8)

In this way we would iterate on the vector holding the list of particles in
each cell, and then inside each cell we would perform only a small amount of
serial operations, without any data race.

  * •

On the High Performance Computing side, avenues of development include the
porting of the existing code to different backends, such as HIP and ROCm, the
support of multi-GPU and multi-node architectures via MPI, and a grid-
particles adaptive refinement through a quad-tree approach [47, 13]. Moreover,
currently the code employs Structure-of-Arrays for the memory, since,
differently from arrays of arbitrary structures, it provides coalescing, which
improves memory efficiency [27]. However, recent developments on GPUs
computing literature [54, 20] show that the best performance is achieved by
adopting an Array of SoA (AoSoA). When addressing the MPI-based multi-GPU
implementation, this possibility could be explored: indeed, the current data
structure gives a conceptually straightforward extension towards multi-GPU
implementations, where a further level of parallelisation through an MPI-based
domain decomposition will make each GPU hold a SoA as the current one,
globally giving rise to an AoSoA data structure.

## Acknowledgements

Research was sponsored by the Leonardo Labs, Leonardo S.p.A., through the PhD
scholarship “Multi-GPUs based PIC / MPM methods for Compressible Flows” funded
at Politecnico di Milano and by ICSC - Centro Nazionale di Ricerca in High
Performance Computing, Big Data, and Quantum Computing funded by European
Union - NextGenerationEU. Views and conclusions contained in this document are
those of the author and should not be interpreted as representing the official
policies, either expressed or implied, of Leonardo Labs or ICSC.

## Conflict of interest

The authors declare that they have no conflict of interest.

## References

  * [1] Andersen, S., Andersen, L.: Modelling of landslides with the material-point method.  Computational Geosciences 14(1), 137–147 (2010)
  * [2] Baioni, P.J., Benacchio, T., Capone, L., de Falco, C.: GPUs based material point method for compressible flows.  VIII International Conference on Particle-Based Methods pp. 1–12 (2023)
  * [3] Bardenhagen, S., Brackbill, J., Sulsky, D.: The material-point method for granular materials.  Computer Methods in Applied Mechanics and Engineering 187(3), 529–541 (2000)
  * [4] Bardenhagen, S.G., Kober, E.M.: The generalized interpolation material point method.  Computer Modeling in Engineering & Sciences 5(6), 477–496 (2004)
  * [5] Berzins, M.: Energy conservation and accuracy of some MPM formulations.  Computational Particle Mechanics 9(6), 1205–1217 (2022)
  * [6] Brachet, M.E., Meiron, D.I., Orszag, S.A., Nickel, B.G., Morf, R.H., Frisch, U.: Small-scale structure of the taylor–green vortex.  Journal of Fluid Mechanics 130, 411–452 (1983)
  * [7] Brackbill, J., Ruppel, H.: FLIP: A method for adaptively zoned, particle-in-cell calculations of fluid flows in two dimensions.  Journal of Computational Physics 65(2), 314–343 (1986)
  * [8] Breyer, M., Van Craen, A., Pflüger, D.: A comparison of SYCL, OpenCL, CUDA, and OpenMP for massively parallel support vector machine classification on multi-vendor hardware.  In: International Workshop on OpenCL, IWOCL’22. Association for Computing Machinery, New York, NY, USA (2022)
  * [9] Bruchon, J., Digonnet, H., Coupez, T.: Using a signed distance function for the simulation of metal forming processes: Formulation of the contact condition and mesh adaptation. from a Lagrangian approach to an Eulerian approach.  International Journal for Numerical Methods in Engineering 78(8), 980–1008 (2009)
  * [10] Bull, J.R., Jameson, A.: Simulation of the Taylor–Green vortex using high-order flux reconstruction schemes.  AIAA Journal 53(9), 2750–2761 (2015)
  * [11] Cao, Y., Chen, Y., Li, M., Yang, Y., Zhang, X., Aanjaneya, M., Jiang, C.: An efficient B-spline Lagrangian/Eulerian method for compressible flow, shock waves, and fracturing solids.  ACM Trans. Graph. 41(5), 1–13 (2022)
  * [12] Chen, Z., Brannon, R.M.: An evaluation of the material point method.  Sandia National Lab. Technical Report (2002).  URL https://www.osti.gov/biblio/793336
  * [13] Cheon, Y.J., Kim, H.G.: An adaptive material point method coupled with a phase-field fracture model for brittle materials.  International Journal for Numerical Methods in Engineering 120(8), 987–1010 (2019)
  * [14] Clayton, B., Guermond, J.L., Maier, M., Popov, B., Tovar, E.J.: Robust second-order approximation of the compressible Euler equations with an arbitrary equation of state.  Journal of Computational Physics 478, 111926 (2023)
  * [15] Coombs, W.M., Charlton, T.J., Cortis, M., Augarde, C.E.: Overcoming volumetric locking in material point methods.  Computer Methods in Applied Mechanics and Engineering 333, 1–21 (2018)
  * [16] de Vaucorbeil, A., Nguyen, V.P., Hutchinson, C.R.: A total-Lagrangian Material Point Method for solid mechanics problems involving large deformations.  Computer Methods in Applied Mechanics and Engineering 360, 112783 (2020)
  * [17] de Vaucorbeil, A., Nguyen, V.P., Sinaie, S., Wu, J.Y.: Chapter two - Material point method after 25 years: Theory, implementation, and applications.  In: S.P. Bordas, D.S. Balint (eds.) Advances in Applied Mechanics, vol. 53, pp. 185–398. Elsevier (2020)
  * [18] Dyke, M.V.: An Album of Fluid Motion.  ParabolicPress, Inc. (2008)
  * [19] Elafrou, A., Larkin, J.: No More Porting: GPU Computing with Standard C++ and Fortran.  NVIDIA - GTC Spring 2023 (2023).  URL https://www.nvidia.com/en-us/on-demand/session/gtcspring23-s51043/
  * [20] Fei, Y., Huang, Y., Gao, M.: Principles towards real-time simulation of material point method on modern GPUs.  arXiv preprint (2021).  URL https://arxiv.org/abs/2111.00699
  * [21] Fu, C., Guo, Q., Gast, T., Jiang, C., Teran, J.: A polynomial Particle-in-Cell method.  ACM Trans. Graph. 36(6), 1–12 (2017)
  * [22] Gao, M., Wang, X., Wu, K., Pradhana, A., Sifakis, E., Yuksel, C., Jiang, C.: GPU optimization of material point methods.  ACM Trans. Graph. 37(6), 1–12 (2018)
  * [23] Goktekin, T.G., Bargteil, A.W., O’Brien, J.F.: A method for animating viscoelastic fluids.  ACM Trans. Graph. 23(3), 463–468 (2004)
  * [24] Guermond, J.L., Nazarov, M., Popov, B., Tomas, I.: Second-order invariant domain preserving approximation of the Euler equations using convex limiting.  SIAM Journal on Scientific Computing 40(5), A3211–A3239 (2018)
  * [25] Hariri, F., Tran, T., Jocksch, A., Lanti, E., Progsch, J., Messmer, P., Brunner, S., Gheller, C., Villard, L.: A portable platform for accelerated PIC codes and its application to GPUs using OpenACC.  Computer Physics Communications 207, 69–82 (2016)
  * [26] Harlow F. H.; Evans, M.W., Harris Jr., D.E.: The Particle-in-Cell method for two-dimensional hydrodynamic problems.  Report LAMS-2082 of the Los Alamos Scientific Laboratory (1956)
  * [27] Hoberock, J.: Thrust by example: Advanced features and techniques.  In: GPU Technology Conference. NVIDIA (2010)
  * [28] Hu, Y., Fang, Y., Ge, Z., Qu, Z., Zhu, Y., Pradhana, A., Jiang, C.: A moving least squares material point method with displacement discontinuity and two-way rigid body coupling.  ACM Transactions on Graphics 37(4), 1–14 (2018)
  * [29] Jiang, C.: The material point method for the physics-based simulation of solids and fluid.  Ph.D. thesis, University of California, Los Angeles (2015)
  * [30] Jiang, C., Schroeder, C., Teran, J.: An angular momentum conserving affine-particle-in-cell method.  Journal of Computational Physics 338, 137–164 (2017)
  * [31] Kraus, M., Kormann, K., Morrison, P.J., Sonnendrücker, E.: GEMPIC: geometric electromagnetic particle-in-cell methods.  Journal of Plasma Physics 83(4), 905830401 (2017)
  * [32] Laskowski, G., Kopriva, J., Michelassi, V., Shankaran, S., Paliath, U., Bhaskaran, R., Wang, Q., Talnikar, C., Wang, Z., Jia, F.: Future directions of high fidelity CFD for aerothermal turbomachinery analysis and design.  46th AIAA Fluid Dynamics Conference (2016)
  * [33] Maier, M., Kronbichler, M.: Efficient parallel 3d computation of the compressible Euler equations with an invariant-domain preserving second-order finite-element scheme.  ACM Trans. Parallel Comput. 8(3), 1–30 (2021)
  * [34] Moutsanidis, G., Long, C.C., Bazilevs, Y.: IGA-MPM: The isogeometric material point method.  Computer Methods in Applied Mechanics and Engineering 372, 113346 (2020)
  * [35] NVIDIA: CUDA C++ programming guide (2022).  URL https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html
  * [36] NVIDIA: NVIDIA A100 Tensor Core GPU architecture (2022).  URL https://resources.nvidia.com/en-us-genomics-ep/ampere-architecture-white-paper?xs=169656
  * [37] NVIDIA: Hpc compilers - c++ parallel algorithm (2023).  URL https://docs.nvidia.com/hpc-sdk/compilers/c++-parallel-algorithms/index.html
  * [38] NVIDIA: Thrust release 12.3 (2023).  URL https://docs.nvidia.com/
  * [39] Sadeghirad, A., Brannon, R.M., 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)
  * [40] Samet, H.: The quadtree and related hierarchical data structures.  ACM Computing Surveys 16(2), 187–260 (1984)
  * [41] Sod, G.A.: A Survey of Several Finite Difference Methods for Systems of Nonlinear Hyperbolic Conservation Laws.  Journal of Computational Physics 27(1), 1–31 (1978)
  * [42] Steffen, M., Kirby, R.M., Berzins, M.: Analysis and reduction of quadrature errors in the material point method (MPM).  International Journal for Numerical Methods in Engineering 76(6), 922–948 (2008)
  * [43] Stomakhin, A., Schroeder, C., Jiang, C., Chai, L., Teran, J., Selle, A.: Augmented MPM for phase-change and varied materials.  ACM Transactions on Graphics 33(4), 1–11 (2014)
  * [44] Su, Y.C., Tao, J., Jiang, S., Chen, Z., Lu, J.M.: Study on the fully coupled thermodynamic fluid–structure interaction with the material point method.  Computational Particle Mechanics 7(2), 225–240 (2020)
  * [45] Sulsky, D., Chen, Z., Schreyer, H.: A particle method for history-dependent materials.  Computer Methods in Applied Mechanics and Engineering 118(1-2), 179–196 (1994)
  * [46] Sulsky, D., Gong, M.: Improving the Material-Point Method, pp. 217–240.  Springer International Publishing, Cham (2016)
  * [47] Tan, H., Nairn, J.A.: Hierarchical, adaptive, material point method for dynamic energy release rate calculations.  Computer Methods in Applied Mechanics and Engineering 191(19), 2123–2137 (2002)
  * [48] Taylor, G.I., 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)
  * [49] Tielen, R., Wobbes, E., Möller, M., Beuth, L.: A high order material point method.  Procedia Engineering 175, 265–272 (2017).  Proceedings of the 1st International Conference on the Material Point Method (MPM 2017)
  * [50] Toro, E.F.: Riemann Solvers and Numerical Methods for Fluid Dynamics.  Springer (2009)
  * [51] Tran, L.T., Kim, J., Berzins, M.: Solving time-dependent PDEs using the material point method, a case study from gas dynamics.  International Journal for Numerical Methods in Fluids 62(7), 709–732 (2010)
  * [52] Tran, Q.A., Sołowski, W.: Generalized interpolation material point method modelling of large deformation problems including strain-rate effects – application to penetration and progressive failure problems.  Computers and Geotechnics 106, 249–265 (2019)
  * [53] Wang, S., Ding, M., Gast, T.F., Zhu, L., Gagniere, S., Jiang, C., Teran, J.M.: Simulation and visualization of ductile fracture with the Material Point Method.  Proc. ACM Comput. Graph. Interact. Tech. 2(2), 1–20 (2019)
  * [54] Wang, X., Qiu, Y., Slattery, S.R., Fang, Y., Li, M., Zhu, S.C., Zhu, Y., Tang, M., Manocha, D., Jiang, C.: A massively parallel and scalable multi-GPU material point method.  ACM Trans. Graph. 39(4), 1–15 (2020)
  * [55] Wobbes, E., Möller, M., Galavi, V., Vuik, C.: Conservative Taylor least squares reconstruction with application to material point methods.  International Journal for Numerical Methods in Engineering 117(3), 271–290 (2019)
  * [56] X. Zhang Z. Chen, Y.L.: The Material Point Method.  Academic Press, Elsevier (2017)
  * [57] York, A.: Development of modifications to the material point method for the simulation of thin membranes, compressible fluids, and their interactions.  Ph.D. thesis, University of New Mexico (1997)
  * [58] York, A.R., Sulsky, D., Schreyer, H.L.: Fluid–membrane interaction based on the material point method.  International Journal for Numerical Methods in Engineering 48, 901–924 (2000)