G.-H Cottet

J.-M Etancelin

F Perignon

C Picard

Mathematical Modelling and Numerical Analysis Will be set by the publisher Modélisation Mathématique et Analyse Numérique HIGH ORDER SEMI-LAGRANGIAN PARTICLE METHODS FOR TRANSPORT EQUATIONS: NUMERICAL ANALYSIS AND IMPLEMENTATION ISSUES

Keywords: AMS Subject Classification. -Please, give AMS classification codes - Advection equations, particle methods, semi-lagrangian methods, GPU computing

This paper is devoted to the definition, analysis and implementation of semi-Lagrangian methods as they result from particle methods combined with remeshing. We give a complete consistency analysis of these methods, based on the regularity and momentum properties of the remeshing kernels, and a stability analysis of a large class of second and fourth order methods. This analysis is supplemented by numerical illustrations. We also describe a general approach to implement these methods in the context of hybrid computing and investigate their performance on GPU processors as a function of their order of accuracy.

Introduction

Particle methods are Lagrangian methods that have been designed for advection dominated problems, with applications mostly in plasma physics [START_REF] Hockney | Simulation Using Particles[END_REF], incompressible flows [START_REF] Chorin | Numerical study of slightly viscous flow[END_REF][START_REF] Cottet | Vortex methods[END_REF][START_REF] Leonard | Computing three-dimensional incompressible flows with vortex elements[END_REF], or gas dynamics [START_REF] Monaghan | An introduction to sph[END_REF]. Following [START_REF] Koumoutsakos | Inviscid axisymmetrization of an elliptical vortex[END_REF][START_REF] Koumoutsakos | High resolution simulation of the flow around an impulsively started cylinder using vortex methods[END_REF], particle methods are often associated with remeshing, in order to maintain the regularity of the particle distribution and, as a consequence, to control the accuracy of the method in presence of strong strain in the flow carrying the particles. Remeshing removes the grid free nature of particle methods but maintains some of its features like conservation, locality, and stability. It not only guarantees accuracy but also allows to combine in a seamless fashion particle methods with grid-based techniques, like domain decomposition methods and fast FFT-based Poisson solvers for field evaluations [START_REF] Cocle | Combining the vortex-in-cell and parallel fast multipole methods for efficient domain decomposition simulations[END_REF][START_REF] Ould-Salihi | Blending finite-difference and vortex methods for incompressible flow computations[END_REF], multi-resolution and adaptive mesh refinement [START_REF] Bergdorf | Multilevel adaptive particle methods for convection-diffusion equations[END_REF][START_REF] Bergdorf | A lagrangian particle-wavelet method[END_REF].

In practice particle remeshing occurs every few time steps. In that case, and as long as the remeshing frequency does not increase when the discretization parameters tend to zero, Remeshed Particle Methods (RPM in the sequel) can be viewed as conservative forward semi-lagrangian methods. Classical semilagrangian methods for transport equations combine tracking of trajectories originating at grid points and interpolation from the grid. Because they operate on local point values, semi-lagrangian methods in general do not conserve mass, unless some specific treatment is done. By contrast, particle methods transport masses which makes them inherently conservative.

Forward Semi Lagrangian methods have already been designed and used with success in plasma physics [START_REF] Crouseilles | A forward semi-lagrangian method for the numerical solution of the vlasov equation[END_REF] and geophysical flows [START_REF] Cotter | The remapped particle-mesh semi-lagrangian advection scheme[END_REF]. Semi Lagrangian particles, as they come out of RPM, have some distinctive features. First, the projection on the grid is explicit, based on explicit remeshing kernels, and does not require to solve a linear system to recover grid values. The computational effort can thus be restricted to the support of the solution and efficient parallel algorithms can be implemented. Moreover, remeshed particle methods can be derived at any order in a systematic fashion, whereas the methods derived in [START_REF] Cotter | The remapped particle-mesh semi-lagrangian advection scheme[END_REF][START_REF] Crouseilles | A forward semi-lagrangian method for the numerical solution of the vlasov equation[END_REF] seem to be restricted to first order in space [START_REF] Respaud | Analysis of a new class of forward semi-lagrangian schemes for the 1d vlasov-poisson equations[END_REF]. One goal of this paper is indeed to derive remeshed particle methods of arbitrary order, independently of any CFL condition, and to give a complete numerical analysis of a large class of second and fourth order methods.

RPM only involve local operations and are therefore well adapted to parallel implementations. In [START_REF] Rossinelli | GPU accelerated simulations of bluff body flows using vortex methods[END_REF] an implementation of remeshed particle methods on GPU processors was described for the two-dimensional Navier-Stokes equations with penalization to account for solid obstacles. In the present paper, following [START_REF] Magni | Accurate, non-oscillatory, remeshing schemes for particle methods[END_REF], we consider implementations of two or three-dimensional RPM combined with directional splitting, that is where particles are pushed and remeshed successively in the three directions. Directional splitting allows to use high order remeshing kernels, with large stencils, at a minimal cost. We in particular discuss how directional splitting impacts the memory management to optimize implementations of RPM on GPU.

The outline of this paper is as follows. In Section 2 we give the definition of RPM, recall how remeshing kernels are derived and produce a list of high order kernels. In Section 3 we give the numerical analysis of RPM, based on the regularity and moment properties of the remeshing kernels. In Section 4 we present numerical illustrations and refinement studies of RPM in one to three dimensions. Section 5 is devoted to our implementation of RPM on GPU, in the context of a software library devoted to hybrid computations where different solvers can be implemented on different platforms. Finally Section 6 is devoted to concluding remarks and we give in the appendix the analytical formulas of the remeshing kernels that are considered in the paper.

Remeshed particle methods

In all the sequel we consider advection equations in the following conservation form

∂u ∂t + div (au) = 0, x ∈ R d , t > 0, (2.1)
where a is a given smooth velocity field.

Definition

In this paper we consider the case of uniform grids (we refer to [START_REF] Bergdorf | Multilevel adaptive particle methods for convection-diffusion equations[END_REF][START_REF] Bergdorf | A lagrangian particle-wavelet method[END_REF] for the design and applications of adaptive RPM with variable grid-size). RPM consist of alternating particle motion and remeshing. Particles are at the beginning of every time-step on a regular grid x i = i∆x, i ∈ Z d , then move with the following equation

x n+1 i = x i + ãn i ∆t.
(2.2) In the above equation ãn i denotes an evaluation of the velocity field at time t n = n∆t and location x i which depends on the chosen time-stepping scheme. Remeshing follows, through interpolation with a remeshing kernel Γ which satisfies Γ(-x) = Γ(x). If u n i denotes the value approximating u(x i , t n ), this gives the following formula :

u n+1 i = j u n j Γ x n+1 j -x i ∆x , i ∈ Z d , n 0. (2.3)
In traditional remeshing schemes for multidimensional problems, multidimensional kernels Γ are derived as tensor products of 1D formulas. In the present study we follow the approach initiated in [START_REF] Magni | Accurate, non-oscillatory, remeshing schemes for particle methods[END_REF] which consists of reducing the advection problem (2.1) into one-dimensional advection equations through directional splitting. One the one hand, this approach has the drawback that increasing the order in time of the method is less straightforward than for the original multidimensional problem. On the other hand, it offers the advantage to significantly reduce the computational cost of the method if high order kernels, with large stencils, are used. In short, if the kernel Γ involves a stencil with N s grid points, in three dimensions a multidimensional remeshing method will cost O(N 3 s ) per particle, while a directional splitting will cost O(3N s ) operations. For values of N s 6, a typical case we will consider in the sequel, the ratio is larger than 10. Directional splitting also allows to derive limiting techniques to reduce oscillations [START_REF] Magni | Accurate, non-oscillatory, remeshing schemes for particle methods[END_REF] and, for the numerical analysis of the method, one only needs to consider the case d = 1.

Derivation of remeshing kernels

We consider the equation (2.1) and its approximation formulas (2.2), (2.3) for d = 1. As will be clear from the numerical analysis below, one key feature which controls the accuracy of RPM is the conservation of momentum up to a given order p. More precisely, one seeks to satisfy for the scheme (2.3) the following identity

i u n+1 i x α i = i u n i x α i , 0 ≤ α ≤ p, (2.4)
for all sequences (u n i ). It is readily seen that this is equivalent to the following moments conditions for the remeshing kernel

Γ k∈Z k α Γ(x -k) = x α , 0 ≤ α ≤ p, x ∈ R (2.5)
or, equivalently,

k∈Z (x -k) α Γ(x -k) = 1 if α = 0 0 if 1 ≤ α ≤ p , x ∈ R. (2.6)
Note that for α = 0 these conditions enforce the conservation of mass. Using these identities for a given value of p and assuming that the kernel Γ remeshes particle weights among the p + 1 nearest grid points, one obtains a piecewise polynomial function of degree p. For p = 1 this gives the piecewise linear tent function. With p = 2 one obtains a piecewise quadratic function, the so-called Λ 2 formula, that has been used with success in the first particle simulation of the Navier-Stokes equations using remeshing in a systematic fashion [START_REF] Koumoutsakos | High resolution simulation of the flow around an impulsively started cylinder using vortex methods[END_REF].

This method to derive kernels is straightforward but has the drawback that it does not deliver smooth kernels. The kernel Λ 2 is not even continuous (this is the case more generally for kernels corresponding to even values of p; for odd values, the kernels are continuous but their derivatives are discontinuous). This lack of smoothness may result in a loss of accuracy or in oscillations, in particular if large time-steps are used. In [START_REF] Magni | Accurate, non-oscillatory, remeshing schemes for particle methods[END_REF] local correction techniques where derived to guarantee at least first order for the kernel Λ 2 and third order for the analogous kernel Λ 4 corresponding to p = 4.

To derive smooth kernels, one option is to use extrapolation techniques starting from smooth B-splines. If we denote by M 1 the top-hat filter with support in [-1/2, +1/2], and by M n its successive convolution:

M n = M ( * n) 1 ∈ W n,∞ (R), the idea is to derive linear combinations of M n , xM ′ n , x 2 M ′′ n ,.
. . to cancel the successive continuous moments of Γ. More precisely, given an even integer p ≥ 2 and k > 1 + p/2, one looks for coefficients α 1 , . . . , α p/2+1 such that Γ =

p/2 l=0 α l+1 x l M (l) k satisfies y α Γ(y) dy = 1 if α = 0, 0 if 1 ≤ α ≤ p. (2.7)
For symmetry reasons these conditions need only be enforced for even values of α. This leads to a square linear system of size 1 + p/2 the coefficients of which only involve the even moments of M k . With this method, one for instance readily obtains the following kernels

Γ = 1 2 (3M 4 + xM ′ 4 ), (2.8)
and

Γ = 1 8 (15M 8 + 9xM ′ 8 + x 2 M ′′ 8 ).
(2.9)

The first formula corresponds to a kernel of class C 1 , piecewise cubic, with a support of size 4 and p = 2. It was derived under the name of M ′ 4 in [START_REF] Monaghan | Extrapolating b splines for interpolation[END_REF] and has been and still is extensively used in particle simulations, in particular of vortex flows. The second formula corresponds to a kernel of class C 4 , piecewise polynomial of degree 7, with a support of size 8 and p = 4.

The moment conditions (2.7) concern continuous moments of the kernel. To check that the discrete moment conditions (2.6) are satisfied as well, one can use an equivalent condition using the Fourier transform of the kernel [START_REF] Cottet | Vortex methods[END_REF][START_REF] Schoenberg | Contribution to the problem of approximation of equidistant data by analytic functions[END_REF]. If Γ(ξ) = Γ(y)e -iξy dy, it can be shown ( [START_REF] Cottet | Vortex methods[END_REF]) that the properties (2.5) are satisfied provided Γ fulfills the following conditions: Γ(ξ) -1 has a zero of order p + 1 at ξ = 0 Γ(ξ) has a zero of order p + 1 at all ξ = 2πm, m = 0.

Since Γ (α) is proportional to y α Γ(y), the first condition above is clearly equivalent to the conditions (2.7). Moreover one has

M 1 (ξ) = sin (ξ/2) ξ/2 and thus M n (ξ) = sin (ξ/2) ξ/2 n .
As a result, the Fourier transform of the functions

x l M (l)
k has a zero of order k -l at all non zero multiple of π. From these observations, one can easily check that the kernels (2.8), and (2.9) do satisfy the discrete moment conditions (2.5) with p = 2 and p = 4 respectively.

The approach just presented leads to smooth and high order kernels. However the kernels derived in this fashion do not necessarily satisfy the following interpolation property

Γ(i) = 1 if i = 0, 0 otherwise. (2.10)
Indeed the so-called kernel M ′ 4 given in (2.8) does satisfy this condition, but not the one derived in (2.9) from M 8 . Property (2.10) is natural as it ensures that, if the velocity is zero, the exact solution is algebraically conserved. Although this property does not enter the numerical analysis that follows, in practice it seems to have some importance, in particular to represent accurately the smallest scales in turbulent flows.

One can derive kernels which satisfy simultaneously and to any given order, regularity, moment conditions and the interpolation property. For a sake of simplicity, in the following we restrict ourselves to kernels with a stencil covering an even number 2M s of grid points. We seek kernels Γ that have the following properties:

P1: Γ has support in [-M s , +M s ],
P2: Γ is even and piecewise polynomial of degree M in intervals of the form [i, i + 1], P3: Γ is of class C r , P4: Γ satisfies the moment properties (2.5) for a given value of p, P5: Γ satisfies the interpolation property (2.10).

Such kernels are determined by M s (M + 1) coefficients. The regularity property P3 imposes M s (r + 1) interface conditions at integer values, and [(r + 1)/2] conditions to express that derivatives of odd order vanish at zero. The properties P4 and P5 impose p + 1 + M s conditions. One reasonable constraint under which one can expect to find kernels satisfying these conditions is therefore

M s (M + 1) ≥ (r + 1)M s + [(r + 1)/2] + p + 1 + M s .
The Table 1 lists several kernels that have been obtained through this approach by symbolic calculations. In this table, the kernels have been labelled by 2 indices that refer to the regularity and the order to which moment conditions are satisfied, as we will see that these are the parameters which control the order of accuracy of the RPM: Λ p,r is a kernel in W r,∞ (R) which satisfies (2.5). Λ 4,2 corresponds to the kernel M ′ 4 already mentioned. The kernel Λ 4,2 was derived with this approach and used for the first time in [START_REF] Bergdorf | A lagrangian particle-wavelet method[END_REF] under the name of M ′ 6 . For a sake of completeness, we have provided in the appendix the analytical formulas Moments (p in 2.5) Regularity Nb of grid points in stencil Degree Support

Λ 2,1 2 C 1 4 3 [-2; 2] Λ 2,2 2 C 2 4 5 [-2; 2] Λ 2,3 2 C 3 4

Numerical Analysis

We consider in this section RPM with kernels satisfying the moment properties (2.5) and the following regularity conditions :

Γ ∈ W r+1,∞ (R) and Γ ∈ C ∞ (]l, l + 1[) , l ∈ Z. (3.1)
The RPM is defined by the formulas (2.2), (2.3) and we will denote by T i u(•, t n ) the result of the scheme (2.3), at the grid point x i , starting from grid values u(x j , t n ).

In this section we are interested by the stability and spatial accuracy of the method. For a sake of simplicity we will therefore assume that a does not depend on time and that particles advance with an explicit first-order Euler scheme. In this case we simply have ãn j = a(x j ). A striking feature of RPM, common with all semi-lagrangian methods, is that their stability does not rely on CFL conditions. In this section we prove stability and consistency results under the condition

∆t ≤ M a ′ L ∞ (3.2)
for a given constant M < 1. This condition in particular ensures that particle trajectories cannot intersect. The constant M is often called Lagrangian CFL number.

Consistency

We will prove the following consistency result Proposition 1. Assume that the condition (3.2) is satisfied, and that the moment and regularity conditions (2.5), (3.1) hold for some r, p > 1. Let T > 0 and assume further that a and the solution u to equation (2.1) belong to L ∞ 0, T ; W r+1,∞ (R) . Then, if we set β = inf (r, p), the following estimate holds

u x i , t n+1 = T i (u (•, t n )) + O ∆t 2 + O ∆t∆x β + ∆x β+1 + ∆t β+1 . (3.3)
Moreover if every cell of size ∆x contains exactly one particle after an advection step, then β = p.

Proof. If we use the notation u n j for u(x j , t n ) and set j = i + k we can write the following expansion

u n j = β l=0 k l ∆x l l! u (l) (x i ) + O ∆x β+1 , (3.4)
where u (l) (x i ) stands for ∂ l u/∂x l (x i , t n ). Next, one can write x n+1 j = x i + k∆x + a i+k ∆t and therefore

Γ x n+1 j -x i ∆x = Γ(k + λ i+k ),
where λ j is the local CFL number: λ j = a j ∆t/∆x. We then write λ j = λ i + (a j -a i )∆t/∆x and use the following formula

(f • g) (m) (x) = m |q|=1 c q f q 1 +•••+qm (g(x)) m s=1 g (s) (x)
qs where q = (q 1 , . . . , q m ) is a multi-index, |q| = q 1 + 2q 2 + • • • + mq m and c q are positive coefficients, to obtain

Γ x n+1 j -x i ∆x = β m=0 (k∆x) m m |q|=1 c q Γ (q 1 +•••+qm) (k + λ i ) ν q 1 +•••+qm m s=1 a (s) (x) qs (3.5) + O   (k∆x) β+1 Γ β+1,∞ a β+1,∞ |q|=β+1 ν q 1 +•••+q β+1   .
In the above equation we have denoted by ν the ratio ∆t/∆x. Since q 1 + • • • + q β+1 ≤ |q| we have

|q|=β+1 ν q 1 +•••+q β+1 ≤ C ν + ν β+1 .
The remainder in (3.5) can thus be bounded by O ∆x β+1 + ∆x β ∆t + ∆t β+1 . Combining (3.4) and (3.5) we obtain

T i (u(•, t n )) = j u n j Γ x n+1 j -x i ∆x = (3.6) k 0≤l+m≤β (k∆x) l+m u (l) (x i ) l! m |q|=1 c q Γ (q 1 +•••+qm) (k + λ i ) ν q 1 +•••+qm m s=1 a (s) (x i ) qs + O ∆x β+1 + ∆x β ∆t + ∆t β+1 = k E k (x i ) + O ∆x β+1 + ∆x β ∆t + ∆t β+1 .
By differentiation, the moment conditions (2.5) yield

k k q ′ Γ (q) (k + x) = (-x) q ′ -q q ′ (q ′ -1) • • • (q ′ -q + 1), for 0 ≤ q ≤ q ′ ≤ β.
Using these identities with x = λ i and observing that λ

l+m-(q 1 +•••+qm) i ν q 1 +•••+qm = d i (∆t/
E k (x i ) = 0≤l+m≤β O ∆t l+m + O ∆x β+1 + ∆x β ∆t + ∆t β+1 . (3.7)
One can easily check that the zero and first order terms in the above expansion, corresponding to

0 ≤ l + m ≤ 1, result in u(x i , t n ) -∆t a(x i )∂u/∂x(x i , t n ) -∆t a ′ (x i )u(x i , t n ) = u(x i , t n ) -∆t ∂(a u)/∂x(x i , t n ).
We therefore finally obtain

T i (u(•, t n )) = u n i -∆t ∂(au) ∂x (x i , t n ) + O ∆t 2 + O ∆x β+1 + O ∆x β ∆t (3.8) = u(x i , t n+1 ) + O ∆t 2 + O ∆x β+1 + ∆x β ∆t + ∆t β+1 . (3.9)
To prove the second assertion of our proposition we observe that if, at the end of an advection step, every cell within the stencil centered at the grid point x i contains exactly one particle, since particles cannot cross this means that λ j and λ i lie in the same interval between successive integers. The kernel Γ is C ∞ in such intervals and thus the expansion 3.5 is valid to any order. The coefficient β can therefore be taken equal to p.

Stability

We start with the linear stability analysis and assume that the velocity field a is constant.

Linear stability

In this paragraph, we will show unconditional stability (with respect to λ = a ∆t/∆x) for a certain class of remeshing kernels. We set, for k ∈ Z,

α k (λ) = Γ(k + λ) , A k (λ) = i α i (λ)α i+k (λ). (3.10)
In all the calculations that follow we will write α k and A k for α k (λ) and A k (λ), respectively, when there is no ambiguity on the value of λ. We start with the following result Lemma 2. If the kernel Γ satisfies (2.5) then

k,l (k -l) q α k α l = 0 , k≥1 k q A k = 0, if 1 ≤ q ≤ p, q even. (3.11)
Proof. We have

k,l (k -l) q α k α l = k,l q m=0 C m q k m (-l) q-m α k α l .
By the moment properties (2.5)

k,l k m (-l) q-m α k α l = (-λ) m λ q-m and thus k,l (k -l) q α k α l = q m=0 C m q (-λ) m λ q-m = 0.
The second identity readily follows, since

A -k = A k .
We consider the case of kernels involving stencils with an even number of points. We further restrict our analysis to second order kernels involving 4 points (which means that the support of Γ is the interval [-2, +2]) and fourth order kernels involving 6 points (the support of Γ is the interval [-3, +3]). The key properties of the remeshing kernel that are needed to ensure stability are the moment conditions and a sign and a decay properties for the kernel values. More precisely we will prove the following result Proposition 3. Assume that the velocity field a is constant, and that the kernel Γ satisfies one of the following conditions:

i: Γ satisfies (2.5) with p = 2, Γ has support in [-2, +2] and satisfies for all values of λ ∈ [0, 1]

A 1 ≥ 0 , A 2 ≤ 0 , A 3 ≥ 0 , -A 2 ≥ 6A 3 .
(3.12)

ii: Γ satisfies (2.5) with p = 4, Γ has support in [-3, +3] and satisfies for all values of λ ∈ [0, 1]

A 1 ≥ 0 , A 2 ≤ 0 , A 3 ≥ 0 , A 4 ≤ 0 , A 5 ≥ 0 , A 3 ≥ -8A 4 . (3.13)
Then the remeshed particle method is unconditionally stable.

Before proceeding with the proof of this result, some comments are in order. For all the kernels that are used in practice one has Γ(x) ≥ 0 for x ∈ [-1, 1] and Γ alternates sign in successive integer intervals. The conditions (3.12), (3.13) therefore mean that the kernel decays fast enough. In the case i,

A 2 = α -2 α 0 + α -1 α 1 and A 3 = α -1 α 1 . The condition -A 2 ≥ 6A 3 is thus satisfied as soon has α 0 ≥ 6α 1 .
For the Λ 2,1 kernel, this condition can easily be checked: one has, for λ ∈ [0, 1]

|α 0 /α 1 | = 2 + 2λ -3λ 2 (1 -λ)λ = 2 + 2 -λ 2 (1 -λ)λ ≥ 2 + 1 (1 -λ)λ ≥ 6.
Similar calculations show that the other 2nd order kernels in the Table 1 satisfy this decay condition. For 4th order, 6 points kernels, the analytic calculations are more involved, but it is possible to check by symbolic calculations that the conditions (3.13) hold for the kernels that are indicated in Table 1. Our result only covers 2nd and 4th order kernels. However one can conjecture from the proof that will be given below that it extends to higher order kernels, with larger stencils, under an appropriate decay condition for the coefficients A k .

We now proceed with the proof of Proposition 3 and begin with the case i.

We define by I the integer such that λ = a∆t/∆x ∈ [I, I + 1[ and we set µ = λ -I ∈ [0, 1[. We write j = i -I + k and set v i = u n i-I . The scheme (2.3) now reads

u n+1 i = k v i+k Γ(k + µ). (3.14) Since the support of Γ is in [-2, +2], the coefficients α k (µ) vanish if k ≤ -3 or k ≥ 2. Using (3.10) we obtain i |u n+1 i | 2 = i -2≤k,l≤1 v i+k v i+l α k (µ) α l (µ) = i v 2 i k α 2 k + 2 i -2≤k<l≤1 v i+k v i+l α k α l . (3.15)
We then write

2v i+k v i+l = v 2 i+k + v 2 i+l -|v i+k -v i+l | 2 .
The conservation of the first moment gives

1 = k α k 2 = k α 2 k + 2 -2≤k<l≤1 α k α l ,
and therefore

i u n+1 i 2 = i v 2 i - i S i (3.16)
where

S i = -2≤k<l≤1 α k α l |v i+k -v i+l | 2 . Since obviously i v 2 i = i |u n i | 2 , it remains to prove that i S i ≥ 0.
Using the change of index l = k + m for l > k and the notations in (3.10) allows to write

i S i = k 1≤m≤3 α k α k+m i |v i+m -v i | 2 = 1≤m≤3 A m i |v i+m -v i | 2 .
We set δ i = v i+1 -v i to rewrite the above expression as

i S i = A 1 i δ 2 i + A 2 i (δ i + δ i+1 ) 2 + A 3 i (δ i + δ i+1 + δ i+2 ) 2 .
(3.17)

Expanding the above sums we obtain

i S i = (A 1 + 2A 2 + 3A 3 ) i δ 2 i + (2A 2 + 4A 3 ) i δ i δ i+1 + 2A 3 i δ i δ i+2 .
We again write

2δ i δ i+1 = δ 2 i + δ 2 i+1 -|δ i -δ i+1 | 2 and a similar identity for 2δ i δ i+2 to obtain i S i = (A 1 + 4A 2 + 9A 3 ) i δ 2 i -(A 2 + 2A 3 ) i |δ i -δ i+1 | 2 -A 3 i |δ i -δ i+2 | 2 .
By (3.11) with q = 2 we get

A 1 + 4A 2 + 9A 3 = 0. If we set η i = δ i -δ i+1 we thus have i S i = -(A 2 + 2A 3 ) i |η i | 2 -A 3 i |η i + η i+1 | 2 . Since A 3 > 0 and |η i + η i+1 | 2 ≤ 4(η 2 i + η 2 i+1 ) we can write i S i = -(A 2 + 6A 3 ) i |η i | 2 ≥ 0,
provided the decay property (3.12) is satisfied.

We now turn to the second assertion of the proposition. We start from (3.15) and obtain an estimate similar to (3.17), with m = 5 instead of 3 since the kernel has now a support of size 6, yields:

i S i = k 1≤m≤3 α k α k+m i |v i+m -v i | 2 = 1≤m≤5 A m i |v i+m -v i | 2 (3.18)
We set δ i = v i+1 -v i and obtain:

i S i = A 1 i δ 2 i + A 2 i (δ i + δ i+1 ) 2 + A 3 i (δ i + δ i+1 + δ i+2 ) 2 +A 4 i (δ i + δ i+1 + δ i+2 + δ i+3 ) 2 + A 5 i (δ i + δ i+1 + δ i+2 + δ i+3 + δ i+4 ) 2 = (A 1 + 2A 2 + 3A 3 + 4A 4 + 5A 5 ) i δ 2 i + (2A 2 + 4A 3 + 6A 4 + 8A 5 ) i δ i δ i+1 +(2A 3 + 4A 4 + 6A 5 ) i δ i δ i+2 + (2A 4 + 4A 5 ) i δ i δ i+3 + 2A 5 i δ i δ i+4 .
Rewriting double products as above and introducing

η i = δ i -δ i+1 , we obtain i S i = (A 1 + 4A 2 + 9A 3 + 16A 4 + 25A 5 ) i δ 2 i -(A 2 + 2A 3 + 3A 4 + 4A 5 ) i η 2 i -(A 3 + 2A 4 + 3A 5 ) i (η i + η i+1 ) 2 -(A 4 + 2A 5 ) i (η i + η i+1 + η i+2 ) 2 -A 5 i (η i + η i+1 + η i+2 + η i+3 ) 2 .
By (3.11) with q = 2 we have A 1 + 4A 2 + 9A 3 + 16A 4 + 25A 5 = 0. Expanding the remaining squares in the above identity we get

i S i = -(A 2 + 4A 3 + 10A 4 + 20A 5 ) i η 2 i -2(A 3 + 4A 4 + 10A 5 ) i η i η i+1 -2(A 4 + 2A 5 ) η i η i+2 -2A 5 η i η i+3 .
We again write

2η i η i+1 = η 2 i + η 2 i+1 -|η i -η i+1 | 2
and similar identities for the other double products to obtain

i S i = -(A 2 + 6A 3 + 20A 4 + 50A 5 ) i η 2 i + (A 3 + 4A 4 + 10A 5 ) i (η i -η i+1 ) 2 + (A 4 + 4A 5 ) i (η i -η i+2 ) 2 + A 5 i (η i -η i+3 ) 2
Subtracting (3.11) with q = 4 from the corresponding identity with q = 2 yields A 2 +6A 3 +20A 4 +50A 5 = 0. Moreover, since A 4 ≤ 0 and A5 ≥ 0 we can write

(A 4 + 4A 5 ) i (η i -η i+2 ) 2 ≥ A 4 i (η i -η i+2 ) 2 ≥ 4A 4 i (η i -η i+i ) 2 . Therefore i S i ≥ (A 3 + 8A 4 ) i (η i -η i+i ) 2 ≥ 0,
provided the decay property (3.13) is satisfied. This concludes our proof in the case of a constant velocity field.

General case

We now consider the general case of a smooth, non-constant, velocity field a. Proposition 4. Assume that the velocity field a is in W 1,∞ (R) and that the assumptions of Proposition 3 are satisfied. Assume in addition that the kernel Γ satisfies the interpolation property (2.10). Then there exists a constant C, independent of ∆t and ∆x, such that

i u n+1 i 2 ≤ (1 + C∆t) i |u n i | 2 (3.19)
Proof. We will only give the proof in the case of a second order method corresponding to the case i of Proposition 3. The proof in the case of ii follows along the same lines.

The local Courant number λ can now vary from one particle to the next. To start with we will first consider the case when this number remains in the same integer interval, that is we make the assumption that

∃ I ∈ Z such that, ∀i ∈ Z, λ i = a(x i , t n )∆t/∆x ∈ [I, I + 1].
We proceed like in the proof of Proposition 3. We set µ j = λ j -I ∈ [0, 1[, v i = u n i-I and j = i -I + k, and the scheme (3.14) becomes

u n+1 i = k v i+k Γ(k + µ j ) = k v i+k α k (µ j ). (3.20)
We first observe that, by the regularity of the velocity field a,

|λ j -λ j ′ | = O |j -j ′ |∆t and, since Γ ∈ W 1,∞ (R), for -2 ≤ k < l ≤ 1, j = i -I + k, j ′ = i -I + l α k (µ j ) = α k (µ j ′ ) + O(∆t). (3.21) In particular α k (µ j ) = α k (µ i-I ) + O(∆t) if j = i -I + k with k ∈ [-2, 1]. This allows to write i |u n+1 i | 2 = i   -2≤k≤1 v i+k α k (µ i-I )   2 + O ∆t|v| 2 ,
where we have used the notation

|v| 2 = i |v i | 2 .
Proceeding like in estimates (3.16), (3.17) we can write :

i |u n+1 i | 2 = O ∆t|v| 2 + i k |v i+k | 2 α 2 k (µ i-I ) + k<l |v i+k | 2 + |v i+l | 2 -|v i+k -v i+l | 2 α k (µ i-I )α l (µ i-I ) . (3.22)
Using again (3.21) and setting j = i + k we have

i k |v i+k | 2 α 2 k (µ i-I ) = j k |v j | 2 α 2 k (µ j-I ) + O ∆t|v| 2 , i k<l |v i+k | 2 α k (µ i-I )α l (µ i-I ) = j k<l |v j | 2 α k (µ j-I )α l (µ j-I ) + O ∆t|v| 2 , i k<l |v i+l | 2 α k (µ i-I )α l (µ i-I ) = j k<l |v j | 2 α k (µ j-I )α l (µ j-I ) + O ∆t|v| 2 , i k<l |v i+k -v i+l | 2 α k (µ i-I )α l (µ i-I ) = j m A m (µ j-I )|v j -v j+m | 2 + O ∆t|v| 2 .
This allows us to rewrite the summation in the right hand side of (3.22) as

j |v j | 2 k α 2 k (µ j-I ) + 2 k<l α k (µ j-I )α l (µ j-I ) - j m A m (µ j-I )|v j -v j+m | 2 + O ∆t|v| 2 .
By the conservation of the first moment we have

k α 2 k (µ j-I ) + 2 k<l α k (µ j-I )α l (µ j-I ) = 1
and, like for the constant velocity case, it only remains to check the sign of

S = j m A m (µ j-I )|v j -v j+m | 2 .
We proceed like in the constant velocity case and set δ i = v i+1 -v i and η i = δ i+1 -δ i . We expand the squares in S using the fact that

δ 2 i+1 A 2 (µ i-I ) = δ 2 i+1 A 2 (µ i+1-I ) + O(∆t) |v i+2 | 2 + |v i+1 | 2
and similar expressions for δ 2 i+1 A 3 (µ i-I ) and δ 2 i+2 A 3 (µ i-I ). This leads to

S = |v| 2 O(∆t) + i δ 2 i [A 1 (µ i-I ) + 4A 2 (µ i-I ) + 16A 3 (µ i-I )] - i η 2 i [A 2 (µ i-I ) + 2A 3 (µ i-I )] - i (η i + η i+1 ) 2 A 3 (µ i-I ). (3.23)
We next write

(η i + η i+1 ) 2 A 3 (µ i-I ) ≤ 2 η 2 i A 3 (µ i-I ) + η 2 i+1 A 3 (µ i+1-I ) + C∆t |v i+3 | 2 + |v i+2 | 2 + |v i+1 | 2 .
Due to (3.11), under the condition (3.12), (3.23) finally yields

S ≥ -C∆t|v| 2 and therefore i |u n+1 i | 2 ≤ i |u n i | 2 + C∆t.
We finally turn to the general case, when neighboring particles can have their local CFL numbers in different integer intervals. We can always group particles in subsets where the CFL numbers are in the same integer interval. For I ∈ Z, we set

J I = { j ∈ Z such that λ j ∈ [I, I + 1[ }.
We can rewrite the remeshed particle method as

u n+1 i = I j∈J I u n j Γ(λ j + j -i) = I S I (i) (3.24)
The next step is to take the square of the above identity. The key point is to observe that, by the interpolation property (2.10), double products of the form S I (i)S I ′ (i) with I = I ′ are of order ∆t. We indeed observe that, due to the regularity of a, we have

|λ j -λ j ′ | ≤ C∆t|j -j ′ |
and therefore, if the kernel Γ has a support of size 2M s and if ∆t is small enough

Γ(λ j + j -i)Γ(λ j ′ + j ′ -i) = 0 ⇒ |j -j ′ | ≤ 2M s + 1. (3.25)
Next, again in view of the regularity of a, if two neighboring particles have indices which belong to a different set J I this implies that their local CFL number is close to an integer value. More precisely

j ∈ J I , j ′ ∈ J I ′ , I = I ′ , |j -j ′ | ≤ 2M s + 1 ⇒ ∃(K, K ′ ) ∈ Z 2 , K = K ′ such that j + λ j = K + O(∆t), j ′ + λ j ′ = K ′ + O(∆t).
By the interpolation property (2.10) and the regularity of Γ this implies

j ∈ J I , j ′ ∈ J I ′ , I = I ′ , |j -j ′ | ≤ 2M s + 1 ⇒ Γ(λ j + j -i)Γ(λ j ′ + j ′ -i) = O(∆t).
Combined with (3.25) this shows that, for

I = I ′ , S I (i)S I ′ (i) ≤ C∆t j∈J I J I ′ |u n j | 2 . (3.26)
It remains now to sum over i the above identities. To identify the indices j which, for a given index i, contribute to the sum S I (i), we denote by φ the application such that

φ(x) = x + a(x)∆t.
In view of (3.2) we have 0 < 1 -M ≤ φ ′ (x) ≤ 1 + M . φ is thus one-to-one and its inverse φ -1 is strictly increasing. If we set ψ(i) = φ -1 (i∆x)/∆x, the indices j which contribute to S I (i) must satisfy

i∆x -M s ∆x ≤ j∆x + a j ∆t ≤ i∆x + M s ∆x which yields ψ(i -M s ) ≤ j ≤ ψ(i + M s ).
If a particle j contributes to two different terms S I (i) and S I (i ′ ) with, say, i < i ′ this means that

ψ(i ′ -M s ) ≤ j ≤ ψ(i + M s ) and therefore |i ′ -i| ≤ M s .
As a result, a given particle of index j contributes to at most 2M s sums S I (i). One can thus deduce from (3.26) that, for a certain constant C independent of ∆t and ∆x

i I =I ′ S I (i)S I ′ (i) ≤ C∆t|u n | 2 ,
and, from (3.24),

i |u n+1 i | 2 ≤ i I S 2 I (i) + C∆t|u n | 2 . (3.27)
Finally, if we set

w I j = u n j if j ∈ J I 0
otherwise, by the preceding proof in the case when I takes a constant value, we have i

S 2 I (i) ≤ (1 + C∆t)|w I | 2
and

I i S 2 I (i) ≤ (1 + C∆t) I |w I | 2 .
In view of (3.27) this leads to

i |u n+1 i | 2 ≤ |u n | 2 + C∆t|u n | 2
and our proof is completed.

Remark. For CF L ≤ 1, Remeshed Particle Methods reduce to finite-difference schemes [START_REF] Cottet | Particle methods revisited: a class of high order finite-difference methods[END_REF]. The convergence analysis above contains thus as a by-product a convergence proof of a class of conservative second and fourth order finite-difference schemes.

Numerical illustrations

In this section we compare the accuracy of the various kernels considered above in 1D, 2D and 3D smooth flows. The reference [START_REF] Lagaert | Hybrid spectral particle method for turbulent transport of passive scalar[END_REF] deals with remeshed particle methods using the Λ 4,2 kernel in the context of the transport of passive scalar in turbulent flows. It in particular emphasizes the efficiency of this method in the case of high Schmidt numbers, that is when the diffusivity of the scalar is small compared to the viscosity of the flow. The extension of this work to higher order kernels introduced in the present paper will be reported in a future work.

We first consider a 1D advection equation in the interval [-1, +1] with periodic boundary condition. The velocity a and the initial condition u 0 are given by a(x) = 1 + sin (πx)/2 , u 0 (x) = sin (πx)

This flow induces compression and dilatation which successively increase and decrease the value of the solution. The exact solution is given by the following formula :

u(x(t), t) = u 0 (x 0 ) 2 + sin 2πx 0 2 + sin 2πx(t)
where x(t) is the trajectory with origin x 0 associated to the velocity field a, given by

arctan [tan (πx(t)) + 1/2] = arctan [tan (πx 0 ) + 1/2] + tπ √ 3/2.
The solutions are periodic in time with period T = 4/ √ 3.

Figure 1 shows the solution at initial time and t = √ 3. For this case we compare in Figure 2 the results obtained with the kernels Λ 2,1 , Λ 4,2 , Λ 2,2 and Λ 4,4 . In all cases, and throughout this section, particles where pushed with a fourth order Runge-Kutta time-stepping. The number of grid points ranged from 128 to 4096. For the purpose of this refinement study, we had to choose a constant value for the ratio ∆t/∆x. We fixed a CFL value of 12. This value corresponds, for the coarsest resolution considered in this test, to a constant M in (3.2) equal to 0.7. The error is measured in the maximum norm. The numerical analysis above indicates that the two last kernels should deliver respectively second and fourth order methods, while the two first kernels, because of their lack of regularity, should be limited to first and second order, respectively. One can indeed observe that, for a given number of conserved moments, increasing the regularity does improve the accuracy of the RPM. However it is interesting to note that the kernel Λ 2,1 already gives a second order method. The reason is that, in this example, zones where the solution undergoes large variations, and thus where numerical errors are likely to be larger, only rarely correspond to grid points where the local CFL numbers cross integer values.

To illustrate next the directional splitting which we use to address multidimensional problems, we first consider the classical two-dimensional case of a an off-centered circle strained in an incompressible rotational flow. More precisely, the computational box is the square [0, 1] 2 . The velocity field is given by

a(x 1 , x 2 , t) = f (t) -sin 2 (πx 1 ) sin(2πx 2 ), sin(2πx 1 ) sin 2 (πx 2 ) (4.2)
The initial condition is a sign function whose zero level set is a circle of radius 0.15 and centered at the point (0.5, 0.15). The function f (t) is introduced to allow the solution to eventually return to the initial condition. We take f (t) = cos(πt/12), and compare the computed solution at time t = 12 with the initial condition. In this experiment, like in the 3D case below, we use the classical second order Strang formula for the dimensional splitting, that is we successively push and remesh particles in the x 1 , x 2 , x 2 and x 1 directions for ∆t/2. It is indeed possible to derive and use higher order splitting methods, but our choice of a second order splitting was made for a practical reason -fourth order splitting methods are more demanding in terms of CPU and memory requirement. Furthermore, although second order splitting limits the theoretical order of convergence to 2, one can still see the benefit of using higher order kernels. Using a 6th order kernel in this example indeed yields an effective order of convergence close to 6. In this experiment the number of grid points in each direction ranged form 32 to 512. Like in the previous experiment, for the purpose of this refinement study, the CFL number was kept constant, equal to 12.

The average order of convergence is shown on the right picture of Figure 3. In this experiment the kernels Λ 2,2 , Λ 4,4 and Λ 6,6 would give similar results to the kernels Λ 2,1 , Λ 4,2 and Λ 6,4 respectively. Figure 4 shows the solution obtained with N = 256 grid points in each direction and the kernels Λ 2,1 and Λ 6,4 , at t = 6 and t = 12, compared to the exact solution. At t = 6, which is the time of maximal stretching for this experiment, the exact solution was obtained by a front tracking method using 10000 markers. This study clearly shows the gain obtained in his case by using high order kernels, even though the theoretical order of convergence is limited by that of the dimensional splitting.

As already stressed, one distinctive feature of Semi Lagrangian Particles are that their stability and the spatial order of convergence is not constrained by a CFL condition. To illustrate this property, we repeated the above experiment wit 256 points in each direction and a CFL number equal to 30. These grid resolution and time-step give a value of the constant M in (3.2) equal to 0.35. Figure 5 shows the zero level set at times 6 and 12. Although one cannot expect the same overall accuracy as in the previous experiment with such a large time-step, the interface appears to remain rather well captured.

We now turn to a 2D example where the kernel regularity seems to be more important than in the previous cases. We consider an expansion field given by

a(x 1 , x 2 ) = (x 1 /r, x 2 /r), r = x 2 1 + x 2 2 , (4.3)
and the initial value of u is an axisymmetric function supported in an annulus: The exact solution at time t is given by u(r, t) = u 0 (r -t)(r -t)/r. Figure 6 shows the result of a refinement study for the kernels Λ 2,1 , Λ 2,2 , Λ 4,2 and Λ 4,4 . In that case areas where the local CFL number cross integer values are aligned with the flow directions. As a result the kernel Λ 2,1 only gives first order convergence, while Λ 2,2 is close to the second order accuracy predicted by our analysis. Note that increasing the order and regularity of the kernel improves the convergence but does not allow to reach fourth order.

u 0 (r) = C[(r -r a )(r -r b )] 4 if r a ≤ r ≤ r b , 0
In this example the second order time-splitting limits the observed overall accuracy. We close this section with a 3D version case suggested in [START_REF] Leveque | High-resolution conservative algorithms for advection in incompressible flow[END_REF] and often used to validate level set methods.

In this example the computational box is the unit cube and the initial condition is a sphere of radius 0.15 centered around the point with coordinates (0.35, 0.35, 0.35). This sphere is advected with the velocity field a 1 (x 1 , x 2 , x 3 ) = 2 sin 2 (πx 1 ) sin(2πx 2 ) sin(2πx 3 ), (4.4) a 2 (x 1 , x 2 , x 3 ) = -sin(2πx 1 ) sin 2 (πx 2 ) sin(2πx 3 ), (4.5)

a 3 (x 1 , x 2 , x 3 ) = -sin(2πx 1 ) sin(2πx 2 ) sin 2 (πx 3 ). (4.6)
This motion wraps the sphere into thinner and thinner surfaces difficult to capture. In this case we want also to illustrate the fact that RPM only need particles in the support of the advected quantity. To capture the sphere we use a color function (1 inside the sphere, 0 outside). At the end of the remeshing step, particles are created only when their strength is above 0.001. In this experiment, we use 256 grid points in each direction. The CFL number is equal to 30, which corresponds to the same value of M as in the 2D case (4.2). The number of particles roughly varies between 1.5%, at t = 0, and 10%, at t = 4, of the total number of grid points inside the computational box. Due to the cut-off used in the remeshing step, the total mass was not algebraically conserved. However, in the present experiments the deviation did not exceed 0.1% of the total mass. Figure 7 shows the evolution of the volume inside the interface (this volume should remain constant at all times) obtained with the kernels Λ 2,1 , Λ 4,2 and Λ 8,4 . Figure 8 show the interface u = 0.5 at time t = 4, with the kernels Λ 2,1 and Λ 8,4 . These results confirm the gain in accuracy that can be obtained by using high order kernels. The next section will discuss the computational complexity associated with these kernels.

Remeshed Particle Methods and software engineering

One challenge in High Performance Computing is to develop algorithms that are able to take advantage of clusters made of emerging exascale hardware. These clusters are hybrid architectures which combine many heterogeneous distributed processors and accelerators. A natural idea is to match models, algorithms and physical scales resolved by these algorithms to the specific features of the processors involved in these architectures.

A typical example is the case of turbulent transport in incompressible flows. Hybrid algorithms can be designed to resolve in an optimal fashion the minimal scales present in the flow and in the advected scalar [START_REF] Lagaert | Hybrid spectral particle method for turbulent transport of passive scalar[END_REF]. Moreover in these hybrid approaches the different algorithms can have different parallel scalability and it may be advantageous to distribute these algorithms to different type of processors. The local nature of RPM make these methods well suited for highly parallel processors, like GPUs, and one may envision simulations of transport at high Schmidt numbers in turbulent flows (that is when the scalar smallest scales extend much beyond those of the flows) at very high resolution at a cost that does not exceed that of the flow solver at a much lower resolution.

To be able to implement hybrid algorithms on hybrid architectures, one needs to develop high level frameworks and libraries with a high level description which allows to distribute different solvers and grids to different parts of the clusters in a seamless fashion.

In this section we first outline the general approach followed for this framework. We then go into more details for a specific branch of this framework, namely the implementation on GPU of RPM. We detail the performance of the different kernels associated to the different parts of the algorithm and the overall performance of the GPU implementation of RPM. We in particular focus on the complexity of the RPM using various high order kernels that have been derived and analyzed in the previous section.

Library description

In this library we use object oriented programming techniques to reach a high level of modularity, with a strong focus on usability and flexibility. The goal is to enable the user to launch indifferently sequential, parallel or hybrid numerical simulations. We use Python as an abstraction framework.

Our library currently provides two levels of parallelism, MPI and OpenCL. However, thanks to the modularity it would also be possible to implement other parallelism paradigms such as task parallelism.

In order to achieve good portability, the computational frameworks are written using OpenCL C whenever it will enable good performances on the target architecture. OpenCL is an open standard for parallel programming of heterogeneous systems [START_REF] Munshi | The OpenCL Specification[END_REF]. It provides application programming interfaces to address hybrid platforms containing many CPUs and GPUs and a programming language based on C99 to write instructions executed concurrently on the OpenCL devices. This framework is used through the PyOpenCL python interface to OpenCL [START_REF] Klöckner | PyCUDA and PyOpenCL: A Scripting-Based Approach to GPU Run-Time Code Generation[END_REF].

In OpenCL applications, the program, that executes on the host system, must send OpenCL kernels to available devices in a command queue. The execution model of the application defines the devices and queues management. A given kernel contains an executable code that concurrently runs on a device compute units which are called work-items. A memory model needs to be explicitly defined to manage data layout in the memory hierarchy. Details of these models, such as the number of work-items in a work-group or the memory access patterns, have a very strong impact on the program efficiency.

The high level of abstraction of Python allows to conceal from the user the parallel paradigms and the low level implementations of the numerical algorithms. It will participate in our goal to develop a library that is portable to various kind of hybrid architectures of modern computers. Figure 9a shows the global design of the library from a functional point of view using concepts from [START_REF] Labbé | Méthodologie et environnement de développement orientés objets: de l'analyse mathématique à la programmation[END_REF].

In this work, a methodology is proposed for the design of object oriented codes. It relies on a distinction of components in different fields: abstract mathematical concepts are independent of their implementation and of the numerical methods. In order to use the library, the user only needs to describe his problem with high level components -mathematical operators, problem variables and physical domain description. The resolution is performed using default options (Mathematics level). The user chooses the available numerical methods (Numerical methods level) and algorithms (Algorithms level), or develop his own modules, depending on his needs. Thanks to the flexibility and modularity of Python, new modules are readily available in the library (User interface level) through its programming layers, illustrated on Figure 9b. In order to add his modules written in a common programming language (Computations level), the user needs to develop a small Python interface. The module may use a parallel paradigm (HPC level). This design of the library allowed us to easily verify and optimize all the kernels presented below.

Implementation on GPU

The literature contains several GPU implementations of mesh-free particle methods. In the context of Smoothed Particles Hydrodynamics (SPH) for gas dynamics, [START_REF] Valdez-Balderas | Towards accelerating smoothed particle hydrodynamics simulations for free-surface flows on multi-gpu clusters[END_REF] presents an implementation on multi-GPUs cluster for 32 millions of particles distributed on 4 GPUs. A mesh-free solver for incompressible flows using a Fast Multipole Method accelerated on a large number GPU is presented in [START_REF] Yokota | Petascale turbulence simulation using a highly parallel fast multipole method[END_REF]. These authors report reference a 0.5 petaflop/s computation for 2048 3 particles distributed on 2048 cards for homogeneous isotropic turbulence simulations.

The advection equation is a key component in the simulation of gas dynamics. In [START_REF] De | Gpu-accelerated numerical simulations of the knudsen gas on time-dependent domains[END_REF], the authors show that porting the resolution of Knudsen gases by particle methods to the GPGPU results in an important speedup in the simulations. Remeshed vortex methods have also been ported on GPU within the Parallel Particle Mesh library [START_REF] Sbalzarini | PPM-a highly efficient parallel particle-mesh library for the simulation of continuum systems[END_REF]. The implementation presented in [START_REF] Rossinelli | Vortex methods for incompressible flow simulations on the GPU[END_REF] is devoted to two-dimensional incompressible flows in vorticity formulation, using a multigrid finite-difference solver for the Poisson equation and the Λ 2,1 remeshing formula for the advection. This reference reports performances in single precision floating point from 25 frames/iterations per seconds (FPS) with 1024 2 particles -for their fastest solver, Euler time stepping and coarse multigrid -to 3 FPS -for their most accurate, fourth order Runge-Kutta time stepping and accurate multigrid. In [START_REF] Rossinelli | GPU accelerated simulations of bluff body flows using vortex methods[END_REF], penalization was added to deal with complex boundaries and a FFT solver was used instead of the multi-grid solver. The hand coded interpolation between the Cartesian grid and the particles was dropped in favor of OpenGL shaders with improved performances. For 1024 2 particles the full Navier-Stokes solver reached 22 FPS in single precision, which represents a speedup against a similar CPU solver of about 30.

As particle-mesh interpolations are core functions in RPM, some works have been devoted to these operations. In [START_REF] Büyükkeçeci | A portable opencl implementation of generic particle-mesh and mesh-particle interpolation in 2d and 3d[END_REF][START_REF] Rossinelli | Mesh-particle interpolations on graphics processing units and multicorecentral processing units[END_REF] an implementation based on OpenCL Images objects, allowed to obtain, for a 4096 × 2048 problem size, a maximum speedup of 45 against a single-threaded CPU implementation, in double precision, and 150 in single precision. The speedup drops to 2.8 when comparing to a multi-threaded CPU version (9.4 in single precision). In [START_REF] Büyükkeçeci | A portable opencl implementation of generic particle-mesh and mesh-particle interpolation in 2d and 3d[END_REF], an efficient data structure enables interesting performances on GPU for tensorial particle-mesh interpolations in 2D and 3D. This reference reports that linear interpolation exhibits lower speedups than higher order schemes and that 3D interpolations are more efficient than their 2D counterparts.

The method presented in this article differs from the previous works by the fact that we use a dimensional splitting and higher order remeshing formulas. The dimensional splitting reduces 2D and 3D problems to a collection of 1D problems. The complexity will therefore remain linear with the stencil width, but this strategy has some implications that we will see below in the memory management.

In this work, we intend to produce an efficient GPU implementation of the particle advection and remeshing steps that will be used, in future work, as part of a fluid solver. The following results will demonstrate that high order remeshing formulas can be implemented in an efficient way. As in the previous section, particles are advanced with a fourth order Runge-Kutta method. The velocity field is given at grid points at the beginning of each time-step and is linearly interpolated at particle locations corresponding to the intermediate Runge-Kutta sub-steps. We use a second order time-splitting which for each time-step alternates advection successively along the directions x, y, z, y, x. The time-step of the advection is ∆t/2 in the x and y directions and ∆t in the z direction.

From method to algorithm

In order to better understand the library performances and behavior with respect to the number of particles, we analyze the computational complexity of the algorithm. The complexity C is given for one advection and one remeshing step for N d particles, where N is the number of particles in one direction and d is the number of directions. C depends on several parameters: c, the number of scalar quantities or vector components carried by particles; the order of the dimensional splitting (second order in our case); N s , the remeshing stencil width; and p, the polynomial degree of the formula :

C(N s , p, c) = (C A + C R (N s , p, c))(2d -1)N d = O N d [Operations], (5.1)
where C A and C R represent the advection and remeshing complexity for one particle. Advection with a fourth order Runge-Kutta scheme and linear velocity interpolation leads to the following estimate

C A = 4 × 2 + 3 + 3 × 9 + 5 = 43 [Operations].
positions computation linear interpolations weights combination

(5.

2)

The complexity of remeshing algorithms takes into account the polynomial evaluations, performed with the Horner's rule rearrangements. The total amount of operations corresponds to N s polynomials of degree p, c linear combinations of the N s elements and the N s grid point indices.

C R (N s , p, c) = N s (2p + 1) + 2N s c + 3 + N s = 2N s (p + c + 1) + 3 [Operations].
Horner's rule linear combinations grid indices

(5.3)

Similarly, we study the memory access complexity. We obtain the following evaluations

M A = (3 × 2 + 1 + 1 )P = 8P [Bytes],
read velocity write particle position

(5.4) M R (N s , c) = ( 1 + (1 + N s )c )P = ((1 + N s )c + 1)P [Bytes],
read particle position read and write particle quantity

(5.5)

D(c) = 2cP [Bytes],
read and write

(5.6)

M (N s , c) = (D(c)+M A + M R (N s , c))(2d -1)N d = P ((1 + N s )c + 5)(2d -1)N d [Bytes], (5.7)
where P equals to the size of a floating point number in bytes (4 bytes for single precision or 8 bytes for double precision). D(c) denotes the memory complexity of a copy or a transposition. In this expression, reading and writing access are supposed to be equivalent.

These complexities do not take into account the algorithms optimization for targeted platforms. For instance, vectorized access or fma operations are not considered here.

As already mentioned semi-Lagrangian particle methods can be viewed as an alternative to forward semi-Lagrangian methods developped in [START_REF] Crouseilles | A forward semi-lagrangian method for the numerical solution of the vlasov equation[END_REF]. The use of explicit remeshing kernels significantly increase the parallel efficiency of the method. In the method described in [START_REF] Crouseilles | A forward semi-lagrangian method for the numerical solution of the vlasov equation[END_REF], one needs to solve for each 1D problem a periodic tridiagonal system. [START_REF] Zhang | Fast tridiagonal solvers on the gpu[END_REF] compares various implementations for solving tridiagonal systems on GPU. They show that a cyclic reduction algorithm requires at least log 2 N -2 algorithmic steps and 5N accesses to global memory. For the Λ 4,2 formula, with stencils using 6 points, the remeshing algorithm only requires 1 algorithmic step, and 2N accesses to global memory. For problems with at least 64 particles in one direction, the semi-Lagrangian particle methods proposed in this paper require minimal access to global and shared memories.

From algorithm to GPU implementation

The choice of the OpenCL execution model constrained us to an in-order OpenCL kernels call sequence.

The main bottleneck of this model is the host-device data transfer rate. The data should stay on the device as long as possible. However, host-device transfers are mandatory for regular checkpointing and shared computations between CPU and GPUs or multiple GPUs. The index space for calling OpenCL kernels follows the intrinsic dimension decomposition of the method: each work-group is in charge of a 1D problem. The kernels performances are strongly connected to the layout of work-items in each work-groups. Section 5.3 will detail the different strategies for the work-items layout.

The computation of one splitting step requires only the velocity component and the 1D coordinate of the particles in the corresponding direction. The velocity field needs to be stored as structured arrays. The maximum bandwidth for accessing global memory is reached for contiguous data. Arrays are thus stored contiguously in the direction of use. Figure 10 illustrates the storage of the velocity components with respect to the grid layout for a 2D problem. This setup imposes the advected quantities to be transposed when alternating splitting directions.

Algorithm 1 summarizes the GPU implementation, where each part of the inner loop is performed using an OpenCL kernel.

(1,1) (1,2)

• • • (1,N ) (2,1) (2,2) • • • (2,N ) . . . . . . . . . . . . (N ,1) (N ,2) • • • (N ,N ) x y 2D grid Vx (1,1) (2,1) • • • (N ,1) (1,2) (2,2) • • • (N ,2) • • • (1,N ) • • • (N ,N ) Vy (1,1) (1,2) • • • (1,N ) (2,1) (2,2) • • • (2,N ) • • • (N ,1) • • • (N ,

Description of the OpenCL kernels

All our tests were performed on a NVIDIA Tesla K20m. The main specifications of this device are given in Table 2. Obviously, the global and local memory sizes and work-group size set some limits in our performance analysis. Each kernel used in Algorithm 1 must be studied separately to obtain the most efficient execution model characteristics such as e.g. work-group size, work-item workload and memory levels usage. We considered single and double precision implementation of the various remeshing kernels considered in the previous sections. Depending on the type of kernels, different optimizations -concerning OpenCL vector types, built-in functions, temporary private variables, bank conflicts avoidance, .

Memory management related to method and hardware

From the data on the grid, we update the advected quantities on particles. Data initialization depends on the alternating splitting direction in which advection and remeshing are performed. Computation efficiency requires the data to be contiguous in the working direction. Therefore, one needs to change the data layout in memory for each direction. Data initialization consists in transpositions except in the cases of unchanged successive directions where a simple copy is sufficient. For a 2D problem, efficient transposition algorithms on GPU have been designed in [START_REF] Ruetsch | Optimizing matrix transpose in cuda[END_REF]. This algorithm reaches 96 percent performance of the simple copy bandwidth for a 2048 2 square matrix in single precision. We implemented this algorithm and extended it to 3D. In order to achieve good performances, global memory access should be contiguous, fit the bus width and avoid memory camping. Memory camping occurs when successive work-groups operate on the same global memory page. It affects global efficiency of the algorithm and should be avoided as much as possible. Each work-group is handling a subset of the data, called a tile, which is stored in local memory.

The tile is transposed locally, requiring bank conflict free access.

We denote by XY transposition the transposition operating on the first and second directions in which the data are contiguous. Similarly, XZ transposition operates on the first and third directions.

The XY transposition follows [START_REF] Ruetsch | Optimizing matrix transpose in cuda[END_REF] except that, in 3D, it is performed on every slices along Z direction. Figure 11a illustrates this transformation. Each work-group loads a tile (grey square) into local memory (red). Each work-item transposes its own subset of data by writing data to their correct positions (green). Contiguous access in global memory is ensured by use of the tile and the bus width is filled up by wrapping elements into OpenCL vector types. Memory camping is avoided by roaming the global memory arrays diagonally -in dashed lines in the Figure . Finally a one row padding is used in tiles to avoid bank conflicts.

Similarly, XZ transpositions are performed along slices in Y direction. The work of [START_REF] Ruetsch | Optimizing matrix transpose in cuda[END_REF] is extended to deal with 3D data. Using same conventions as previously, Figure 11b sketches this transformation. To avoid memory camping, the global memory is roamed in XZ planes. Because of the larger strides when going along Z direction than in Y, we used a cubic tile to reduce the number of large strides. Performances of the kernels associated to copy and transposition are presented on Figures 12a and12b. These results concern fully optimized kernels and the best configurations. From Figure 12a we conclude that using 256 work-items in work-groups lead to performances close to 72 percent of the theoretical bandwidth. The remaining gap can be explained by the presence of the NVIDIA Error Correction Code (ECC) that checks and eventually corrects bit errors in data storage and transmissions. The performances should increase when disabling ECC because it uses a significant part of the bandwidth and 12.5 percent of the memory1 to work. One can notice that the performance of the XY transposition is better in 3D than in 2D, very close to the performance of the copies for work-group size up to 512. This is due to smaller stride when moving work-group diagonally in 3D. In the optimal configuration, with respect to the work-group size, the XZ transposition exhibits performances similar to the 2D XY transposition.

Advection

This kernel computes the particles positions by means of a Runge-Kutta scheme and a 2D or 3D linear interpolation of velocities on particles from known grid values. As suggested by (5.4), the fourth order Runge-Kutta method requires 8 reads in the velocity field. We bring it to a single read using local

Remeshing

Remeshing particles on the grid involves particle positions and strengths that are read once in global memory. Each work-group computes new quantities on the grid in a local buffer to avoid multiples writes in global memory, as shown on Figure 15. Again, the number of global memory access decreases and (5.5) becomes:

M R,GPU = (2c + 1)P [Bytes] (5.9)

In case of large velocity gradients, it may happen that several particles have the same remeshing points.

To avoid concurrent writings, compute load is organized so that any pair of work-items do not remesh contiguous particles.

As for the advection benchmarks, Figure 14b shows that performances increase with the work-group size until it reaches a plateau.

Monolithic advection and remeshing

RK + interpolation

Velocity component

Particle quantity

Grid quantity Cache

Cache

Λp,r This kernel combines the advection and the remeshing computations in a single kernel. The particle positions are computed directly by the kernel, which reduces global memory occupation. The number of OpenCL kernels launching is also reduced. However, it requires more local memory, as illustrated on Figure 16. For this kernel, the global memory complexity, given by (5.7) is reduced to: .11) for separate kernels.

M GPU,Monolithic (N s , c) = P (2c + 1)(2d -1)N d [Bytes] (5.10) instead of: M GPU (N s , c) = P (2c + 3)(2d -1)N d [Bytes] ( 5
Again similar performances are obtained for monolithic kernel where results are close to the sum of the separate ones. Because of local memory limitations to 48 kBytes, we were not able to compute the 4096 2 double precision case (which would require 64 kBytes).

Performances

The whole code is profiled with the best configurations determined by the previous benchmarks. In this section, for a sake of simplicity we have restricted the discussion to calculations with double precision arithmetic.

Figure 17 gives the global performances for 2D and 3D problems and different remeshing formulas. The measured time includes kernel execution and OpenCL host code such as queue management, kernels launching, profiling events, . . . . Figure 17a shows that the time needed to compute one step for one particle is nearly constant, especially for small problems. This is consistent with (5.1), since the complexity per particle and per iteration, (C A +C R (N s , p, c))(2d-1), does not depend on the problem size. The variations are due to built-in mechanisms of the device depending on the memory occupancy. In 3D, Figure 17a shows that, except for the 64 3 resolution, the computational time is also constant.

Similarly, the global complexity, C(N s , p, c), in (5.1) grows linearly with the problem size, N d . This linearity is observed for the time per iteration on Figures 17b and17d. For the 256 3 and 4096 2 resolutions, performances range from 10 to 20 FPS. We provide on Figure 18 profiling details for the largest cases, corresponding to the 256 3 and 4096 2 particles. As expected, most of the time is devoted to particle remeshing. Remeshing times are grouped with respect to the stencil width. Table 3 exhibit the remeshing computational time per particle averaged -from data on Figures 17a and17c -on different problem sizes together with the arithmetic complexityfrom (5.3) -and the stencil width for the different remeshing formulas. In the Ratio columns, we compute the costs using the Λ 2,1 formula as a reference. We observe that the cost of remeshing depends mostly on the stencil width and only slightly on the polynomial computational complexity. The reason is that memory access are a bottleneck for remeshing kernels and this confirms the efficiency of GPU to accelerate algorithms with high computational intensity. The number of arithmetic operation plays only a marginal role in performances. Note that if the traditional tensor-product formulas were used in 3D, instead of the directional splitting used here, the ratio between the kernels Λ 6,6 and Λ 2,1 would be of order 17b for the 2D cases and Figure 17c for the 3D cases.

The dimensional splitting is also advantageous from the point of view of memory access. The number of operations required for a tensorial formula is of order N 3 s instead of 5N s for the splitting. In a tensorproduct remeshing algorithm, 2N 3 s particles would need to be transferred to and from the global memory. Cache optimization allows to reduce the number of transfer per particle, but due to the limited size of the shared memory, this reduction is limited. For instance, with the GPU described in Table 2, one can only load 12 3 particles in shared memory. This limitation would cause the number of transfers with global memory performed at each particle location to increase. This would significantly impact the overall performance of the remeshing step.

Our benchmarks showed that the computational time remains linear with the stencil width. Increasing the order of the numerical method enables us to leverage the processing power of the GPU. From an applicative point of view, the accuracy of the results is improved either by increasing the problem size or increasing the numerical order of the method while keeping the problem size constant, which allows to relax the constraints due to the global memory size on the GPU. Note that the remeshing formulas involve algebraic fractions for the polynomials coefficients whose numerator and denominator are increasing with the order. In the present implementation, we kept this coefficients under their rational form but this could introduce numerical errors if high order kernels are used.

+

Figure 1 .

 1 Figure 1. Advection equation with data (4.1). Initial condition (solid curve) and solution at t = √ 3 (dashed curve).

Figure 4 .

 4 Figure 4. Interface u = 0 for the advection field (4.2) with N x = N y = 256. Top-left: t = 6, kernel Λ 2,1 ; Top-middle: t = 6, kernel Λ 6,4 ; Top-right: exact front-tracking solution; Bottom-left: t = 12, kernel Λ 2,1 ; Bottom-right: t = 12, kernel Λ 6,4 . On the bottom pictures the red curve represents the exact (initial) interface at t = 12.

  otherwise. with r a = 0.1, r b = 0.25 and C = [2/(r a -r b )] 8 .

Figure 5 .

 5 Figure 5. Same experiment as in Figure 4 with N x = N y = 256, the kernel Λ 6,4 and CF L = 30. Left picture : interface at t = 6, right picture : interface at t = 12 (red curve is the exact solution).

Figure 8 .

 8 Figure 8. Isosurface u = 0.5 at t = 4.0 for the test case of the sphere in the flow (4.4)-(4.6) with N = 256 and CF L = 30. Left picture: kernel Λ 2,1 ; right picture: kernel Λ 8,4 .

Figure 9 .

 9 Figure 9. Application layout.

Figure 11 .

 11 Figure 11. Transposition kernels.

Figure 14 .

 14 Figure 14. Benchmarks for computational kernels in single (SP) and double precision (DP).

Figure 15 .

 15 Figure 15. Remeshing algorithm.

Figure 16 .

 16 Figure 16. Monolithic advection and remeshing algorithm.

Figure 17 .

 17 Figure 17. Overall performances for different problem sizes in double precision. (lower is better). Times are averaged over 20 independent executions.

Table 1 .

 1 Kernels of various regularity, moment properties and complexity. In bold, the kernels that are considered in the numerical experiments of Section 4 and for which analytical formulas are given in the appendix.

	7	[-2; 2]

  ∆x) l+m where d i are coefficients only depending on a, we get

	k

  N )

	Storage
	Figure 10. Velocity field storage.
	Algorithm 1: Particle method
	while t T end do
	Computation of velocity field a G {on grid G} ;
	{Solve 1D problems in each direction};
	foreach splitting direction do
	Initialize particles ;
	Compute particle positions {from corresponding velocity component} ;
	Compute particle remeshing {in corresponding direction to grid G};
	t ← t + dt

Table 2 .

 2 NVIDIA Tesla K20m main features.

	. . -were

  Benchmarks for this kernel are shown on Figure14a. The maximum work-group size allowed on the device is 1024, as given in Table2. Performances increase with the work-group size until they reach a plateau around 128 work-items in 3D and 512 in 2D. Note that constraints on the global memory size and the bus width limits 3D cases to small work-group sizes.

				150									150
			Bandwidth (GBytes/sec)	100								Bandwidth (GBytes/sec)	100	2D SP 2D DP	3D XY SP 3D XY DP
						2D SP	3D SP				3D XZ SP
						2D DP	3D DP				3D XZ DP
				32	64	128	256	512	1024			32	64	128	256	512	1024
						Work-Group size						Work-Group size
					(a) Copy kernels						(b) Transposition kernels
		Figure 12. Performances of Copy, XY and XZ transpositions kernels for 2D and 3D,
		single (SP) and double precision (DP).
	memory as a cache as illustrated on Figure 13. After filling up the cache, each work-group performs the
	computations from this local buffer. Results are directly stored contiguously in global memory. Using
	local memory, (5.4) becomes:
										M A,GPU = 2P [Bytes]	(5.8)
												RK +
											interpolation
					Velocity component		Cache		Particle position
							Figure 13. Advection algorithm.
	Time per particle (ns/part)	0.1 1 10					Time per particle (ns/part)	0.1 1 10				Time per particle (ns/part)	0.1 1 10
		4	8 16 32 64 128 256 512 1024			4	8 16 32 64 128 256 512 1024	4	8 16 32 64 128 256 512 1024
				Work-Group size						Work-Group size	Work-Group size
			4096 2 SP 4096 2 DP 2048 2 DP	256 3 SP 256 3 DP					4096 2 SP 4096 2 DP 2048 2 DP	256 3 SP 256 3 DP	4096 2 SP 2048 2 DP	256 3 SP 256 3 DP
		(a) Advection kernels (RK4)			(b) Remeshing kernels (Λ 2,1 )	(c) Monolithic kernels (RK4 + Λ 2,1 )

Table 3 .

 3 8 instead of 2. Figure 18. Code profiling of one time step for different problem sizes in double precision. R (N s , p, c) Ratio Points Ratio Average Ratio Average Ratio Influence of remeshing formulas on performances. Times are averaged among the different runs reported in Figure

		0.1					copy
							transpose XY
	Computational time (sec)	0.02 0.04 0.06 0.08					transpose XZ advection RK4 remeshing Λ2,1 remeshing Λ2,2 remeshing Λ4,2 remeshing Λ4,4 remeshing Λ6,4 remeshing Λ6,6
		0	2D 4096 2	3D 256 3
		Arithmetics	Stencil		2D Time	3D Time
	Formula C Λ 2,1	43	1	4	1	3.13	1.00	3.88	1.00
	Λ 2,2	59	1.37	4	1	3.17	1.01	3.99	1.03
	Λ 4,2	87	2.02	6	1.5	4.22	1.35	5.81	1.50
	Λ 4,4	135	3.14	6	1.5	4.49	1.43	6.28	1.62
	Λ 6,4	179	4.16	8	2	5.22	1.67	7.88	2.03
	Λ 6,6	243	5.65	8	2	5.62	1.79	8.78	2.26

From http://www.nvidia.com/object/tesla-servers.html

Acknowledgements

The authors are grateful to M. Bergdorf, D. Rossinnelli and P. Koumoutsakos for enlightening discussions on GPU implementation and the derivation of high order kernels. The first author gratefully acknowledges support from Institut Universitaire de France. This research has been partially supported by the Agence Nationale pour la Recherche (ANR) under Contracts No. ANR-2010-COSI-0009 and ANR-2010-JCJC-091601.

Unlike in [START_REF] Büyükkeçeci | A portable opencl implementation of generic particle-mesh and mesh-particle interpolation in 2d and 3d[END_REF], for an equivalent number of particles our 3D and 2D implementations exhibit similar performances, as they both rely on 1D problems. The global efficiency of our implementations could be further improved by gathering 1D problems to better fit the device characteristics. To conclude this section and give an idea of the compared efficiency of GPU and multi-threaded CPU implementations of the particle method, let us mention that the algorithm just described yields speedups ranging between 20 and 25, depending on the remeshing formula, against the optimized Fortran MPI code used in [START_REF] Lagaert | Hybrid spectral particle method for turbulent transport of passive scalar[END_REF] running on 8 Xeon E5-2640 cores.

Conclusion

In this paper we have discussed and analyzed remeshed particle methods, from the point of view of forward semi-Lagrangian methods. We have given a systematic consistency analysis of these methods, based on the regularity and moment properties of the remeshing kernels. These consistency results are supplemented by a stability analysis valid for a class of second order and fourth order methods, under a Lagrangian CFL condition which is independent of the grid size.

Remeshed particle methods are designed to deal with advection dominated problems and we have outlined a general computational framework in which they can be combined with other numerical methods for hybrid calculations. In this context, GPU are ideally suited to the local nature of RPM. We have presented a strategy to optimize the efficiency of RPM on GPU and discussed the complexity of theses algorithms in function of the kernel accuracy. Our results show in particular that, for a given number of conserved moments, increasing the smoothness of the kernels can improve the accuracy of the RPM without added computational complexity. The performance of the algorithms for a linear transport equation in double precision ranges between 10 and 20 FPS for grid resolutions of 4096 2 in 2D or 256 3 in 3D resolutions and kernels with order of accuracy varying between 1 and 6.

Although the focus in this paper was the implementation on GPU of particle methods, the framework we have described is designed with hybrid architectures in mind. Multi-CPUs and multi-GPUs implementations of semi-Lagrangian particle methods are the topic of ongoing work.

Appendix

For a sake of completeness, we give here the analytical formulas for the kernels used in this work.