# Stability properties of the Discontinuous Galerkin Material Point Method for hyperbolic problems in one and two space dimensions

Adrien Renaud, Thomas Heuzé, Laurent Stainier

► To cite this version:

Adrien Renaud, Thomas Heuzé, Laurent Stainier. Stability properties of the Discontinuous Galerkin Material Point Method for hyperbolic problems in one and two space dimensions. International Journal for Numerical Methods in Engineering, 2019, 10.1002/nme.6239. hal-02378909

HAL Id: hal-02378909

https://hal.science/hal-02378909

Submitted on 25 Nov 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Stability properties of the Discontinuous Galerkin Material Point Method for hyperbolic problems in one and two space dimensions

Adrien Renaud¹,², Thomas Heuzé², Laurent Stainier²

¹MSSMat Laboratory, CentraleSupélec, France

²Research Institute in Civil and Mechanical Engineering, Ecole Centrale de Nantes, France

## Abstract

In this paper, stability conditions are derived for the Discontinuous Galerkin Material Point Method [1, 2] on the scalar linear advection equation for the sake of simplicity and without loss of generality for linear problems. The discrete systems resulting from the application of the DGMPM discretization in one and two space dimensions are first written. For these problems a second-order Runge-Kutta and the forward Euler time discretizations are respectively considered. Moreover, the numerical fluxes are computed at cell faces by means of either the Donor-Cell Upwind [3] or the Corner Transport Upwind [4] methods for multi-dimensional problems. Second, the discrete scheme equations are derived assuming that all cells of a background grid contain at least one particle. Although a Cartesian grid is considered in two space dimensions, the results can be extended to regular grids. The von Neumann linear stability analysis then allows the computation of the critical Courant number for a given space discretization. Though the DGMPM is equivalent to the first-order finite volume method if one particle lies in each element[1], so that the Courant number can be set to unity, other distributions of particles may restrict the stability region of the scheme. The study of several configurations is then proposed.

**Keywords** — Discontinuous Galerkin Material Point Method, Stability analysis, Critical CFL number, Hyperbolic problems, Corner Transport Upwind method

## 1 Introduction

The numerical simulation of physical problems modeled by means of hyperbolic partial differential equations involves the solution of possibly discontinuous waves. The precise tracking of those waves is of major importance in solid mechanics, especially for history-dependent constitutive equations, as it allows to properly assess the residual states. More specifically, high-speed forming techniques such that electromagnetic forming [5] also requires the tracking of solid interfaces. Nevertheless, the accurate simulation of such problems can be prevented by several reasons. First, Lagrangian formulations of mesh-based approaches such as the widely used Finite Element (FEM) [6] and Finite Volume (FVM) [3] methods become less accurate for very large deformations due to severe distortions of the mesh. Although those instabilities can be avoided by resorting to Eulerian or Arbitrary Lagrangian-Eulerian formulations, other difficulties arise owing to additional diffusive projections of fields. Second, it is well-known that explicit finite element schemes can exhibit numerical noise near sharp solutions that can be hard to remove with artificial viscosity without loss of accuracy. Such oscillations can nevertheless be removed from FVM solutions due to numerical fluxes involved in the formulation, allowing the building of Total Variation Diminishing (TVD) schemes [7]. The Discontinuous Galerkin approximation in space [8], combined with FEM (DGFEM), enables to take advantage of similar interface fluxes so as to construct Total Variation Diminishing in the Means (TVDM) finite element procedures [9]. While the introduction of the DG approximation within FEM schemes enables to avoid non-physical oscillations, providing the use of suitable limiters [10], these approaches are constrained by a restrictive CFL condition [11] and hence, suffer from numerical diffusion. Space-time DGFEM formulations [12] enable to circumvent this drawback but are nonetheless subject to mesh entanglement.

One possibility to avoid mesh tangling instabilities while providing a material description of the motion is to use mesh-free approaches such as those of the Particle-In-Cell (PIC) family [13] and, in particular, the Material Point Method [14]. The MPM is based on a dual discretization of a domain made of a collection of material points lying in an arbitrary grid. A discrete system is solved on the grid, whereas the loading history is stored at particles during the motion so that field projections, which introduce some freedom into the scheme, are required. Indeed, the updated velocity at the grid level can be directly interpolated to the particles according to the original PIC procedure. Alternatively, the particle velocities can be updated by interpolating the nodal acceleration resulting from the solution of the discrete system, as introduced in PIC by the FLuid Implicit Particle method (FLIP) [15]. The latter allows to reduce numerical diffusion at the cost of spurious oscillations [16]. Recently, a tunable mapping procedure, based on a parameter *m* that completely eliminates the noise in MPM solutions, has been proposed in the
Extended PIC of order $m$ (XPIC(m)) [17]. A classical interpolation is selected for $m = 1$ whereas the mapping tends to FLIP one as $m \to \infty$. Nevertheless, the numerical diffusion still prevents from capturing (discontinuous) waves.

Note also that the numerical diffusion exhibited by the PIC can be limited by reducing the domain of influence of nodes rather than modifying the projections themselves. This approach is followed in the Discontinuous Galerkin Material Point Method (DGMPM) [1, 2]. The introduction of the DG approximation within the MPM, combined with the use of the PIC projection, thus aims at providing non-oscillating discontinuous solutions with low numerical diffusion due to the support of the shape functions that reduces to one cell. Therefore, the DGMPM enables to take advantage of space-DGFEM and MPM in order to accurately follow waves in finite-deforming media.

In that method, the weak form of a system of conservation laws is written on an arbitrary grid and numerical fluxes arising from the DG approximation are computed at cell faces by means of an approximate Riemann solver. Those intercell fluxes allow to take into account the characteristic structure of hyperbolic problems, and in particular the transverse propagation of waves through the use of the Corner Transport Upwind method (CTU) [4] developed for finite volumes. The CTU is however reformulated in order to fit the DGMPM approximation, in which fields are edge-wise constant in Riemann problems rather than cell-wise constant as it is the case for FVM. Furthermore, as in MPM, all the fields are stored at material points moving in the arbitrary grid, the mapping between nodes and particles being made with PIC projection so that non-oscillating solutions are provided. Hence, the geometry is updated at the particles level based on a single-valued velocity field. As a first development step, the DGMPM has been constructed upon a total Lagrangian formulation. The numerical results provided by the method for a plane wave problem in a finite-deforming hyperelastic material showed excellent agreement with the exact solution consisting of either a shock or a rarefaction wave [1]. It is worth noticing that similar results can be obtained with the FVM written in the reference configuration, the key point being however that the update of the geometry is straightforward in the DGMPM. Furthermore, an interesting feature of the method consists in allowing the employment of mesh adaption strategies so that waves can be accurately captured in the current configuration. On the other hand, only linear shape functions leading to a first-order accuracy [2] have been considered, the extension of the method to higher-order approximations being the object of ongoing works. The low-order approximation is however not seen as a shortcoming for now since the development of the method has been focused so far on capturing discontinuous solutions for which the accuracy of any numerical scheme falls to one [3].

Nonetheless, the stability of the DGMPM highly depends on the distribution of particles inside the computational grid. Indeed, the stability analysis of the one-dimensional DGMPM scheme coupled to the forward Euler time integration [1] yields a stability condition that depends on the space discretization and the CFL number. Conversely, one can ensure the stability for a given distribution of material points by finding the maximum CFL number satisfying the aforementioned relation. Such a condition, which does not exist for the MPM and other DGMPM discretizations, is crucial since it allows to: (i) ensure the stability of the scheme while minimizing the numerical diffusion; (ii) adapt the Courant number when the grid is reconstructed; (iii) adapt the grid so that a given CFL condition is satisfied. It is the purpose of this paper to provide the stability conditions for the one-dimensional DGMPM scheme combined with the two-stage second-order Runge-Kutta (RK2) time discretization and for the two-dimensional DGMPM coupled to the forward Euler algorithm. Although the solution of linear equations is considered here, the results presented must be put into the context of the problems aimed by the method, involving large deformation, and for which the DGMPM enables plenty of perspectives.

In the following, the DGMPM discrete system for the multi-dimensional scalar linear advection equation is derived and the computation of interface fluxes, as well as the solution scheme, are recalled in section 2. In particular, we shall see that the adaptation of the Corner Transport Upwind method (CTU) [4] to DGMPM leads to the same corrections of intercell fluxes as for finite volumes. Second, the system resulting from the combination of DGMPM and RK2 discretization is specialized to one-dimensional problems in section 3 so that the von Neumann linear stability analysis is carried out. At last, the same approach is followed in section 4 for the DGMPM scheme coupled to the forward Euler time integration applied to the two-dimensional problem.

# 2 Discrete system for the scalar linear advection equation

Consider a solid domain with constant volume $\Omega \in \mathbb{R}^3$, bounded by the surface $\partial\Omega$ within the time interval $t \in \tau = [0, T]$. In that domain, a material particle is located in the Cartesian coordinates system by the vector $\vec{x} = x_\alpha \vec{e}_\alpha$, following the convention of implicit summation over repeated indices.

## 2.1 The model equation

We focus here on the transport of an arbitrary scalar quantity $q$ in that domain, governed by the advection equation:

$$ \frac{\partial q}{\partial t}(\vec{x}, t) + \vec{\nabla} \cdot \vec{F}(\vec{x}, t) = 0 \quad \forall \vec{x}, t \in \Omega \times \tau \qquad (1) $$

where $\vec{\nabla} \cdot (\bullet)$ is the right divergence operator, and $\vec{F} = q(\vec{x}, t)s_\alpha \vec{e}_\alpha$, the flux vector. In the expression of the flux vector, $s_\alpha \in \mathbb{R}$ is the speed at which the quantity $q$ is advected in direction $\vec{e}_\alpha$. For the linear advection equation considered here, these
celerities are constant and equation (1) can be rewritten as:

$$ \frac{\partial q}{\partial t}(\vec{x}, t) + s_{\alpha} \frac{\partial q}{\partial x_{\alpha}}(\vec{x}, t) = 0 \quad \forall \vec{x}, t \in \Omega \times \tau $$

## 2.2 DGMPM space discretization

The continuum body $\Omega$ is discretized into a set of $N_p$ material points in a Cartesian grid made of $N_n$ nodes and $E$ non-overlapping cells of volume $\Omega^e$. The boundary of the computational domain is defined by the set of faces separating empty cells from those containing particles (see figure 1 for a two-dimensional example). Even though MPM, and in turn DGMPM, have been developed

<figure><img src="image_2.png" /><figcaption>Figure 1: Representation of a continuum body by a set of material points in a regular grid in $\mathbb{R}^2$.</figcaption></figure>

for solid dynamics based on the balance equation of linear momentum, the formulation can be extended to equation (1) by means of a fictitious mass density $\rho$. Indeed, the methods rely on the representation of the mass density that weights the time partial derivative of the hyperbolic equation by means of a delta Dirac characteristic function:

$$ \rho(\vec{x}) = \sum_{I=1}^{N_p} m_I \delta(\vec{x}^I - \vec{x}) \qquad (2) $$

where $\vec{x}^I$ is the position of the $I$th material point and $m_I$ the associated mass. In the following, a unit mass density is assumed to derive the DGMPM. Analogously to FEM and MPM, the quantity $q$ is approximated on the background grid as:

$$ q(\vec{x}, t) = \sum_{i=1}^{N_n} S_i(\vec{x}) q^i(t) \qquad (3) $$

with $q^i$ the quantity at node $i$ and $S_i(\vec{x})$ the shape function attached to that node. Note that the convention of denoting particle and nodal fields by uppercase and lowercase indices respectively is used in the remainder of the paper.

The key idea of DG methods is to allow jump of fields across mesh element faces by using broken polynomial spaces for the approximate solution [18]:

$$ \mathcal{V}^k = \{\mathcal{V} \in H^k(\Omega^e)\} ; \quad \mathcal{V}_h^k = \{\mathcal{V} \in \mathcal{P}^k(\Omega^e)\} \subset \mathcal{V}^k $$

$H^k(\Omega^e)$ being the Sobolev space, and $\mathcal{P}^k(\Omega^e)$ the space of polynomials of degree $k$ in $\Omega^e$. We restrict our attention here to linear polynomials ($k=1$). Multiplying equation (1) by a test function $\mathcal{V}$ yields the element-wise weak formulation of the problem:

$$ \text{Find } q \in \mathcal{V}_h^1 \text{ such that } \forall e \qquad \int_{\Omega^e} \mathcal{V} \frac{\partial q}{\partial t} d\Omega + \int_{\partial \Omega^e} \mathcal{V} \vec{\nabla} \cdot \vec{F} d\Omega = 0 \quad \forall \mathcal{V}, t \in \mathcal{V}_h^1 \times \tau $$

Then, the use of Gauss' divergence theorem leads to:

$$ \text{Find } q \in \mathcal{V}_h^1 \text{ such that } \forall e \qquad \int_{\Omega^e} \left[ \frac{\partial q}{\partial t} \mathcal{V} - \vec{F} \cdot \vec{\nabla} \mathcal{V} \right] d\Omega + \int_{\partial \Omega^e} \left( \vec{F} \cdot \vec{n} \right) \mathcal{V} d\Gamma = 0 \quad \forall \mathcal{V}, t \in \mathcal{V}_h^1 \times \tau $$
where $\partial\Omega^e$ is the boundary of the $\text{th}$ element with outward normal vector $\vec{n}$, and $\vec{\nabla}(\bullet)$ is the gradient operator. The inner product $\vec{F} \cdot \vec{n}$, corresponding to intercell flux, is written $F_n$ for simplicity. Next, the introduction of specific fields:

$$q = \rho \vec{q} \quad ; \quad \vec{F} = \rho \vec{F} = \rho \vec{s} \vec{q}$$

combined with the definition of mass density (2), leads to the following weak form:

$$\sum_{I=1}^{N_p} m_I \left[ \frac{\partial \vec{q}}{\partial t} \nabla - \vec{s} \cdot \vec{q} \cdot \vec{\nabla} \nabla \right]_{|\vec{x}=\vec{x}^I} + \int_{\partial\Omega^e} F_n \nabla d\Gamma = 0 \quad \forall \nabla, t \in \mathcal{Y}_h^1 \times \tau \quad (4)$$

Introduction of the DGMPM approximation (3) and arbitrariness of the test field in the weak form (4) finally provide the semi-discrete system that must be solved on the grid:

$$\sum_{I=1}^{N_p} \left[ S_{iI} m_I S_{jI} \frac{\partial \vec{q}^j}{\partial t} - \frac{\partial S_{iI}}{\partial x_\alpha} m_I S_{jI} s_\alpha \vec{q}^j \right] + \int_{\Gamma_e} S_i(\vec{x}) F_n d\Gamma = 0 \quad \forall t \in \tau \quad (5)$$

or, in matrix form:

$$M_{ij} \frac{\partial \vec{q}^j}{\partial t} - K_{ij}^{\alpha} s_{\alpha} \vec{q}^j + \hat{F}^i = 0$$

**Remark 1** As for MPM, the particle-based quadrature rule yields a consistent mass matrix $M_{ij}$ that may be singular due to reduced integration. This can be circumvented by resorting to the diagonally lumped mass matrix [19] $M_i^L = \sum_{j=1}^{N_n} M_{ij}$. The aforementioned numerical trick, which is also employed in explicit FEM [6] for efficiency purposes, can be avoided within the MPM owing to recent works. Indeed, the use of moving least squares [20] or spline interpolation [21, 22] techniques enables the evaluation of a reconstructed function at Gauss' point locations. In these versions of the MPM, the particles are no longer seen as integration points. Similar reconstructions can be considered within the DGMPM, which would amount to provide an arbitrary background grid to a DG scheme. This is, however, out of the scope of the present work.

**Remark 2** The DGMPM discretization can be constructed upon regular grids that allow the straightforward building of an orthonormal approximation basis. In that case, the mass matrix is naturally diagonal provided the derivation of a modal formulation as can be done for DGFEM [23].

The DGMPM discrete system is finally derived by dividing the time interval $\tau$ into $N_t$ subintervals and using the explicit forward Euler method:

$$M_i^L \frac{\vec{q}^{i,k+1} - \vec{q}^{i,k}}{\Delta t^k} = K_{ij}^{\alpha} s_{\alpha} \vec{q}^{j,k} - \hat{F}^{i,k} \quad (6)$$

where the superscripts $(\bullet)^{i,k}$ denote a field evaluated at node $i$ and time step $k$. Alternatively, a second-order Runge-Kutta (RK2) explicit time discretization may be employed, leading to the following two-stage discrete form:

$$M_i^L \frac{\vec{q}^{i,k+1/2} - \vec{q}^{i,k}}{\Delta t^k} = \frac{1}{2} (K_{ij}^{\alpha} s_{\alpha} \vec{q}^{j,k} - \hat{F}^{i,k}) \\ M_i^L \frac{\vec{q}^{i,k+1} - \vec{q}^{i,k}}{\Delta t^k} = K_{ij}^{\alpha} s_{\alpha} \vec{q}^{j,k+1/2} - \hat{F}^{i,k+1/2} \quad (7)$$

## 2.3 Computation of intercell fluxes

The relaxation of the continuity of fields across cell interfaces introduced by the DG approximation allows to define Riemann problems at element faces in the direction $x_n = \vec{x} \cdot n$:

$$\begin{gathered} \frac{\partial q}{\partial t} + \frac{\partial F_n}{\partial x_n} = 0 \\ q(x_n, 0) = \begin{cases} q^{-} & \text{if } x_n < 0 \\ q^{+} & \text{if } x_n > 0 \end{cases} \end{gathered} \quad (8)$$

where the $q^{\pm}$ are downwind and upwind states. The use of DG approximation can be seen as a duplication of nodes so that those states can be obtained by averaging nodal fields on each side of the interface as depicted in figure 2 for a two-dimensional case. By doing so, only one Riemann problem per interface rather than nodal ones is considered, thus avoiding a dramatic increase in computational time. As for the original formulation of the DGMPM [1], the numerical flux at a given interface is based on the
<figure><img src="image_3.png" /><figcaption>Figure 2: Duplication of nodes at an interface and building of initial conditions of the Riemann problem (2D).</figcaption></figure>

stationary solution $q^*$ of the Riemann problem (8) according to Godunov's method [24]. For the linear advection, the stationary solution is simply that of the upwind side, namely:

$$ \begin{cases} q^* = q^{-} & \text{if } s_n > 0 \\ q^* = q^{+} & \text{if } s_n < 0 \end{cases} \qquad (9) $$

where $s_n = \vec{s} \cdot \vec{n}$, so that the intercell flux at a given interface reads $F_n = s_n q^*$.

The method derived above for the computation of normal fluxes can be seen as the Donor-Cell Upwind method (DCU) [3] in which only the influence of upwind neighbor cells is considered. For multi-dimensional problems, however, waves can travel in several directions so that contributions coming from corner cells must be taken into account in order to improve accuracy and stability of the numerical scheme. The Corner Transport Upwind method (CTU) [4] consists in considering contributions propagating in bias and coming from upwind cells sharing only a node (in two dimensions) with another. This approach has been developed for finite volumes in which fields are cell-wise constant and allows to improve the CFL condition. The approach is now reformulated for DGMPM, based on edge-wise constant fields within Riemann problems.

At each cell interface, one defines left-going and right-going fluctuations as:

$$ A^{-}(\Delta q) = F_{n}(q^{*}) - F_{n}(q^{-}) \quad ; \quad A^{+}(\Delta q) = F_{n}(q^{+}) - F_{n}(q^{*}) \qquad (10) $$

From equation (10) and the definition of the stationary state $q^*$ (9), one of the fluctuation obviously vanishes for the scalar linear advection equation. The non-zero fluctuation, on the other hand, characterizes a wave that carries a jump discontinuity $\Delta q$ at speed $s_n$ in the direction $\vec{n}$. This jump influences the state vector $q$ in the neighbor cell so that the initial data in the Riemann problems formulated on adjacent edges are changed. To illustrate the above discussion, let's consider the patch of two-dimensional regular grid cells of length $\Delta x$ shown in figure 3. Assuming that the speed vector $\vec{s}$ is such that $s_1, s_2 > 0$, the

<figure><img src="image_4.png" /><figcaption>Figure 3: Illustration of the influence of fluctuations on fluxes computed at adjacent cells interfaces.</figcaption></figure>

solution of the Riemann problem on edge (i) yields the fluctuation $A^+(\Delta q^{(i)})$ which travels the distance $s_1 \Delta t$ in cell C. On the other hand, the flux through edge (j) reads [25]:

$$ F_{n}^{(j)} = \frac{1}{\Delta t} \int_{0}^{\Delta t} s_{2} \tilde{q} dt \qquad (11) $$
where $\tilde{q}$ is the upwind state at edge (j). When using DCU, this constant state is built from nodal values in cell C and is written $q^{-(j)}$. Alternatively, $\tilde{q}$ is taken in CTU as the average of states $q^{-(i)}$ and $q^{-(j)}$, weighted by the portion of edge (j) along which they respectively act at time t, namely:

$$ \tilde{q} = \frac{s_1 t q^{-(i)} + (\Delta x - s_1 t) q^{-(j)}}{\Delta x} $$

Thus, integral (11) yields:

$$ F_n^{(j)} = s_2 \frac{s_1 \Delta t}{2 \Delta x} q^{-(i)} + s_2 \frac{2 \Delta x - s_1 \Delta t}{2 \Delta x} q^{-(j)} $$

which can be rewritten:

$$ F_n^{(j)} = s_2 q^{-(j)} - \frac{1}{2} \frac{\Delta t}{\Delta x} s_2 s_1 (q^{-(j)} - q^{-(i)}) = F_n^{(j)}(q^{-(j)}) - \frac{1}{2} \frac{\Delta t}{\Delta x} \mathcal{B}^+ \mathcal{A}^+ (\Delta q^{(i)}) \quad (12) $$

where $F_n^{(j)}(q^{-(j)})$ is the flux through the edge (j) resulting from the DCU, and $\mathcal{B}^+\mathcal{A}^+(\Delta q^{(i)})$ is a transverse correction coming from edge (i). It then comes out that the CTU leads to the same corrections in DGMPM than it does in FVM [3].

## 2.4 Solution scheme

Suppose that the quantity $q$ is known at every material points that discretize a solid domain $\Omega$ at a time increment $t^k$. Since the volume is assumed constant, the lumped mass and pseudo-stiffness matrices $M_i^L$ and $K_{ij}^\alpha$ can be computed once and for all at the beginning of the computation. Then, the DGMPM procedure followed to update the field between $t^k$ and $t^{k+1}$ consists in the following steps:

(a) project the field onto the grid by solving for the $\tilde{q}^i$ the conservation of volume quantities [1]:

$$ M_i^L \tilde{q}^i = \sum_{I=1}^{N_p} S_{iI} m_I \tilde{q}^I \quad (13) $$

(b) compute intercell fluxes by using either DCU or CTU approach, as well as volume fluxes.

(c) advance the solution in time by solving the discrete system (6) or (7)

(d) interpolate the updated fields $\tilde{q}^i$ to material points according to equation (3)

In what follows, the scheme equations associated with the above numerical procedure are derived for one and two-dimensional problems. The von-Neumann stability analysis of these equations are then carried out in order to determine an upper bound of the Courant number that ensures the stability of the method. The aforementioned upper bound is hereinafter referred to as the critical Courant number and is denoted by CFL so that the Courant number must satisfy:

$$ \max_{\alpha} (s_{\alpha}) \frac{\Delta t}{\Delta x} \leq \text{CFL} $$

## 3 Stability properties of the one-dimensional schemes

The discrete system (7) is now specialized to a one-dimensional problem for which $s_1 > 0$. Thus, a domain of length $l$ is divided with $N_p$ material points arbitrarily distributed in $E$ two-node elements of constant length $\Delta x$ (figure 4). The grid mesh is such that at least one particle lies in every cell during the computation in order to ensure that there is no hole. Moreover, periodic boundary conditions are considered to simplify the analysis. In that mesh, the cell containing the particle $I$ is denoted by $c(I)$ so

<figure><img src="image_5.png" /><figcaption>Figure 4: One-dimensional mesh made of $E$ elements of constant length $\Delta x = \frac{l}{E}$.</figcaption></figure>

that the nodes interacting with this particle are $2c(I) - 1$ and $2c(I)$. Since only scalar quantities are considered here, subscript can be used to denote nodal or particle values without ambiguity. Therefore, the linear shape functions defined in element $c(I)$ are:

$$ S_{2c(I)-1}(x) = \frac{x_{2c(I)} - x}{\Delta x} ; S_{2c(I)}(x) = \frac{x - x_{2c(I)-1}}{\Delta x} \quad x \in [x_{2c(I)-1}, x_{2c(I)}] $$

and $S_{iI}$ or $S_{i,I}$ correspond to the shape function of node $i$ evaluated at the position of the $I$th material point. In order to better distinguish nodal and particle fields, uppercase symbols are used for material point quantities.
## 3.1 The scheme equation

The method followed to write the scheme equation is to trace backward the numerical procedure described in section 2 to get an expression of the form:

$$ \bar{Q}_I^{k+1} = H(\bar{q}_J^k) \quad J = 1, \dots, N_p $$

where H stands for the DGMPM discrete solution operator.

The two stages of the RK2 algorithm can be written as:

$$ \bar{q}_i^{k+\frac{p+1}{2}} = \bar{q}_i^k + \frac{p+1}{2} \frac{\Delta t}{M_i^L} s_1 \left( \sum_{j=1}^{2E} K_{i,j} \bar{q}_j^{k+\frac{p}{2}} - \hat{F}_i^{k+\frac{p}{2}} n_i \right), \quad \text{no sum on } i \quad (14) $$

in which $\hat{F}_i^{k+\frac{p+1}{2}}$ and $n_i = \pm 1$ are respectively the intercell flux and the outward unit normal at node i, and $p = \{0, 1\}$ refers to the two stages of the scheme. Moreover, the mass density is approximated on the grid as:

$$ \rho(x) = \frac{M_{2c-1}^L + M_{2c}^L}{\Delta x} = \sum_{J=1}^{N_p^c} m_J, \quad x \in [x_{2c-1}, x_{2c}] $$

where $N_p^c$ is the number of particles in cell c. The mass is uniformly distributed between particles so that the previous definition reduces to $\rho = N_p^c m^c / \Delta x$, with $m^c$ the mass carried by particles lying in c.

First, quantities at time $t^{k+1}$ at particles are obtained by interpolating nodal solutions of the discrete system (14):

$$ \bar{Q}_I^{k+1} = S_{2c(I)-1,J} \bar{q}_{2c(I)-1}^{k+1} + S_{2c(I),J} \bar{q}_{2c(I)}^{k+1} \quad (15) $$

Second, provided linear shape functions, the lumped mass and the pseudo-stiffness matrices are:

$$ \begin{align*} M_i^L &= \sum_{J=1}^{N_p} S_{iJ} m_J = m^{c(i)} \sum_{J=1}^{N_p} S_{iJ} \\ K_{2c(I)-1,j} &= \sum_{J=1}^{N_p} \frac{\partial S_{2c(I)-1,J}}{\partial x} m_J S_{jJ} = -m^{c(i)} \sum_{J=1}^{N_p} \frac{S_{jJ}}{\Delta x} \\ K_{2c(I),j} &= \sum_{J=1}^{N_p} \frac{\partial S_{2c(I),J}}{\partial x} m_J S_{jJ} = m^{c(i)} \sum_{J=1}^{N_p} \frac{S_{jJ}}{\Delta x} \end{align*} $$

The discontinuous approximation basis moreover yields a bloc diagonal pseudo-stiffness matrix so that one can write:

$$ K_{ij} \bar{q}_j^k = K_{i,2c(i)-1} \bar{q}_{2c(i)-1}^k + K_{i,2c(i)} \bar{q}_{2c(i)}^k $$

Third, a right-going wave leads to a stationary solution of Riemann problems equals to the state of the upwind node of an interface, that is:

$$ \begin{align*} q_{2c(I)-1}^{c(I)} &= \rho \bar{q}_{2c(I)-2}^{k} = N_p^{c(I)} \frac{m^{c(I)}}{\Delta x} \bar{q}_{2c(I)-2}^{k} \\ q_{2c(I)}^{c(I)} &= \rho \bar{q}_{2c(I)}^{k} = N_p^{c(I)} \frac{m^{c(I)}}{\Delta x} \bar{q}_{2c(I)}^{k} \end{align*} $$

Therefore, gathering all the previous considerations, equation (14) reads for each node of cell c(I):

$$ \begin{align} \bar{q}_{2c(I)-1}^{k+\frac{p+1}{2}} &= \bar{q}_{2c(I)-1}^{k} - \frac{p+1}{2} \frac{s_1 \Delta t}{\Delta x} \left( \frac{f_{c(I)}^{k+\frac{p}{2}} - N_p^{c(I)} \bar{q}_{2c(I)-2}^{k+\frac{p}{2}}}{\sum_{K=1}^{N_p^{c(I)}} S_{2c(I)-1,K}} \right), \quad p = \{0,1\} \tag{16} \\ \bar{q}_{2c(I)}^{k+\frac{p+1}{2}} &= \bar{q}_{2c(I)}^{k} + \frac{p+1}{2} \frac{s_1 \Delta t}{\Delta x} \left( \frac{f_{c(I)}^{k+\frac{p}{2}} - N_p^{c(I)} \bar{q}_{2c(I)}^{k+\frac{p}{2}}}{\sum_{K=1}^{N_p^{c(I)}} S_{2c(I),K}} \right) \end{align} $$

with the volume fluxes contributions $f_c^k = \sum_{J=1}^{N_p^c} [S_{2c-1,J} \bar{q}_{2c-1}^k + S_{2c,J} \bar{q}_{2c}^k]$. Note that equation (16) involves the Courant number $s_1 \Delta t / \Delta x$.
The first stage of equation (16) (i.e. $p=0$) yields the intermediate nodal fields in cell $c(I)$:

$$ 
\begin{aligned}
\bar{q}_{2c(I)-1}^{k+\frac{1}{2}} &= \bar{q}_{2c(I)-1}^{k} - \frac{s_1 \Delta t}{2\Delta x} \left( \frac{f_{c(I)}^{k} - N_{p}^{c(I)} \bar{q}_{2c(I)-2}^{k}}{\sum_{K=1}^{N_{p}^{c(I)}} S_{2c(I)-1,K}} \right) \\
\bar{q}_{2c(I)}^{k+\frac{1}{2}} &= \bar{q}_{2c(I)}^{k} + \frac{s_1 \Delta t}{2\Delta x} \left( \frac{f_{c(I)}^{k} - N_{p}^{c(I)} \bar{q}_{2c(I)}^{k}}{\sum_{K=1}^{N_{p}^{c(I)}} S_{2c(I),K}} \right)
\end{aligned}
\quad (17) $$

and the second one (i.e. $p=1$) leads to the expression of nodal quantities at the end of the time step:

$$ 
\begin{aligned}
\bar{q}_{2c(I)-1}^{k+1} &= \bar{q}_{2c(I)-1}^{k} - \frac{s_1 \Delta t}{\Delta x} \left( \frac{\sum_{J=1}^{N_p^{c(I)}} \left[ S_{2c(I)-1,J} \bar{q}_{2c(I)-1}^{k+\frac{1}{2}} + S_{2c(I),J} \bar{q}_{2c(I)}^{k+\frac{1}{2}} \right] - N_p^{c(I)} \bar{q}_{2c(I)-2}^{k+\frac{1}{2}}}{\sum_{K=1}^{N_p^{c(I)}} S_{2c(I)-1,K}} \right) \\
\bar{q}_{2c(I)}^{k+1} &= \bar{q}_{2c(I)}^{k} + \frac{s_1 \Delta t}{\Delta x} \left( \frac{\sum_{J=1}^{N_p^{c(I)}} \left[ S_{2c(I)-1,J} \bar{q}_{2c(I)-1}^{k+\frac{1}{2}} + S_{2c(I),J} \bar{q}_{2c(I)}^{k+\frac{1}{2}} \right] - N_p^{c(I)} \bar{q}_{2c(I)}^{k+\frac{1}{2}}}{\sum_{K=1}^{N_p^{c(I)}} S_{2c(I),K}} \right)
\end{aligned}
\quad (18) $$

Then, introduction of equations (18) in the interpolation from nodes to particles (15) leads to the solution at material point I and time step $k+1$:

$$ 
\begin{aligned}
\bar{Q}_I^{k+1} &= S_{2c(I)-1,I} \bar{q}_{2c(I)-1}^{k} + S_{2c(I),I} \bar{q}_{2c(I)}^{k} + N_p^{c(I)} \frac{s_1 \Delta t}{\Delta x} \frac{S_{2c(I)-1,I}}{\sum_K S_{2c(I)-1,K}} \bar{q}_{2c(I)-2}^{k+\frac{1}{2}} \\
&\quad - \left( \frac{s_1 \Delta t}{\Delta x} \left[ S_{2c(I)-1,I} - S_{2c(I),I} \frac{\sum_J S_{2c(I)-1,J}}{\sum_K S_{2c(I),K}} \right] \right) \bar{q}_{2c(I)-1}^{k+\frac{1}{2}} \\
&\quad + \frac{s_1 \Delta t}{\Delta x} \left[ S_{2c(I),I} - S_{2c(I)-1,I} \frac{\sum_J S_{2c(I),J}}{\sum_K S_{2c(I)-1,K}} - N_p^{c(I)} \frac{S_{2c(I),I}}{\sum_K S_{2c(I),K}} \right] \bar{q}_{2c(I)}^{k+\frac{1}{2}}
\end{aligned}
$$

Nodal values $\bar{q}_i^{k+\frac{1}{2}}$ are provided by the first stage of RK2 algorithm and can be substituted the above equation:

$$ 
\begin{aligned}
\bar{Q}_I^{k+1} = S_{2c(I)-1,I} \bar{q}_{2c(I)-1}^{k} + S_{2c(I),I} \bar{q}_{2c(I)}^{k} + N_p^{c(I)} \frac{s_1 \Delta t}{\Delta x} \frac{S_{2c(I)-1,I}}{\sum_K S_{2c(I)-1,K}} \left( \bar{q}_{2c(I)-2}^{k} + \frac{s_1 \Delta t}{2\Delta x} \left( \frac{f_{c(I)-1}^{k} - N_p^{c(I)-1} \bar{q}_{2c(I)-2}^{k}}{\sum_K S_{2c(I)-2,K}} \right) \right) \\
- \frac{s_1 \Delta t}{\Delta x} \left[ S_{2c(I)-1,I} - S_{2c(I),I} \frac{\sum_J S_{2c(I)-1,J}}{\sum_K S_{2c(I),K}} \right] \left( \bar{q}_{2c(I)-1}^{k} - \frac{s_1 \Delta t}{2\Delta x} \left( \frac{f_{c(I)}^{k} - N_p^{c(I)} \bar{q}_{2c(I)-2}^{k}}{\sum_K S_{2c(I)-1,K}} \right) \right) \\
+ \frac{s_1 \Delta t}{\Delta x} \left[ S_{2c(I),I} \left( 1 - \frac{N_p^{c(I)}}{\sum_K S_{2c(I),K}} \right) - S_{2c(I)-1,I} \frac{\sum_J S_{2c(I),J}}{\sum_K S_{2c(I)-1,K}} \right] \left( \bar{q}_{2c(I)}^{k} + \frac{s_1 \Delta t}{2\Delta x} \left( \frac{f_{c(I)}^{k} - N_p^{c(I)} \bar{q}_{2c(I)}^{k}}{\sum_K S_{2c(I),K}} \right) \right)
\end{aligned}
\quad (19) $$

Note that the solution of the downstream node of the adjacent cell $\bar{q}_{2c(I)-2}^{k+\frac{1}{2}}$ results from the second equation of the set (17). This gives rise to the number of particles in the cell $c(I)-1$ which is involved in the mass density and the solution of the discrete system at the time step $k+1/2$ in that element. By rearranging formula (19), one gets:

$$ 
\begin{aligned}
\bar{Q}_I^{k+1} = & \left( S_{2c(I)-1,I} - \frac{s_1 \Delta t}{\Delta x} \left[ S_{2c(I)-1,I} - S_{2c(I),I} \frac{\sum_J S_{2c(I)-1,J}}{\sum_K S_{2c(I),K}} \right] \right) \bar{q}_{2c(I)-1}^{k} \\
& + \left( S_{2c(I),I} + \frac{s_1 \Delta t}{\Delta x} \left[ S_{2c(I),I} \left( 1 - \frac{N_p^{c(I)}}{\sum_K S_{2c(I),K}} \right) - S_{2c(I)-1,I} \frac{\sum_J S_{2c(I),J}}{\sum_K S_{2c(I)-1,K}} \right] \right) \bar{q}_{2c(I)}^{k} \\
& + \frac{1}{2} \left( \frac{s_1 \Delta t}{\Delta x} \right)^2 \left( N_p^{c(I)} \left[ \frac{S_{2c(I)-1,I}}{\sum_K S_{2c(I)-1,K}} - \frac{S_{2c(I),I}}{\sum_K S_{2c(I),K}} \right] + S_{2c(I),I} \left( \frac{N_p^{c(I)}}{\sum_K S_{2c(I),K}} \right)^2 \right) \bar{q}_{2c(I)}^{k} \\
& + N_p^{c(I)} \frac{s_1 \Delta t}{\Delta x} \left[ \frac{S_{2c(I)-1,I}}{\sum_K S_{2c(I)-1,K}} \left( 1 - \frac{s_1 \Delta t}{2\Delta x} \left( 1 + \frac{N_p^{c(I)-1}}{\sum_K S_{2c(I)-2,K}} \right) \right) + \frac{s_1 \Delta t}{2\Delta x} \frac{S_{2c(I),I}}{\sum_K S_{2c(I),K}} \right] \bar{q}_{2c(I)-2}^{k} \\
& - \frac{1}{2} \left( \frac{s_1 \Delta t}{\Delta x} \right)^2 N_p^{c(I)} \frac{S_{2c(I),I}}{\left( \sum_K S_{2c(I),K} \right)^2} f_{c(I)}^{k} + \frac{1}{2} \left( \frac{s_1 \Delta t}{\Delta x} \right)^2 N_p^{c(I)} \frac{S_{2c(I)-1,I}}{\sum_K S_{2c(I)-1,K}} \frac{f_{c(I)-1}^{k}}{\sum_K S_{2c(I)-2,K}}
\end{aligned}
\quad (20) $$
Next, the nodal solutions at time step k in equation (20) result from the projection between particles and the grid (13):

$$ \bar{q}_i^k = \frac{\sum_L S_{iL} m_L \bar{Q}_L^k}{\sum_K S_{iK} m_K} = \frac{\sum_L S_{iL} \bar{Q}_L^k}{\sum_K S_{iK}} \qquad (21) $$

In particular, volume fluxes contributions can be written:

$$ \begin{aligned} f_c^k &= \sum_{J=1}^{N_p} \left[ S_{2c-1,J} \frac{\sum_L S_{2c-1,L} \bar{Q}_L^k}{\sum_K S_{2c-1,K}} + S_{2c,J} \frac{\sum_L S_{2c,L} \bar{Q}_L^k}{\sum_K S_{2c,K}} \right] \\ f_c^k &= \sum_{L=1}^{N_p} (S_{2c-1,L} + S_{2c,L}) \bar{Q}_L^k \\ f_c^k &= \sum_{L \in c} \bar{Q}_L^k \end{aligned} \qquad (22) $$

where the partition of unity yields the last equality.

The use of mapping equations (21) and (22) allows to write after some simplifications:

$$ \begin{aligned} \bar{Q}_I^{k+1} = & \sum_L \bar{Q}_L^k \left\{ \sum_i S_{iL} \frac{S_{iI}}{\sum_K S_{iK}} + \frac{s_1 \Delta t}{\Delta X} \left[ \frac{S_{2c(L),I}}{\sum_K S_{2c(L),K}} - \frac{S_{2c(L)-1,I}}{\sum_K S_{2c(L)-1,K}} \right] \right. \\ & + \frac{s_1 \Delta t}{\Delta X} N_p^{c(I)} \left[ \frac{S_{2c(I)-1,I}}{\sum_K S_{2c(I)-1,K}} \frac{S_{2c(I)-2,L}}{\sum_K S_{2c(I)-2,K}} - \frac{S_{2c(I),I} S_{2c(I),L}}{(\sum_K S_{2c(I),K})^2} \right] \\ & + \frac{1}{2} \left( \frac{s_1 \Delta t}{\Delta X} \right)^2 N_p^{c(I)} \left( \frac{S_{2c(I),L}}{\sum_K S_{2c(I),K}} - \frac{S_{2c(I)-2,L}}{\sum_K S_{2c(I)-2,K}} \right) \left[ \frac{S_{2c(I)-1,I}}{\sum_K S_{2c(I)-1,K}} - \frac{S_{2c(I),I}}{\sum_K S_{2c(I),K}} \right] \\ & + \frac{1}{2} \left( \frac{s_1 \Delta t}{\Delta X} \right)^2 N_p^{c(I)} \frac{S_{2c(L),I}}{(\sum_K S_{2c(I),K})^2} \left[ N_p^{c(I)} \frac{S_{2c(I),L}}{\sum_K S_{2c(I),K}} - 1 \right] \\ & \left. + \frac{1}{2} \left( \frac{s_1 \Delta t}{\Delta X} \right)^2 \frac{S_{2c(L)+1,I} N_p^{c(I)}}{\sum_K S_{2c(I)-1,K} \sum_K S_{2c(I)-2,K}} \left[ 1 - N_p^{c(I)-1} \frac{S_{2c(I)-2,L}}{\sum_K S_{2c(I)-2,K}} \right] \right\} \end{aligned} \qquad (23) $$

The three first terms of the latter scheme equation correspond to that obtained for the Euler algorithm [1], whereas the second-order terms are provided by the two-stage time integration. Alternatively and given the compact support of shape functions, equation (23) can be split into two sums, one over the material points contained in the cell c(I), and one over these of cell c(I) − 1:

$$ \begin{aligned} \bar{Q}_I^{k+1} = & \sum_{L \in c(I)} \bar{Q}_L^k \left\{ \sum_i S_{iL} \frac{S_{iI}}{\sum_K S_{iK}} + \frac{s_1 \Delta t}{\Delta X} \left[ \frac{S_{2c(L),I}}{\sum_K S_{2c(L),K}} - \frac{S_{2c(L)-1,I}}{\sum_K S_{2c(L)-1,K}} - N_p^{c(I)} \frac{S_{2c(I),I} S_{2c(I),L}}{(\sum_K S_{2c(I),K})^2} \right] \right. \\ & + \frac{1}{2} \left( \frac{s_1 \Delta t}{\Delta X} \right)^2 N_p^{c(I)} \left( \frac{S_{2c(I),L}}{\sum_K S_{2c(I),K}} \left[ \frac{S_{2c(I)-1,I}}{\sum_K S_{2c(I)-1,K}} - \frac{S_{2c(I),I}}{\sum_K S_{2c(I),K}} \right] + \frac{S_{2c(L),I}}{(\sum_K S_{2c(L),K})^2} \left[ N_p^{c(I)} \frac{S_{2c(I),L}}{\sum_K S_{2c(I),K}} - 1 \right] \right) \\ & + \sum_{L \in c(I)-1} \bar{Q}_L^k N_p^{c(I)} \left\{ \frac{s_1 \Delta t}{\Delta X} \frac{S_{2c(I)-1,I}}{\sum_K S_{2c(I)-1,K}} \frac{S_{2c(L),L}}{\sum_K S_{2c(L),K}} \right. \\ & \left. + \frac{1}{2} \left( \frac{s_1 \Delta t}{\Delta X} \right)^2 \left( \frac{S_{2c(L)+1,I}}{\sum_K S_{2c(I)-1,K}} \left[ \frac{1}{\sum_K S_{2c(L),K}} - \frac{S_{2c(L),L} N_p^{c(I)-1}}{(\sum_K S_{2c(L),K})^2} \right] - \frac{S_{2c(L),L}}{\sum_K S_{2c(L),K}} \left[ \frac{S_{2c(I)-1,I}}{\sum_K S_{2c(I)-1,K}} - \frac{S_{2c(I),I}}{\sum_K S_{2c(I),K}} \right] \right) \right\} \end{aligned} $$

or, for simplicity:

$$ \bar{Q}_I^{k+1} = \sum_{L=1}^{N_p} H_{IL} \bar{Q}_L^k \qquad (24) $$

When only one point lies in each cell of the grid, the sum of shape functions over particles reduces to one single term (e.g. $\sum_K S_{2c(J)-1,K} = S_{2c(J)-1,J}$). Thus, one can show that the terms under brackets in the above equation all vanish so that the scheme reduces to:

$$ \bar{Q}_I^{k+1} = \left(1 - \frac{s_1 \Delta t}{\Delta X}\right) \bar{Q}_I^k + \frac{s_1 \Delta t}{\Delta X} \bar{Q}_{I-1}^k $$

It then appears that, as for the DGMPM combined with the Forward Euler time integration [1], the use of one particle per cell and RK2 time discretization within the DGMPM leads, for one-dimensional problems, to the well-known first-order upwind method.
## 3.2 Validation of the discrete operator

Given the mathematical complexity of the scheme equation (24), the matrix operators $H$ derived above are validated with numerical results before studying the stability properties of the DGMPM. Let us then consider a one-dimensional medium $x \in [0, l]$ with $l = 1 m$, submitted to the Dirichlet condition $q(0, t) = q^d$, to an outflow boundary condition at $x = l$, and to the initial condition $q(x, 0) \equiv 0$. The combination of the scalar linear advection equation in that domain and the aforementioned boundary and initial conditions leads to a Picard problem whose exact solution is [25]:

$$q_{\text{exact}}^{1D}(x, t) = \begin{cases} q^d & \text{if } x \le s_1 t \\ 0 & \text{otherwise} \end{cases}$$

The linear advection equation is discretized by means of the DGMPM coupled with the Euler or RK2 time integrators. In parallel, the DGMPM discrete operators $H$ are explicitly constructed so that the problem is solved at each time step as in equation (24). For each numerical result, the $L^1$ norm of the relative error with the exact solution, which is suitable for discontinuous solutions [3], is computed as:

$$\epsilon_q(t^n) = \Delta x \frac{\sum_{I=1}^{N_p} |q_{\text{exact}}^{1D}(x_I, t^n) - Q_I^n|}{\sum_{J=1}^{N_p} |q_{\text{exact}}^{1D}(x_J, t^n)|}$$

Both errors are expected to be identical from one to another and to decrease with the size of the background mesh for a fixed distribution of material points in the cells. In what follows two or four particles placed symmetrically with respect to the cell centers are considered in every element. The material points are moreover regularly spaced in the whole domain. In addition, the CFL number is set to 0.43 for two particles per cell and to 0.23 for four particles per cell when the Euler algorithm is used, and to 0.50 in both cases for the RK2 time discretization. As we shall see in section 3.4, those values ensure the stability of the scheme.

It must be emphasized that the results presented below aim at showing the consistency of the discrete operators derived in section 3.1 since the study of the accuracy of DGMPM schemes is out of the scope of the present paper.

<figure><img src="image_6.png" /><figcaption>Figure 5: Convergence curves of the DGMPM using either the Euler or the RK2 time integrator for two distributions of material points in a one-dimensional mesh. Dashed and dotted curves: solution of the DGMPM scheme; Markers: solutions provided by the discrete operators.</figcaption></figure>

Figure 5 shows the evolution of the $L^1$ norm of the relative error computed at time $t = l/s_1$ ($s_1 = 10 \text{ m.s}^{-1}$) with respect to
the cell size for the two aforementioned distributions of particles. First, it can be seen that both DGMPM results converge to the exact solution with a rate of approximately one, which is an already observed result [2]. Second, the error resulting from the use of the discrete operators perfectly fits that of the DGMPM computations, which validates the expressions derived in the previous section.

## 3.3 The von Neumann linear stability analysis

The computational domain is now repeated periodically by mapping it to the domain $[-l, 0]$ so that the solution at material point $I$ and time step $k$ is expanded into a discrete Fourier basis of $2E + 1$ harmonics over the domain $x \in [-l, l]$:

$$ \bar{Q}_I^k = \sum_{j=-E}^{E} A_j^k e^{iI\lambda_j\Delta x} $$

with $A_j^k$, the magnitude of the $j$th harmonic at time step $k$, $i = \sqrt{-1}$, and $\lambda_j$ the wave number. Introduction of this expansion in equation (24) yields:

$$ A_j^{k+1} e^{iI\lambda_j \Delta x} = \sum_{L=1}^{N_p} A_j^k H_{IL} e^{iL\lambda_j \Delta x} \quad \forall j = -E, \dots, E $$

The amplification factor between two time steps at a given point is defined as:

$$ \frac{A_j^{k+1}}{A_j^k} = \sum_{L=1}^{N_p} e^{i(L-I)\lambda_j \Delta x} H_{IL} \quad \forall j = -E, \dots, E \qquad (25) $$

A necessary condition to ensure the stability of a numerical scheme is that the amplification factor must be lower than or equal to one in modulus: $|A^{k+1}/A^k| \le 1$. This upper bound allows to prevent an increasing error during the computation. For expression (25), this leads to:

$$ \left| \sum_{L=1}^{N_p} e^{i(L-I)\lambda_j \Delta x} H_{IL} \right| \le \sum_{L=1}^{N_p} \left| e^{i(L-I)\lambda_j \Delta x} H_{IL} \right| = \sum_{L=1}^{N_p} |H_{IL}| \quad \forall j = -E, \dots, E $$

where the triangle inequality, and the unit modulus of the complex number $e^{i(L-I)\lambda_j \Delta x}$ have been used. Hence, the DGMPM scheme is stable for a given discretization if the Courant number is set so that the following condition is satisfied for all material points:

$$ \sum_{L=1}^{N_p} |H_{IL}| \le 1 \quad \forall I = 1, \dots, N_p \qquad (26) $$

The stability of the scheme is thus ensured by using the lowest CFL number satisfying (26). Note however that the use of the triangle inequality leads to a more severe constraint than that really admissible. As a consequence, the Courant number can be set in practice to higher values than that resulting from the solution of (26).

## 3.4 Evaluation of the critical Courant number for particular space discretizations

The stability condition (26) can be very hard to solve analytically for general discretizations. Nevertheless, the Courant number can be set to unity for one particle lying in every cells, regardless of their positions, since in that case the DGMPM scheme identifies to the first-order upwind method. On the other hand, though the infinity of possible distributions of material points prevents the explicit derivation of a general stability condition, the optimal CFL satisfying the equality in (26) can be found numerically. To do so, the discrete DGMPM operator $H$ in equation (24) is constructed numerically. Then, the highest CFL number ensuring that all the sums of the absolute value of the row entries are lower than 1 is sought. Some configurations are studied in table 1 where the critical Courant number resulting from the two time discretizations studied above are compared. Those results have been obtained by using the same particles distribution in every elements of a one-dimensional regular mesh. The following configurations are considered:

(i) particles are positioned symmetrically with respect to cell centers and regularly spaced in the mesh. This space discretization is referred to as the natural configuration in the following.

(ii) particles in the natural configuration are all shifted of $u = \Delta x/10$.

(iii) the same than (ii) with $u$ so that one material point overlaps every left nodes of cells.

(iv) the same than (iii) for right nodes.

(v) particles are placed symmetrically with respect to cell centers but not regularly spaced in the mesh. Material points in the left half of cells are shifted of $u_1 = -\Delta x/10$ while those in the right half are shifted of $u_2 = \Delta x/10$.
<table><caption>Table 1: DGMPM critical Courant number values for Euler or RK2 time integration for several one-dimensional discretizations. Black circles denote material points while white ones represent grid nodes.</caption><thead><tr><th>Number of particles</th><th>Position of particles in cell c</th><th>Euler CFL</th><th>RK2 CFL</th></tr></thead><tbody><tr><td>1</td><td>(-)</td><td>x<sub>1</sub> &isin; [x<sub>2c-1</sub>, x<sub>2c</sub>]</td><td>1.00</td><td>1.00</td></tr><tr><td>2</td><td>(i)</td><td>&#9633;●&#9633;+●&#9633;○</td><td>0.43</td><td>1.00</td></tr><tr><td>2</td><td>(ii)</td><td>&#9633;●&#9633;+●&#9633;○</td><td>0.40</td><td>0.55</td></tr><tr><td>2</td><td>(iii)</td><td>●&#9633;●&#9633;○</td><td>0.50</td><td>0.61</td></tr><tr><td>2</td><td>(iv)</td><td>&#9633;●&#9633;+●&#9633;●</td><td>0.30</td><td>0.36</td></tr><tr><td>2</td><td>(v)</td><td>&#9633;●&#9633;+●&#9633;○</td><td>0.27</td><td>1.00</td></tr><tr><td>2</td><td>(vi)</td><td>●&#9633;+●&#9633;●</td><td>0.00</td><td>1.00</td></tr><tr><td>3</td><td>(i)</td><td>&#9633;●&#9633;●&#9633;○</td><td>0.29</td><td>1.00</td></tr><tr><td>3</td><td>(ii)</td><td>&#9633;●&#9633;+●&#9633;○</td><td>0.26</td><td>0.32</td></tr><tr><td>3</td><td>(iii)</td><td>●&#9633;●&#9633;+●&#9633;○</td><td>0.33</td><td>0.73</td></tr><tr><td>3</td><td>(iv)</td><td>&#9633;●&#9633;+●&#9633;●</td><td>0.22</td><td>0.26</td></tr><tr><td>3</td><td>(v)</td><td>&#9633;●&#9633;●&#9633;○</td><td>0.13</td><td>1.00</td></tr><tr><td>3</td><td>(vi)</td><td>●&#9633;●&#9633;●&#9633;●</td><td>0.00</td><td>1.00</td></tr><tr><td>4</td><td>(i)</td><td>&#9633;●&#9633;+●&#9633;●○</td><td>0.23</td><td>1.00</td></tr><tr><td>4</td><td>(ii)</td><td>&#9633;●&#9633;●&#9633;●○</td><td>0.19</td><td>0.22</td></tr><tr><td>4</td><td>(iii)</td><td>●&#9633;●&#9633;○</td><td>0.25</td><td>0.79</td></tr><tr><td>4</td><td>(iv)</td><td>&#9633;●&#9633;●&#9633;●</td><td>0.18</td><td>0.21</td></tr><tr><td>4</td><td>(v)</td><td>●&#9633;●&#9633;+●&#9633;○</td><td>0.05</td><td>1.00</td></tr><tr><td>4</td><td>(vi)</td><td>●&#9633;●&#9633;+●&#9633;●</td><td>0.00</td><td>1.00</td></tr></tbody></table>

(vi) the same than (v) with the first and last particles overlapping left and right nodes of cells respectively.

First, table 1 shows that the natural configuration leads to the optimal stability bound for the RK2 integrator while the CFL number allowed by using Euler time discretization decreases with increasing number of particles per cell. Second, moving every points rightward in the mesh (i.e. cases (ii) and (iv)) causes a drop in the critical Courant number for both RK2 and Euler algorithms. Nonetheless, the stability condition provided by the RK2 always remains less restrictive than that of the Euler. Third, the leftward shifting (i.e. case (iii)) leads on the other hand to an improvement of the stability condition for Euler time integration compared to that enabled in the natural configuration. This property illustrates the upwind nature of the scheme since $s_1 > 0$. The CFL number provided by RK2 integration also decreases due to the shifting while remaining higher than Euler one. At last, particle distributions conserving the symmetry with respect to cells centers (i.e. cases (v) and (vi)), yield the optimal stability condition for the RK2 integrator while the Euler CFL depends on the spacing between material points. More specifically, Euler algorithm leads to a vanishing CFL for the case (vi), thus preventing any simulation.

Generally speaking, one sees that for the examples treated here, the critical Courant number of the DGMPM-RK2 is higher than that of the DGMPM-Euler. This results is however expected since the stability region of the RK2 algorithm contains that of the forward Euler time integrator [26], and since these two approaches apply here to the same space discretization.

The same approach is now repeated with random distributions of particles. First, the number of particles lying in the cells is randomly selected in the range [1, 2, 3, 4]. Then, the positions of the material points in the cells, which is assumed to be repeated
in the whole mesh, are also randomly generated following a uniform distribution, that is: $x^p \sim \text{unif}(x_{2c(I)-1}, x_{2c(I)}) \forall p$. The space discretizations thus constructed are referred to as the *periodic* configurations hereinafter. Second, two different distributions in terms of number and position of particles are randomly generated in a cell and in its upwind neighbor. Those configurations are denoted by *non periodic* in what follows. In both cases, 1000 samples are considered for which the critical Courant number of the DGMPM-Euler and the DGMPM-RK2 are computed. The present study aims at giving some mean values of the stability bound as well as at highlighting the effects of the non-periodicity of the particle distributions between adjacent cells, which has been omitted so far.

Figure 6 shows the distributions of critical CFL for the situations mentioned above. It can first be seen that the DGMPM-

<figure><img src="image_7.png" /><figcaption>Figure 6: Critical Courant number for randomly generated distributions of material points in a regular one-dimensional mesh. The vertical lines stand for the median of the distributions.</figcaption></figure>

Euler leads to similar results for the periodic and non-periodic configurations whereas the DGMPM-RK2 is more sensitive to the difference between adjacent cells, although the impact is rather slight. This has an influence on the mean value of the critical CFL for the whole samples, which is 0.59 in any case with the forward Euler integrator, and 0.77 for periodic configurations and 0.76 for non-periodic ones with the RK2. On the other hand, the first-order time integration leads to more occurrences in the low values of the critical Courant number (i.e. CFL≤ 0.5), which is shown by the median of the distributions depicted with vertical lines. As a result, the above results confirm that in general, the RK2 time integrator provides a better stability bound to the DGMPM. At last, even though some prohibitive configurations for which the CFL is severely restricted are highlighted in table 1, figure 6 shows that such cases are rather rare.

Notice that although the stability depends on the position of particles, the optimal Courant number can be reached even with the Euler discretization, whereas the classical DGFEM scheme presented in [27] is restricted to condition $\Delta t/\Delta x = O(\sqrt{\Delta x})$. This limitation has been addressed by introducing slope limiters in order to remove non-physical oscillations in the vicinity of sharp solutions while providing high-order accuracy in smooth regions [28]. However, the stability of the method is still bounded by $CFL \le 1/2$ and the scheme is first-order accurate. Alternatively, the use of a second-order Runge-Kutta [29] allows to increase the stability bound to $CFL \le 1/3$ for first-order polynomials. The combination of the RK2 time discretization and slope limiters moreover yields second-order accuracy for linear polynomials while ensuring the stability for $CFL \le 1/3$ [11]. It is noteworthy that space-time DGFEM formulations [12, 30] provided more recently the ability to relax constraints of pure space DGFEM and obtained a critical Courant number set at 1.

# 4 Stability properties of the two-dimensional schemes

## 4.1 The scheme equation

We now move on to the scalar linear advection equation with positive coefficients $s_1, s_2 > 0$. Moreover, the space coordinates $x_1, x_2$ in the $(\vec{e}_1, \vec{e}_2)$ plane are written $x$ and $y$. Then, the physical domain under consideration, that is $x, y \in [0, l] \times [0, h]$, is discretized with $N_p$ material points arbitrarily distributed in a Cartesian grid made of $E$ four-node bilinear elements with
constant size $\Delta x \times \Delta y$. With positions of nodes of cell C denoted by $\vec{x}_i^C = [x_i^C \ y_i^C]$, as depicted in figure 7, the location $\vec{x}'$ of an

<figure><img src="image_8.png" /><figcaption>Figure 7: Parent and current configuration of a rectangular four-node bilinear element</figcaption></figure>

arbitrary point in cell C maps to the parent coordinates $(\xi, \eta)$ in the domain $[-1, 1] \times [-1, 1]$ according to:

$$ \xi = 2 \frac{x - x_1^C}{\Delta x} - 1 \quad ; \quad d\xi = 2 \frac{dx}{\Delta x} \qquad (27) $$

$$ \eta = 2 \frac{y - y_1^C}{\Delta y} - 1 \quad ; \quad d\eta = 2 \frac{dy}{\Delta y} $$

Horizontal and vertical edge lengths are distinguished here in spite of the Cartesian nature of the grid in order to easily extend the following study to rectilinear grids. Again, there is no empty cell inside the physical domain so that no hole is generated, and periodic boundary conditions are considered to simplify the analysis. Though one can imagine to combine the DGMPM discretization with a multi-stage time integration as proposed for one-dimensional problems, the equations one has to deal with become complicated. Hence, the analysis of the DGMPM scheme for two-dimensional problems carried out here only considers the Euler time discretization.

The updated solution at material point I is obtained by interpolation of updated nodal quantities:

$$ \bar{Q}_I^{k+1} = \sum_{i=1}^{4E} S_{iI} \bar{q}_i^{k+1} \qquad (28) $$

where the $\bar{q}_i^{k+1}$ result from the discrete system:

$$ \bar{q}_i^{k+1} = \bar{q}_i^k + \frac{\Delta t}{M_i^L} (K_{ij}^x s_1 \bar{q}_j^k + K_{ij}^y s_2 \bar{q}_j^k - \hat{F}_i^k) $$

The lumped mass matrix in the above equation has the same expression as in the one-dimensional case that depends on the shape functions of the four-node bilinear element: $M_i^L = \sum_K m_K S_{iK}$. Making use of parent coordinates (27), the pseudo-stiffness matrices read:

$$ K_{ij}^x = \sum_j \frac{\partial S_{iJ}}{\partial x} m_J S_{jJ} = \frac{2}{\Delta x} \sum_j \frac{\partial S_{iJ}}{\partial \xi} m_J S_{jJ} $$

$$ K_{ij}^y = \sum_j \frac{\partial S_{iJ}}{\partial y} m_J S_{jJ} = \frac{2}{\Delta y} \sum_j \frac{\partial S_{iJ}}{\partial \eta} m_J S_{jJ} $$

As for one-dimensional cases, the homogeneous medium yields the same mass for every particles so that, by writing $\frac{\partial(\bullet)}{\partial\xi} = \partial_{\xi}(\bullet)$, one gets:

$$ \frac{K_{ij}^x}{M_i^L} = \frac{2}{\Delta x} \frac{\sum_J \partial_{\xi} S_{iJ} S_{jJ}}{\sum_K S_{iK}} $$

$$ \frac{K_{ij}^y}{M_i^L} = \frac{2}{\Delta y} \frac{\sum_J \partial_{\eta} S_{iJ} S_{jJ}}{\sum_K S_{iK}} $$
The nodal solutions in cell C at time k being given by the projection $\bar{q}_i^{C,k} = \frac{\sum_L S_{iL} \bar{Q}_L^k}{\sum_K S_{iK}}$, volume fluxes of the discrete form can be rewritten as:

$$ 
\begin{aligned}
s_1 \frac{K_{ij}^x}{M_i^L} \bar{q}_j^k &= \sum_L \bar{Q}_L^k \frac{2}{\Delta x} \frac{s_1 \sum_J \partial_\xi S_{iJ} \sum_j S_{jJ} S_{jL}}{\sum_K S_{iK} \sum_K S_{jK}} \\
s_2 \frac{K_{ij}^y}{M_i^L} \bar{q}_j^k &= \sum_L \bar{Q}_L^k \frac{2}{\Delta y} \frac{s_2 \sum_J \partial_\eta S_{iJ} \sum_j S_{jJ} S_{jL}}{\sum_K S_{iK} \sum_K S_{jK}}
\end{aligned} 
\quad (29) $$

Then, the nodal interface flux $\hat{F}_i^k$ results from the integration of Godunov's fluxes along edges the node belongs to, according to the weak form (5). Referring to a quantity defined at an interface by means of parenthesis superscripts, the Godunov's flux at interface (i) is:

$$ F_n^{(i)} = s_n q^{-,(i)} = \underbrace{s_n q^{+,(i)}}_{F_n(q^{+,(i)})} - \underbrace{s_n (q^{+,(i)} - q^{-,(i)})}_{A_{U/D}^+} \quad (30) $$

where $s_n$ is the speed in the normal direction to the interface (i.e. $s_2$ for horizontal edges and $s_1$ for vertical edges). Equation (30) further involves state vectors $q^{-,(i)}$ and $q^{+,(i)}$ obtained by averaging nodal values connected to interface (i) on Upwind and Downwind sides respectively, and the right-going fluctuation $A_{U/D}^+$. The CTU is adopted by adding transverse corrections to fluxes (30), based on those fluctuations according to equation (12):

$$ \mathcal{B}^+ A_{U/D}^+ = s_t s_n (q^{+, (i)} - q^{-, (i)}) $$

with $s_t$ the speed in the tangent direction to the interface. The final expression of intercell fluxes is hence:

$$ F_n^{(i)} = s_n q^{-,(i)} - s_t s_n \frac{\Delta t}{2\Delta x} (q^{+, (i)} - q^{-, (i)}) \quad (31) $$

Figure 8 shows transverse corrections in the cell C based on fluctuations coming from Bottom (B), Left (L) and Bottom Left

<figure><img src="image_9.png" /><figcaption>Figure 8: Two-dimensional patch of cells of constant size $\Delta x \times \Delta x$.</figcaption></figure>

(BL) neighbor elements. The use of the numbering of interfaces and nodes adopted in figure 8 allows the specialization of equation (31) to the intercell fluxes of cell C:

$$ F_n^{(1)} = s_2 \frac{q_3^{B,k} + q_4^{B,k}}{2} - s_1 s_2 \frac{\Delta t}{2\Delta y} \left( \frac{q_1^{B,k} + q_4^{B,k}}{2} - \frac{q_2^{BL,k} + q_3^{BL,k}}{2} \right) $$

$$ F_n^{(2)} = s_1 \frac{q_2^{C,k} + q_3^{C,k}}{2} - s_1 s_2 \frac{\Delta t}{2\Delta x} \left( \frac{q_1^{C,k} + q_2^{C,k}}{2} - \frac{q_3^{B,k} + q_4^{B,k}}{2} \right) $$

$$ F_n^{(3)} = s_2 \frac{q_3^{C,k} + q_4^{C,k}}{2} - s_1 s_2 \frac{\Delta t}{2\Delta y} \left( \frac{q_1^{C,k} + q_4^{C,k}}{2} - \frac{q_2^{L,k} + q_3^{L,k}}{2} \right) $$

$$ F_n^{(4)} = s_1 \frac{q_2^{L,k} + q_3^{L,k}}{2} - s_1 s_2 \frac{\Delta t}{2\Delta x} \left( \frac{q_1^{L,k} + q_2^{L,k}}{2} - \frac{q_3^{BL,k} + q_4^{BL,k}}{2} \right) $$
Denoting the number of particles in cell C and the mass they carry by $N_p^C$ and $m^C$ respectively, the mass density reads $\rho = \frac{N_p^C m^C}{\Delta x \Delta y}$. Thus, introduction of the particle fields projection yields the following expressions for interface fluxes:

$$
\begin{align*}
F_n^{(1)} &= \sum_{J=1}^{N_p} \bar{Q}_J^k \frac{s_2 N_p^C m^C}{2 \Delta x \Delta y} \left[ \left( \frac{S_{3J}^B}{\sum_L S_{3L}^B} + \frac{S_{4J}^B}{\sum_L S_{4L}^B} \right) - s_1 \frac{\Delta t}{2 \Delta y} \left( \frac{S_{1J}^B}{\sum_L S_{1L}^B} + \frac{S_{4J}^B}{\sum_L S_{4L}^B} - \frac{S_{2J}^{BL}}{\sum_L S_{2L}^{BL}} - \frac{S_{3J}^{BL}}{\sum_L S_{3L}^{BL}} \right) \right] \\
F_n^{(2)} &= \sum_{J=1}^{N_p} \bar{Q}_J^k \frac{s_1 N_p^C m^C}{2 \Delta x \Delta y} \left[ \left( \frac{S_{2J}^C}{\sum_L S_{2L}^C} + \frac{S_{3J}^C}{\sum_L S_{3L}^C} \right) - s_2 \frac{\Delta t}{2 \Delta x} \left( \frac{S_{1J}^C}{\sum_L S_{1L}^C} + \frac{S_{2J}^C}{\sum_L S_{2L}^C} - \frac{S_{3J}^B}{\sum_L S_{3L}^B} - \frac{S_{4J}^B}{\sum_L S_{4L}^B} \right) \right] \\
F_n^{(3)} &= \sum_{J=1}^{N_p} \bar{Q}_J^k \frac{s_2 N_p^C m^C}{2 \Delta x \Delta y} \left[ \left( \frac{S_{3J}^C}{\sum_L S_{3L}^C} + \frac{S_{4J}^C}{\sum_L S_{4L}^C} \right) - s_1 \frac{\Delta t}{2 \Delta y} \left( \frac{S_{1J}^C}{\sum_L S_{1L}^C} + \frac{S_{4J}^C}{\sum_L S_{4L}^C} - \frac{S_{2J}^L}{\sum_L S_{2L}^L} - \frac{S_{3J}^L}{\sum_L S_{3L}^L} \right) \right] \\
F_n^{(4)} &= \sum_{J=1}^{N_p} \bar{Q}_J^k \frac{s_1 N_p^C m^C}{2 \Delta x \Delta y} \left[ \left( \frac{S_{2J}^L}{\sum_L S_{2L}^L} + \frac{S_{3J}^L}{\sum_L S_{3L}^L} \right) - s_2 \frac{\Delta t}{2 \Delta x} \left( \frac{S_{1J}^L}{\sum_L S_{1L}^L} + \frac{S_{2J}^L}{\sum_L S_{2L}^L} - \frac{S_{3J}^{BL}}{\sum_L S_{3L}^{BL}} - \frac{S_{4J}^{BL}}{\sum_L S_{4L}^{BL}} \right) \right]
\end{align*}
$$

written for simplicity:

$$F_n^{(i)} = \sum_{J}^{N_p} \bar{Q}_J^k \frac{s_n N_p^C m^C}{2 \Delta x \Delta y} [\phi_J^{(i)} + \phi_J^{(i),T}] \quad (32)$$

In the latter expressions, $\phi^{(i)}$ is devoted to normal contributions while $\phi^{(i),T}$ stands for transverse corrections at interface (i). Numerical fluxes considered above are based on normal vectors oriented in the direction of the stream (see figure 8). Nodal interface fluxes on the other hand, as defined in the semi-discrete system:

$$\hat{F}_i = \int_l S_i(\vec{x}) F_n dl$$

are based on the outgoing flux to an element so that $F_n^{(1)}$ and $F_n^{(4)}$ must be counted negatively. The integral for cell C is then:

$$\hat{F}_i = - \int_{x_1^C}^{x_2^C} S_i(x, y_1^C) F_n^{(1)}(dx) + \int_{y_2^C}^{y_3^C} S_i(x_2^C, y) F_n^{(2)}(dy) + \int_{x_2^C}^{x_3^C} S_i(x, y_3^C) F_n^{(3)}(dx) - \int_{y_1^C}^{y_4^C} S_i(x_1^C, y) F_n^{(4)}(dy)$$

which can be computed analytically by using parent coordinates (27):

$$\hat{F}_1 = - \frac{1}{2} [\Delta x F_n^{(1)} + \Delta y F_n^{(4)}] \quad ; \quad \hat{F}_2 = - \frac{1}{2} [\Delta x F_n^{(1)} - \Delta y F_n^{(2)}]$$

$$\hat{F}_3 = \frac{1}{2} [\Delta x F_n^{(3)} + \Delta y F_n^{(2)}] \quad ; \quad \hat{F}_4 = \frac{1}{2} [\Delta x F_n^{(3)} - \Delta y F_n^{(4)}]$$

A condensed way of writing those fluxes is adopted by means of the middle point of edge (j) with coordinates $\vec{x}^{(j)}_{1/2}$, at which the shape functions are:

$$S_i(\vec{x}^{(j)}_{1/2}) = \begin{cases} \frac{1}{2} & \text{if node i belongs to edge (j)} \\ 0 & \text{otherwise.} \end{cases}$$

In addition, components of the outward normal vector to edges $n_x^{(i)}$ and $n_y^{(i)}$ allow to take into account different signs of intercell fluxes in the Cartesian grid. One thus writes:

$$\hat{F}_i = \frac{1}{2} \sum_{j=1}^{\text{Nedges}} 2S_i(\vec{x}^{(j)}_{1/2}) (\Delta yn_x^{(j)} + \Delta xn_y^{(j)}) F_n^{(j)}$$

which, combined with equation (32) leads to:

$$\hat{F}_i^k = \sum_J^{N_p} \bar{Q}_J^k \sum_{j=1}^{\text{Nedges}} S_i(\vec{x}^{(j)}_{1/2}) (s_1 \Delta yn_x^{(j)} + s_2 \Delta xn_y^{(j)}) \frac{N_p^C m^C}{2 \Delta x \Delta y} [\phi_J^{(i)} + \phi_J^{(i),T}]$$

These terms are divided by the lumped mass matrix in the discrete form:

$$\frac{\hat{F}_i^k}{M_i^L} = \sum_J \frac{\bar{Q}_J^k}{\sum_K S_{iK}} \sum_{j=1}^{\text{Nedges}} \frac{1}{2} S_i(\vec{x}^{(j)}_{1/2}) N_p^C \left( s_1 \frac{n_x^{(j)}}{\Delta x} + s_2 \frac{n_y^{(j)}}{\Delta y} \right) [\phi_J^{(j)} + \phi_J^{(j),T}] \quad (33)$$
At last, gathering the mapping of updated nodal quantities to the particles (28), expressions of volume fluxes (29) and intercell ones (33), the updated value at material point I contained in cell C reads:

$$ \bar{Q}_I^{k+1} = \sum_{L=1}^{N_p} \bar{Q}_L^k \sum_{i=1}^{4E} \frac{S_{iI}}{\sum_K S_{iK}} \left\{ S_{iL} + 2 \sum_{j=1}^{4E} \frac{S_{jL}}{\sum_K S_{jK}} \sum_{J=1}^{N_p} S_{JJ} \left[ s_1 \frac{\Delta t}{\Delta x} \partial_{\xi} S_{iJ} + s_2 \frac{\Delta t}{\Delta y} \partial_{\eta} S_{iJ} \right] \right. \\ \left. - \frac{1}{2} \sum_{p=1}^{\text{edges}} S_i (\vec{x}_{1/2}^{(p)}) N_p^C \left( s_1 \frac{\Delta t}{\Delta x} n_x^{(j)} + s_2 \frac{\Delta t}{\Delta y} n_y^{(j)} \right) \left[ \phi_L^{(p)} + \phi_L^{(p),T} \right] \right\} \quad (34) $$

Recall that transverse contributions $\phi_L^{(p),T}$ depend on $\Delta t$, thus providing second-order corrections in the two-dimensional scheme equation (34), which can also be rewritten as:

$$ \bar{Q}_I^{k+1} = \sum_{L=1}^{N_p} \bar{Q}_L^k H_{IL} \quad (35) $$

## 4.2 Validation of the discrete operator

As for one-dimensional problems, the DGMPM discrete operators $H$ derived in section 4.1 are validated on the scalar linear advection equation. Therefore, a two-dimensional medium $(x, y) \in [0, l] \times [0, l]$, with $l = 4m$, is considered. This domain is submitted to the Dirichlet conditions $q(x = 0, y, t) = q^d$ and $q(x, y = 0, t) = 0$, the other ends being outflow boundaries, and to the initial condition $q(x, y, 0) \equiv 0$, in such a way that the exact solution reads [3]:

$$ q_{\text{exact}}^{2D}(x, y, t) = \begin{cases} q^d & \text{if } x \le s_1 t \text{ and } y \ge s_2 t \\ 0 & \text{otherwise} \end{cases} $$

The linear advection equation is solved by means of the DGMPM scheme based on the DCU or the CTU on the one hand, and by using the discrete operators as in equation (35) on the other hand. Then, for all the numerical solutions the $L^1$ norm of the relative error is computed as:

$$ \epsilon_q(t^n) = \Delta x^2 \frac{\sum_{I=1}^{N_p} |q_{\text{exact}}^{2D}(x_I, y_I, t^n) - Q_I^n|}{\sum_{J=1}^{N_p} |q_{\text{exact}}^{2D}(x_J, y_J, t^n)|} $$

The medium is discretized with two distributions of particles satisfying vertical and horizontal symmetries in every cell, as depicted in figure 9, are considered. Both of those space discretizations lead to overlapping regular computational and material points grids. Thus, the evolution of the error between numerical solutions and the exact one is studied at $t = l/s_1$ when the background grid is refined, while keeping the same distribution of particles in the elements.

<figure><img src="image_10.png" /><figcaption>Figure 9: Material point distributions used for validating the two-dimensional DGMPM discrete operators $H$. The circles denote the nodes of the computational grid while the crosses denote material points.</figcaption></figure>

In figure 10, convergence curves are depicted in the arbitrary case of $s_1 = s_2 = 10 \text{ m.s}^{-1}$. As can be seen, the DGMPM-DCU and the DGMPM-CTU yield solutions that converge to the exact solution with a rate of approximately one. In addition, the results provided by the discrete operators $H$, which fit the solutions of the complete DGMPM procedures, show that the developments of section 4.1 are valid.

## 4.3 The von Neumann linear stability analysis

Analogously to the one-dimensional case, the solution at a material point can be expanded into a discrete Fourier basis over the domain $[-l, l] \times [-h, h]$. We consider here a structured distribution of particles made of $N_p = N_p^x \times N_p^y$ material points so that one
<figure><img src="image_11.png" /><figcaption>Figure 10: Convergence curves of the DGMPM using either the DCU or the CTU method for two distributions of material points in a two-dimensional mesh. Dashed and dotted curves: solution of the DGMPM scheme; Markers: solutions provided by the discrete operators.</figcaption></figure>

can denote the solution at particles by $\bar{Q}_{IL}$, where I and L are the row and column indices of material points. For one arbitrary Fourier mode, one has [3]:

$$ \bar{Q}_{IL}^k = A_{jq}^k e^{i(I\lambda_j\Delta y + L\lambda_q\Delta x)} $$

where $i = \sqrt{-1}$, and $\lambda_j$ and $\lambda_q$ are wave numbers in directions x and y respectively. Then, the amplification factor reads:

$$ \frac{A_{jq}^{k+1}}{A_{jq}^k} = \sum_{K=1}^{N_p^y} \sum_{J=1}^{N_p^x} e^{i([I-K]\lambda_j\Delta y + [L-J]\lambda_q\Delta x)} H_{IL,KJ} $$

The requirement that the absolute value of the amplification factor is lower than or equal to one leads to the following stability condition:

$$ \left| \frac{A_{jq}^{k+1}}{A_{jq}^k} \right| = \left| \sum_{K=1}^{N_p^y} \sum_{J=1}^{N_p^x} e^{i([I-K]\lambda_j \Delta y + [L-J]\lambda_q \Delta x)} H_{IL,KJ} \right| \le 1 \Leftrightarrow \sum_{K=1}^{N_p^y} \sum_{J=1}^{N_p^x} |H_{IL,KJ}| \le 1 $$
$$ \sum_{L=1}^{N_p} |H_{IL}| \le 1 \quad \forall I = 1, \dots, N_p \qquad (36) $$

## 4.4 Evaluation of the critical Courant number for particular space discretizations

Again the study of stability conditions for general discretizations is complicated given the scheme equation (34). Nevertheless, it can first be noticed that the single particle-per-cell discretization leads to a piece-wise constant reconstruction of the field on the computational grid after the projection from material points to nodes whatever their locations in cells. Therefore, the first-order
upwind method is once again retrieved. This method is known to be characterized by [3]:

$$|s_1| \frac{\Delta t}{\Delta x} + |s_2| \frac{\Delta t}{\Delta y} \le 1 \quad \text{for DCU} \qquad (37a)$$

$$\max \left( |s_1| \frac{\Delta t}{\Delta x}, |s_2| \frac{\Delta t}{\Delta y} \right) \le 1 \quad \text{for CTU} \qquad (37b)$$

From now on, it is assumed that $s_1 \ge s_2 > 0$ and regular cells $\Delta x = \Delta y$ are considered. In that case, the critical time steps resulting from conditions (37a) and (37b) read:

$$\Delta t \le \frac{\Delta x}{s_1 + s_2} \quad \text{for DCU} ; \quad \Delta t \le \frac{\Delta x}{s_1} \quad \text{for CTU}$$

In terms of Courant number, the previous conditions are:

$$s_1 \frac{\Delta t}{\Delta x} \le \frac{s_1}{s_1 + s_2} \quad \text{for DCU} \qquad (38a)$$

$$s_1 \frac{\Delta t}{\Delta x} \le 1 \quad \text{for CTU} \qquad (38b)$$

Hence, the CFL number can be set at one by using the CTU regardless of the speeds values. Conversely, the maximal Courant number governing the DCU tends to one if and only if $s_1 \gg s_2$.

Configurations involving more particles in the computational grid cells are now studied numerically by assuming the same distribution of material points in every elements. Since the coefficients $H_{IJ}$ depend on both horizontal and vertical wave speeds, the scheme equation (35) can be written as a function of the CFL number by means of the speed ratio $s_1/s_2$. Hence, the maximal Courant number satisfying the stability condition (36) also depends on the speed ratio. Evolutions of the CFL numbers corresponding to several distributions of particles in a two-dimensional grid are gathered in tables 2, 3 and 4 for the DGMPM scheme using DCU and CTU methods. The first column of these tables shows the positions of material points inside cells for discretizations based on two, four or nine particles per element respectively.

<figure><img src="image_12.png" /><figcaption>Table 2: Values of critical Courant number $s_1 \frac{\Delta t}{\Delta x}$ for two-dimensional DGMPM scheme using either DCU or CTU with respect to the locations of the two material points lying in every cells, as a function of the speeds ratio $s_1/s_2$.</figcaption></figure>

The space discretization leading to two particles lying in every cells of the mesh is such that within an element, the two material points are both either on the horizontal axis or on the vertical axis of the cell. The aforementioned cases respectively correspond to the results reported in first and second rows of table 2. Two situations are then to be distinguished:
* Material points are regularly-spaced within the grid and placed symmetrically two-by-two with respect to cells centers. Those distributions are drawn in the first column of table 2 by using blue plus signs to represent particles.

* Material points still satisfy symmetry in cells, but are no longer regularly-spaced in the mesh. In that case, particles are drawn with red crosses.

First, the results of table 2 show that the CFL number exhibits a non-linear dependence on the speed ratio $s_1/s_2$, which asymptotically approaches some value depending on the particles distribution. The configuration in which the particles are regularly spaced in the grid is of particular interest. Indeed, this discretization is equivalent to the vertical repetition of the natural configuration for two particles per cell studied in section 3. As a result, the one-dimensional result for Euler discretization CFL≈ 0.43 should be recovered for $s_1 \gg s_2$. This value is depicted with black dashed lines in the plots of first row in table 2. It is thus seen that the one-dimensional value of the critical Courant number corresponds to the asymptotic limit for two-dimensional problems as $s_1/s_2 \to \infty$.

Next, the spacing between particles is reduced while keeping the symmetry with respect to cells centers. It is then seen that, as for the one-dimensional case, the critical Courant number increases for both DCU and CTU approaches. Third, whether particles lie on the horizontal axis or the vertical axis of cells has a great influence on the critical Courant number one can expect. Hence, the configurations of the second row of table 2 yield higher CFL numbers for given speed ratios. It then appears that in order to improve the stability of the scheme, one must use a lower number of material points in the direction of the dominating wave speed than in the perpendicular one. For a Cartesian distribution of particles $N_p^x \times N_p^y$ this corresponds to $N_p^y > N_p^x$ if $s_1 > s_2$, and $N_p^x > N_p^y$ if $s_2 > s_1$. At last, it is worth noticing that the improvement brought by the CTU is much less significant than in the case of one single particle-per-cell discretization for which the Courant number can be set to one according to equations (38a) and (38b).

<figure><img src="image_13.png" /><figcaption>Table 3: Values of critical Courant number $s_1 \frac{\Delta t}{\Delta x}$ for two-dimensional DGMPM scheme using either DCU or CTU with respect to the material points distribution, as a function of the speeds ratio $s_1/s_2$. Four particles per cell.</figcaption></figure>
We now move on to cases for which grid cells each contains four material points, by considering a square shaped distribution of particles in every elements. This square distribution, which has already been used for the convergence study, is then deformed so that the following situations occur:

* the particles are displaced while keeping the center of the square of material points on the cell centroids. Three configurations are then depicted in the first row of table 3.

* the square of particles is translated horizontally without distortion as depicted in the second row of the table.

* the square of particles is translated vertically without distortion (third row of the table).

First, we see in the first row of table 3 that the contraction of the square shaped distribution of particles yields an improvement of the critical CFL number. These observations are similar to those made for one-dimensional problems, especially for the configuration in which the particles coincide with nodes that leads to a vanishing Courant number (see green curves with star marks). On the other hand, the increase in CFL number enabled by the use of the CTU approach is again less important that in the case of one particle per cell.

<figure><img src="image_14.png" /><figcaption>Table 4: Values of critical Courant number $s_1 \frac{\Delta t}{\Delta x}$ for two-dimensional DGMPM scheme using either DCU or CTU with respect to the material points distribution, as a function of the speeds ratio $s_1/s_2$. Nine particles per cell.</figcaption></figure>

Second, the last two rows of table 3 show that the translation of the square of particle inside elements may have a significant influence on the evolution of Courant's number with respect to the speed ratio. On the one hand, an upstream horizontal shifting leads to an increased critical CFL number in the asymptotic limit compared to the natural configuration, and a downstream one yields a more restrictive stability bound (see the horizontal dashed line in the figures). Note that a similar behavior was observed for the one-dimensional scheme. On the other hand, both upstream and downstream vertical translation of the material points lead to a more significant decrease in the critical Courant number for the cases studied here in which $s_1 \ge s_2$. At last,
the curves corresponding to DCU and CTU do not exhibit significant differences, meaning that the CTU seems to be mainly
efficient for the one particle per cell discretization.

Similar distributions of nine particles with square shapes are now looked at. Analogously to the four particles per cell discretization, the square pattern is either distorted or translated inside the cells and the resulting CFL numbers are gathered in table 4. The observations are similar to the four particle cases, namely, we see for the contraction of the square that the closer the particles are from cell centers, the higher the CFL number is. Moreover, upstream translations of particles yield once again a slight increase in CFL number for a given speed ratio. In addition, it can be seen that the critical CFL number is in general more restrictive than for four particles per cell, unless a vertical shifting is made. Nevertheless, since this phenomenon is not obvious for a horizontal shifting, it is not well understood for now and requires additional investigations.

As for one-dimensional problems, let us now look at randomly generated space discretizations. The number of particles lying in one cell and in its upwind neighbors is then randomly selected in the range [1, 2, 3, 4]. Next, the vertical and horizontal positions of the material points in each element follow a uniform distribution, namely: $(\xi^p, \eta^p) \sim \text{unif}(-1, 1) \ \forall p$. Periodic configurations are considered if the distribution of particles is identical in the whole patch of cells while non periodic ones may lead to as much set of particles as there is elements in the stencil (i.e. three for the DCU and four for the CTU). Since the stability bound of the scheme depends on the wave speeds involved, two different cases are considered: $s_1 = s_2 = 2 \text{ m.s}^{-1}$ and $s_1 = 10s_2 = 20 \text{ m.s}^{-1}$. The first situation has already been shown to be more penalizing for the DGMPM whereas the second one is expected to be less restrictive. The critical Courant number has been computed for 1000 samples in each case, the resulting distributions being depicted in figure 11.

<figure><img src="image_15.png" /><figcaption>Figure 11: Critical Courant number for randomly generated distributions of material points in a regular two-dimensional mesh. The vertical lines stand for the median of the distributions.</figcaption></figure>

It can be seen that increasing the stencil of the scheme by using the CTU makes the critical Courant number more sensitive to the periodicity of the collection of material points regardless of the speed ratio. Indeed, the results provided by the DCU are very close in both cases whereas the CTU leads to noticeable differences. Although similar observations have been made for the one-dimensional scheme (see figure 6), the periodicity has a more significant influence on the stability of the CTU than on DGMPM-RK2 one-dimensional scheme due to the number of cells in the stencil. On the other hand, figure 11a shows that the medians of the distributions are lower than 0.20, meaning that low values of the critical Courant number occur more frequently (see the vertical lines in the figure). Note also that the periodicity of the material points set has mainly an influence on the CTU results, for which high CFL values occur with periodic configurations but much less for non periodic ones. Indeed, the mean values of the critical Courant number for the DCU is 0.24 in any case while periodic and non periodic configurations lead
respectively to a mean values of 0.37 and 0.26 for the CTU.

Increasing the speed ratio leads to rather different distribution of the critical Courant number as can be seen in figure 11b. In that case, more high values occur (i.e. CFL≥ 0.7) whatever the periodicity of the material points set for the DCU and the CTU. Nevertheless, the median values raise to approximately 0.3 for a higher speed as can be seen in figure 11b. As a consequence, the mean values are 0.41 for the DCU, and 0.44 in case of a periodic configuration and 0.40 for non-periodic ones for the CTU. In addition, the median values provided by the DCU and the CTU undergo an increase as depicted in the figure.

The above discussion therefore highlights that one can expect a non prohibitive time-stepping with the two-dimensional DGMPM.

# 5 Concluding remarks

The formulation of the Discontinuous Galerkin Material Point Method has been recalled in this paper for the scalar linear advection equation. Numerical fluxes arising in the weak form allow to take into account the transverse propagation of waves for multi-dimensional problems, by means of the reformulation of the well-known CTU method. The semi-discrete system resulting from the combination of the material point discretization and the DG approximation has been discretized by means of both forward Euler and second-order Runge-Kutta time discretizations. Starting from those discrete systems, the scheme equations have next been derived in a finite difference sense for one and two-dimensional problems. At last, the von Neumann linear stability analysis of those scheme equations have been carried out. As a result, conditions which enable to ensure the stability of DGMPM schemes for a given discretization have been written. Based on these conditions, the critical Courant number can be numerically calculated for a given distribution of particles. The method has been applied to several one-dimensional and two-dimensional cases assuming that the distribution of particles was the same in every cells of the computational grid. Thus, it has been shown the RK2 time discretization provides higher critical CFL to the DGMPM. Moreover, for both the forward Euler and the RK2 time integrators, moving all the particles upstream (*resp. downstream*) leads to less restrictive (*resp. more restrictive*) stability conditions. In addition, for distributions of particles satisfying the symmetry with respect to the cell centers, the RK2 algorithm yields the optimal CFL number whereas it is not the case for the forward Euler. On the other hand, similar results have been emphasized for two-dimensional problems. Namely, losing the symmetries of the distribution of particles with respect to the vertical and horizontal axes of the cells has an influence on the critical Courant number. The effect of this depends on whether the translation is upstream or downstream and on the ratio between the horizontal and vertical speeds involved in the problem. Finally, the statistical analyzes performed on 1000 samples randomly generated show that the DGMPM admits satisfactory mean critical CFL values even for non periodic distributions.

The stability conditions derived in this paper allow to fully exploit the arbitrariness of the grid in order to build it so as to accurately capture waves for complex geometries or finite deformations, as suggested in the founding paper of the DGMPM [1]. The present work is also a necessary step in order to build an updated Lagrangian formulation of the method.

Although the derivation of the DGMPM has been done with first-order accuracy for the moment, its extension to higher-order approximation is possible by means of high-order polynomials. Following the development of the original MPM [20, 21, 22, 31], such an extension to high-order approximation can be carried out by modifying the reconstruction of fields on the computational grids and the quadrature rule. In that context, the lumping of the mass matrix as well as the computation of face-constant intercell fluxes should be reconsidered. This would, however, have an impact on the stability properties of the method which will have to be reevaluated for that case.

# References

[1] A. Renaud, T. Heuzé, and L. Stainier, "A discontinuous galerkin material point method for the solution of impact problems in solid dynamics," *Journal of computational physics*, 2018.

[2] A. Renaud, *The Discontinuous Galerkin Material Point Method: Application to hyperbolic problems in solid mechanics*. PhD thesis, École Centrale de Nantes, 2018.

[3] R. J. Leveque, *Finite volume methods for hyperbolic problems*. Cambridge university press, 2002.

[4] P. Colella, "Multidimensional upwind methods for hyperbolic conservation laws," *Journal of Computational Physics*, vol. 87, no. 1, pp. 171-200, 1990.

[5] V. Psyk, D. Risch, B. Kinsey, A. Tekkaya, and M. Kleiner, "Electromagnetic forming – a review," *Journal of Materials Processing Technology*, vol. 211, no. 5, pp. 787-829, 2011.

[6] T. Belytschko, W. K. Liu, B. Moran, and K. Elkhodary, *Nonlinear finite elements for continua and structures*. John wiley & sons, 2nd ed., 2014.

[7] A. Harten, "High resolution schemes for hyperbolic conservation laws," *Journal of computational physics*, vol. 49, no. 3, pp. 357-393, 1983.
[8] W. Reed and T. Hill, "Triangular mesh methods for the neutron transport equation," *Los Alamos Report LA-UR-73-479*, 1973.

[9] B. Cockburn, "Discontinuous galerkin methods for convection-dominated problems," in *High-order methods for computational physics*, pp. 69–224, Springer, 1999.

[10] B. van Leer, "Towards the ultimate conservative difference scheme. ii. monotonicity and conservation combined in a second-order scheme," *Journal of Computational Physics*, vol. 14, no. 4, pp. 361–370, 1974.

[11] B. Cockburn and C.-W. Shu, "The runge-kutta local projection $P^1$-discontinuous-galerkin finite element method for scalar conservation laws," *ESAIM: Mathematical Modelling and Numerical Analysis*, vol. 25, no. 3, pp. 337–361, 1991.

[12] G. R. Richter, "An explicit finite element method for the wave equation," *Applied Numerical Mathematics*, vol. 16, no. 1, pp. 65–80, 1994.

[13] F. H. Harlow, "The particle-in-cell method for numerical solution of problems in fluid dynamics," *Proc. Symp. Appl. Math.*, vol. 15, 03 1962.

[14] D. Sulsky, Z. Chen, and H. L. Schreyer, "A particle method for history-dependent materials," *Computer methods in applied mechanics and engineering*, vol. 118, no. 1-2, pp. 179–196, 1994.

[15] J. Brackbill and H. Ruppel, "Flip: A method for adaptively zoned, particle-in-cell calculations of fluid flows in two dimensions," *Journal of Computational Physics*, vol. 65, no. 2, pp. 314–343, 1986.

[16] A. Nishiguchi and T. Yabe, "Second-order fluid particle scheme," *Journal of Computational Physics*, vol. 52, no. 2, pp. 390–413, 1983.

[17] C. C. Hammerquist and J. A. Nairn, "A new method for material point method particle updates that reduces noise and enhances stability," *Computer Methods in Applied Mechanics and Engineering*, vol. 318, pp. 724–738, 2017.

[18] D. A. Di Pietro and A. Ern, *Mathematical aspects of discontinuous Galerkin methods*, vol. 69. Springer Science & Business Media, 2011.

[19] E. Love and D. Sulsky, "An unconditionally stable, energy-momentum consistent implementation of the material-point method," *Computer Methods in Applied Mechanics and Engineering*, vol. 195, no. 33-36, pp. 3903–3925, 2006.

[20] G. Ming, *Improving the Material Point Method*. PhD thesis, University of New Mexico, 2015.

[21] R. Tielen, "High-order material point method," Master's thesis, Delft University of Technology, 2016.

[22] R. Tielen, E. Wobbes, M. Möller, and L. Beuth, "A high order material point method," *Procedia Engineering*, vol. 175, pp. 265–272, 2017.

[23] J. S. Hesthaven and T. Warburton, *Nodal discontinuous Galerkin methods: algorithms, analysis, and applications*. Springer Science & Business Media, 2007.

[24] S. K. Godunov, "Finite difference method for numerical computation of discontinuous solutions of the equations of fluid dynamics," *Mathematicheskii Sbornik*, vol. 47(89), no. 3, pp. 271–306, 1959.

[25] E. F. Toro, *Riemann solvers and numerical methods for fluid dynamics: a practical introduction*. Springer Science & Business Media, 2013.

[26] C. Hirsch, *Numerical computation of internal and external flows: The fundamentals of computational fluid dynamics*. Butterworth-Heinemann, 2007.

[27] G. Chavent and G. Salzano, "A finite-element method for the 1-d water flooding problem with gravity," *Journal of Computational Physics*, vol. 45, no. 3, pp. 307–344, 1982.

[28] G. Chavent and B. Cockburn, "The local projection $P^0 - P^1$-discontinuous-galerkin finite element method for scalar conservation laws," *ESAIM: Mathematical Modelling and Numerical Analysis*, vol. 23, no. 4, pp. 565–592, 1989.

[29] B. Cockburn and C.-W. Shu, "Runge-kutta discontinuous galerkin methods for convection-dominated problems," *Journal of Scientific Computing*, vol. 16, pp. 173–261, Sep 2001.

[30] R. Abedi, B. Petracovici, and R. Haber, "A space-time discontinuous galerkin method for linearized elastodynamics with element-wise momentum balance," *Computer Methods in Applied Mechanics and Engineering*, vol. 195, no. 25, pp. 3247–3273, 2006. Discontinuous Galerkin Methods.

[31] G. Yong, S. Zheng, C. Zhen, Z. Xiong, and L. Yu, "Enhancement of the material point method using b-spline basis functions," *International Journal for Numerical Methods in Engineering*, vol. 113, no. 3, pp. 411–431, 2018.