Michel Mehrenberger

Laura S Mendoza

Charles Prouveur

Eric Sonnendrücker

Solving the Guiding-Center model on a regular hexagonal mesh *

This paper introduces a Semi-Lagrangian solver for the Vlasov-Poisson equations on a uniform hexagonal mesh. The latter is composed of equilateral triangles, thus it doesn't contain any singularities, unlike polar meshes. We focus on the guiding-center model, for which we need to develop a Poisson solver for the hexagonal mesh in addition to the Vlasov solver. For the interpolation step of the Semi-Lagrangian scheme, a comparison is made between the use of Box-splines and of Hermite finite elements. The code will be adapted to more complex models and geometries in the future.

Introduction

In magnetic fusion applications the embedded closed magnetic flux surfaces play an important role and introduce an important anisotropy [START_REF] Angelino | Role of plasma elongation on turbulent transport in magnetically confined plasmas[END_REF]. For this reason one gets favorable numerical properties when grid points align on the concentric magnetic flux surfaces. When trying to do this with a mapped cartesian grid, one ends up with a polar coordinates mesh (when the flux surfaces are circles) or something topologically equivalent. This yields smaller and smaller cells when getting closer to the center as well as a singularity at the center. This is numerically far from optimal.

Different strategies have been implemented to avoid these singularities, we can cite among others: the iso-parametric analysis approach done by J. Abiteboul et al. [START_REF] Abiteboul | Solving the vlasov equation in complex geometries[END_REF] and A. Ratnani [START_REF] Ratnani | Isogeometric analysis in plasmas physics and electromagnetism[END_REF] or N. Besse and E. Sonnendrücker's work with unstructured meshes [START_REF] Besse | Semi-lagrangian schemes for the vlasov equation on an unstructured mesh of phase space[END_REF]. The methods presented in these papers are particularly interesting as not only they avoid singularities but also they are extremely flexible and can be easily adapted to more complex geometries. However, even if these two approaches are different, they suffer from additional cost due to the numerical complexity coming either from the computation of the inverse of a mapping, the advection of the derivatives or the localization of the feet of the characteristics. As for the Poisson equation, a recent study 1 relevant to our problem was made by T. Nguyen et al. [START_REF] Nguyen | A comparative study of several classical, discrete differential and isogeometric methods for solving poisson's equation on the disk[END_REF] to compare different solvers on the disk, where the Iso-geometric approach stood out as a competitive method.

There are three kinds of regular pavings of the plane: using squares, equilateral triangles or hexagons. When considering meshes, the dual mesh of a square mesh, i.e. the mesh generated when taking the Voronoi cells of every point of the original lattice, is a shifted square mesh and the regular triangle mesh is the dual of the regular hexagonal mesh (See Figure 1). Figure 1: A square mesh and its dual (left), and a hexagonal mesh, in solid blue, and its dual mesh, the triangular tessellation, in dashed black (right).

Tiling a regular hexagon into triangles yields a mesh of equilateral triangles having all the same area. Such a mesh was first introduced for numerical simulations in [START_REF] Sadourny | Integration of the nondivergent barotropic vorticity equation with an icosahedral-hexagonal grid for the sphere[END_REF]. An application to particle methods is proposed in [START_REF] Chatelain | Isotropic compact interpolation schemes for particle methods[END_REF]. This grid can be easily mapped to a circle by slightly stretching the edges of the hexagon. Indeed, the lattice is actually composed of concentric hexagons, thus the transformation is a simple scaling of the hexagons points to their circumscribed circle. The scaling coefficient is given by the ratio between the radius of circumscribed circle and the distance from the origin to the point. This yields a nice mesh of a disk with slightly stretched triangles of almost the same size and there is no singularity in any point of the domain. Additionally, such a mesh has a structure with three privileged directions, and uniform steps in each direction, thus it is completely straightforward to localize points within this mesh. The derivatives along the three directions can also be nicely computed using the regular finite difference method along the three directions. And last but not least, there is a spline construction on this mesh, called Box-spline [START_REF] Condat | New optimized spline functions for interpolation on the hexagonal lattice[END_REF]. These splines have a hexagonal support and are invariant by translations along the three directions of the mesh.

The difficulties mentioned for previous approaches, are not present when dealing with this uniform hexagonal mesh. The simulations of our first model, a circular advection, included a comparison with more common methods and a study on the computational efficiency. Moreover a simple and non singular mapping from this mesh can be used to handle more complex settings like the surface aligned meshes needed for tokamak simulations. This point will be left for further studies.

In this work, we focus on adapting the Semi-Lagrangian scheme to this hexagonal mesh. This scheme consists basically of two steps: computing the characteristics' origins and interpolating at these points. For the latter, we compare two different approaches: one using Box-splines and the second approach using Hermite Finite Elements. Both interpolation methods, as well as the mesh, are presented in Section 1. In Section 2, we present a simple finite difference Poisson solver adapted to the hexagonal mesh. We introduce a guiding-center approximation of the 2D Vlasov Poisson system [START_REF] Golse | L'approximation centre-guide pour l'équation de vlasovpoisson 2d[END_REF], and the Semi-Lagrangian scheme to solve it, in Section 3. Finally, in Section 4, we compare the results of the scheme using Box-splines with the ones using Hermite finite elements. Besides the guiding-center model, the circular advection model is used to compare the numerical methods.

1 Interpolation on regular hexagonal mesh

The hexagonal mesh

The hexagonal mesh is obtained by tiling a regular hexagon into equilateral triangles. The mesh obtained can be generated by three vectors. These unit vectors are

r 1 = √ 3/2 1/2 , r 2 = - √ 3/2 1/2 , r 3 = 0 1 . (1)
The 2D lattice sites are obtained by the product Rk where R = (r 1 r 2 ) and k = (k 1 , k 2 ) T ∈ Z. To obtain exactly the mesh as in Figure 2, we need to define a few extra parameters: an origin, denoted by P 0 = (x 0 , y 0 ), a radius L which is the distance between the origin and any external vertex of the hexagon and the number of cells N c on any radius.

P 0 r 3 r 1 r 2 P 5 P 6 P 1 P 2 P 3 P 4 Figure 2:
The hexagonal lattice and the vectors r 1 , r 2 , r 3 that generate such a mesh

The mesh is based on uniform hexagons of the first type (see [START_REF] Ulichney | Digital halftoning[END_REF]). For local and global indexing we will use the following convention: the point at the center will be the point of index 0. Following the direction r 1 the next point will be indexed 1, and the indexing will follow in a counter-clockwise motion. And so on, until all the points of the domain have been indexed. See Figure 2. We will denote P i the point of global index i, of cartesian coordinates x i = (x i , y i ) and hexagonal coordinates

k i = R -1 x i .
Besides the fact that the hexagonal mesh contains no singularities, its regularity allows us to localize the characteristics' origins for the Semi-Lagrangian scheme by taking three integer values, similarly to what is done on cartesian grids for which only two integer values are needed. Nevertheless, the accuracy of the method depends heavily on the interpolation method chosen. For example, for a Cartesian grid, it is common to use cubic splines which have shown to give accurate results in an efficient manner [START_REF] Sonnendrücker | The semi-lagrangian method for the numerical resolution of the vlasov equation[END_REF]. In our problem, with the hexagonal lattice, Bsplines do not exploit the isotropy of the mesh (for more information see [START_REF] Mersereau | The processing of hexagonally sampled two-dimensional signals[END_REF]) and are defined by a convolution in 2D, which can't be done for our mesh. Therefore, we need to use another approach. In the following two subsections we present two different strategies: the first one using Box-splines and a second one using Hermite Finite Elements.

Box-splines quasi-interpolation

There are mainly two families of splines that take advantage of the geometry's properties: hexsplines, first introduced in [START_REF] Van De Ville | Hexsplines: a novel spline family for hexagonal lattices[END_REF], and the three directional Box-splines. For a detailed comparison between these two types of splines we will refer to [START_REF] Condat | Quasi-interpolating spline models for hexagonally-sampled data[END_REF]. Based on the latter, we chose to use Box-splines, as the results are more stable. And lastly, also based on the previously cited paper, we decided to use a quasi-interpolation method [START_REF] Lyche | Quasi-interpolation projectors for box splines[END_REF].

Box-splines: General Definition

Box-splines are a generalization of the well known B-splines. They are also piecewise polynomial and they share some properties, such as: compact support, positiveness, symmetry and partition of unity. But, unlike B-splines, Box-splines are defined from a generator matrix Ξ. Therefore, to construct them on the hexagonal lattice, we will use the generator vectors r 1 , r 2 , r 3 . The general definition is [START_REF] De Boor | Box Splines[END_REF][START_REF] Condat | Three-directional box-splines: characterization and efficient evaluation[END_REF]:

Definition 1.1 (Box-splines). Let Ξ be a d × m matrix with non-null columns in R d . A Box- spline χ Ξ associated to the matrix Ξ, is a multivariate function χ Ξ : R d -→ R.
If Ξ is a square invertible matrix, i.e. when m = d and det(Ξ) = 0, we define a Box-spline with the formula below

χ Ξ (x) =    1 |det(Ξ)| if Ξ -1 x ∈ [0, 1) 2 , 0 otherwise.
(

) 2
If Ξ ∪ v is a d × (m + 1) matrix, composed by the m column vectors from Ξ to which we append the vector v, we define the Box-spline χ Ξ∪v by recursion Remark 1.2. The Box-splines can have different degrees in each direction. Thus, there are different definitions of the degree. We consider the definition below, which is specific to the kind of Box-splines we use in this paper.

χ Ξ∪v (x) = 1 0 χ Ξ (x -tv) dt. ( 3
Definition 1.2 (Three-directional Box-spline of degree N ). Let Ξ be a 2 × 3 matrix with nonnull columns in R 2 such that they form a generating set of R 2 . Then, the three-directional Box-spline of degree N of generating matrix Ξ, χ N Ξ , is the Box-spline associated to Ξ, where all three generating vectors have multiplicity N .

In Figure 3, the difference between the first two Box-splines, χ [r1,r2] and χ [r1,r2,r3] lies on the generating matrix. We notice that only χ [r1,r2,r3] and χ 2 [r1,r2,r3] are adapted to the hexagonal mesh and have an hexagonal support. In fact, these are the three-directional Box-splines of degree 1 and 2 in the hexagonal-mesh that will be mentioned further in the following sections.

The quasi-interpolation scheme

The problem is the following: we are given an initial sample on the grid s[k] = f (Rk), where the points Rk belong to our hexagonal mesh, and we need to know the values f (x) where x / ∈ Rk. To this aim we define a spline surface f (x) = c[k]χ N (x -Rk), where χ N are the Box-splines of degree N of matrix Ξ = [r 1 , r 2 , r 3 ] and c[k] are the coefficients associated to them. These coefficients are defined such that f (x) approximates f (x) to a certain order M = 2N or, in other words, the approximation is exact only if f (x) is a polynomial of degree M -1 or less [START_REF] Condat | Quasi-interpolating spline models for hexagonally-sampled data[END_REF]. This is different from the classical interpolation method, where the reconstruction is exact on grid points for all smooth functions. The c[k] coefficients are the Box-splines coefficients, to compute them we are no longer able to solve a matrix-vector system because of the extra degree of freedom given by the quasi-interpolation method. Thus, the c[k] coefficients are obtained, for a grid point x j of hexagonal coordinates k j (i.e. such that x j = Rk j ) by discrete filtering [START_REF] Condat | Efficient reconstruction of hexagonally sampled data using three-directional box-splines[END_REF]

c[k j ] = K-1 =0 s[k V j, ]p[ ] = K-1 =0 s[k + k j ]p[ ], ( 4
)
where s is the initial sample data and p[ ] are K pre-filters which will be defined later on. V j is the set of global indices of the K neighbouring points of P j , thus V j, is the global index of the -th point in our stencil, and k V j, its hexagonal coordinates. We notice that (k + k j ) also designates the same hexagonal coordinates. This is due to the properties of the hexagonal mesh, and the notations we chose in Section 1.1.

Box-splines coefficients

We recall we have formula [START_REF] Bernadou | Méthodes d'éléments finis pour les problèmes de coques minces[END_REF]. Based on the literature available, notably [START_REF] Condat | Quasi-interpolating spline models for hexagonally-sampled data[END_REF], we have chosen for Box-splines of degree 1 (i.e. N = 1) the quasi-interpolation pre-filters p IIR2 which seem to give better results within a competitive time. The pre-filter p IIR2 [ ] of the point of local index , for splines of degree 1, is defined by

p IIR2 [ ] =            1775
Here K = 19, and for Box-splines of degree 2, K = 61. For higher degrees, we refer to the previously mentioned papers (particularly [START_REF] Condat | Efficient reconstruction of hexagonally sampled data using three-directional box-splines[END_REF]). To get the Box-splines coefficients, we use (4), where p[ ] = p IIR2 [ ].

Optimizing the evaluation

Finally, we have all the elements for the approximation of a function f with Box-splines of degree N

f (x) = k∈Z 2 c[k]χ N (x -Rk). (6)
Even if we limit our sum to the vector k that defines our domain, we would like to take advantage of the fact that the splines χ N are only non-zeros in a limited number of points. Therefore we need to know the indices k such that χ N (x -Rk) = 0. For this purpose we will use the strategy suggested in [START_REF] Condat | Quasi-interpolating spline models for hexagonally-sampled data[END_REF]: to start we need to obtain the indices on the coordinate system generated by R: k = [ u v ] where (u v) T = R -1 x. Thus, for example, in the case N = 1, we only need 4 terms associated to the encapsulating rhomboid's vertices: R k, R k + r 1 , R k + r 2 and R k + r 1 + r 2 . Finally we obtain:

f (x) = c[ k] χ 1 (x -R k) + c[ k + [1, 0]] χ 1 (x -R k -r 1 ) + c[ k + [0, 1]] χ 1 (x -R k -r 2 ) + c[ k + [1, 1]] χ 1 (x -R k -r 1 -r 2 ). ( 7
)
Remark 1.3. The χ 1 spline has a support of radius the unity, thus one of the elements of ( 7) is null. But this formula allows us to keep a short general formula for all points on the mesh without having to compute the indices of the cell to which x belongs to.

Remark 1.4. For the Box-splines of degree 2, in equation [START_REF] Birdsall | Plasma physics via computer simulation[END_REF] there would be 16 coefficients to compute (see Figure 4) from which 4 would be null terms.

Algorithms

Finally, we give the complete algorithms for quasi-interpolating a function f with Box-splines, which relies on two further algorithms: an algorithm for computing the quasi-interpolating coefficients, and an algorithm for computing the Box-splines.

→ k = k + k ; if k ∈ Ω then Compute global index of k -→ index; Get cartesian coordinates of point at index -→ x = R k; result = result + coeffs(index) * BoxSplineValue(x,
∈ Ω do Initialize coeffs(i ) = 0; for = 0 to K -1 do Compute global index of point at (k i + k ) -→ index; coeffs(i ) = coeffs(i ) + sample(index) * PreFilter( ) end end
As for the third and last algorithm, the procedure to compute a Box-spline of a certain degree on a point, we refer to [START_REF] Condat | Three-directional box-splines: characterization and efficient evaluation[END_REF]. We used the same evaluation for Box-splines. Furthermore, we focused on degree 2 splines as they have been optimized in this same article.

Hermite Finite Elements interpolation

Another type of interpolation method is the Hermite Finite Element interpolation. In which, to interpolate a function f at a point x of barycentric coordinates (λ 1 , λ 2 , λ 3 ) in the triangle T of vertices S 1 , S 2 and S 3 , we need a finite element with a local interpolation operator Π T . This operator can be defined by the duality product of a set of degrees of freedom Σ T and a set of basis functions B which depend on the barycentric coordinates. For the following section, we define the indices i, j, and k, with the following relations:

i ∈ 1; 3 , j = i[3] + 1, k = j[3] + 1.
Here i[3] (respectively j [START_REF] Banks | A new class of nonlinear finite-volume methods for vlasov simulation[END_REF]) is the rest of the euclidean division of i (resp. j) by 3. Using this notation, for any vertex S i of a cell T , the remaining two vertices are S j and S k .

Several elements have been tested here: The Z9 and Z10 Zienkiewicz elements, the Hsieh-Clough-Tocher reduced (HCT-r) and complete (HCT-c), the Ganev-Dimitrov element, and lastly the Mitchell element. These elements can be found in [START_REF] Guscaglia | Finite element methods for the stokes system based on a zienkiewicz type n-simplex[END_REF], [START_REF] Bernadou | Méthodes d'éléments finis pour les problèmes de coques minces[END_REF], and [START_REF] Mitchell | An introductions to the mathematics of the finite element method[END_REF]. We show in the following sections, specifically how the hexagonal structure simplifies the interpolation with these elements. However, before presenting in details each different type of element, we start this section with a brief overview on how we compute the derivatives, which are needed in all methods.

Derivatives reconstruction

In our setting, only values at the vertices are known, whereas other degrees of freedom are reconstructed through finite difference formula, as in [START_REF] Hamiaz | The semi-lagrangian method on curvilinear grids[END_REF] for example. We focus here in the reconstruction of derivatives in a given direction.

Given a function f ∈ R 2 and a direction h ∈ R 2 , the right derivative (symbolized by the sign +) along the direction h of f at a point x i ∈ R 2 is denoted by f h + (x i ) and is approximated to the order p ∈ N * by the formula

hf h + (x i ) s + =r + b + f (x i + h) (8)
where h = |h|, is the euclidean length of h, r + = -p 2 , s + = p+1

2

, and

               b + = s + κ=r + , κ ∈{0, } (-κ) s + κ=r + , κ = ( -κ) , for = r + , . . . , s + , = 0, b + 0 = - s + κ=r + , κ =0 b + κ . (9)
Formula ( 9) yields, for the orders p = 1 . . . 6, the following coefficients.

p = 1 b + = (-1, 1), p = 2 b + = (-1/2, 0, 1/2), p = 3 b + = (-1/3, -1/2, 1, -1/6), p = 4 b + = (1/12, -2/3, 0, 2/3, -1/12), p = 5 b + = (1/20, -1/2, -1/3, 1, -1/4, 1/30), p = 6 b + = (-1/60, 3/20, -3/4, 0, 3/4, -3/20, 1/60). (10)
Using these coefficients in formula ( 8), we can get the derivatives along each direction r 1 , r 2 , r 3 , by taking h ∈ {±r i , i = 1, 2, 3}. We notice that this will be equivalent to obtain the derivatives along the edges of any cell of our mesh.

The Z9 and Z10 Zienkiewicz elements

The Z9 approach [START_REF] Guscaglia | Finite element methods for the stokes system based on a zienkiewicz type n-simplex[END_REF] uses 9 degrees of freedom which are the values of f at the vertices of the triangle (which are given) and the values of the derivatives in the direction of the edges at every vertex (computed with the procedure explained in the previous section 1.3, Formula ( 8), with

∂---→ SiSj f (S i ) = f h + (S i ), h ∈ {±r i , i = 1, 2, 3}) and p = 6, Σ T = {∀i ∈ 1; 3 , f (S i ), ∂---→ SiSj f (S i ), ∂---→ SiS k f (S i )}. ( 11
)
Z10 uses one more degree of freedom which is the value at C, the center of the triangle:

Σ T = {∀i ∈ 1; 3 , f (S i ), ∂---→ SiSj f (S i ), ∂---→ SiS k f (S i ), f (C)}. ( 12
)
The advantage of using the Z10 element is the gain of one order of precision: it reproduces polynomials of total degree ≤ 3 whereas with Z9 only polynomials of degree ≤ 2 are reproduced. Let us note that although adding one degree of freedom seems harmless, adding the cells' centers represents, for a hexagonal mesh, a computational cost three times higher. In fact, there are twice as many centers as vertices in a hexagonal mesh. Therefore the number of computational points is tripled.

Let us define the basis functions needed to interpolate with the element Z9. First denote by

φ = λ 1 λ 2 λ 3 . (13)
We introduce now ξ i and ξ ij for performances purposes:

ξ i = λ 3 i -φ and ξ ij = λ 2 i λ j + φ 2 .
And thus we find that the basis functions associated to the vertex and derivative degrees of freedom write respectively

φ i = 3λ 2 i -2ξ i and φ ij = h ij ξ ij = hξ ij ,
where h ij is the length of [S i S j ], and φ i is the basis function associated with the value of the function at S i while φ ij are associated with the derivatives in the direction of the edges. The fact that T is equilateral is exploited here by replacing h ij with h since the length of [S i S j ] is constant. Finally for Z9 we have:

Π T (f ) = 3 i=1 [f (S i ) φ i + j =i ∂---→ SiSj f (S i ) φ ij ]. (14)
In the same manner, let us define the basis functions needed to interpolate with Z10: Firstly, we introduce φ i the basis function at the vertices,

φ i = 3λ 2 i -2ξ i -9φ.
Secondly, φ ij , the basis functions derivatives along the edges,

φ ij = h ij (ξ ij - 3 2 φ) = h(ξ ij - 3 2 φ).
And lastly, φ 123 , at the cell's center,

φ 123 = 27φ.
Thus, for Z10 we have:

Π T (f ) = 3 i=1 [f (S i ) φ i + j =i ∂---→ SiSj f (S i ) φ ij ] + f (C) φ 123 . (15)

The HCT elements

The HCT elements were tested as well because of their original feature which is to use a division of the triangle into three sub-triangles. This characteristic is the only difference between the interpolation with the HCT-r and the Z9 element as they both use the same 9 degrees of freedom, and reproduce polynomials of degree ≤ 2. Unsurprisingly, they give quasi-identical results which is why we will not detail the HCT-r interpolation and focus solely on the HCT-c elements. These elements use the same degrees of freedom as HCT-r, the values at the vertices and their derivatives along the edges, but it has additionally the values of the derivatives in the normal direction of the edges at the middle of the respective edge. This adds up to twelve degrees of freedom. The HCT-c elements reproduce polynomials of degree ≤ 3. For a detailed description of these elements on an unstructured mesh we refer to [START_REF] Besse | Etude mathématique et numérique de l'équation de Vlasov non linéaire sur des maillages non structurés de l'espace des phases[END_REF], here we focus on applying these schemes to our mesh whose regularity simplifies greatly the computation of the derivatives (see 1.3.1).

Let us now define its interpolation operator. Let S i be a vertex of the triangle T, then we define respectively l i and m i as the length and the middle of the edge opposite to S i . The exterior normal to the edge opposite of S i is noted ν i . Let G be the barycenter of T, then K l is the sub-triangle made with G, S j and S k , we get

Π K l (f ) = (l+1)mod(3)+1 i=l (l+1)mod(3)+1 j=l j =i [f (S i ) φ l,i + ∂---→ SiSj f (S i ) φ l,ij -n i ∂ νi f (m i )φ ⊥ l,i ], (16)
where n i = |m i -S i |, and φ l,i is the basis function in the sub-triangle K l at the vertex S i , φ l,ij the first derivative of φ l,i along the edge S i S j and finally, φ ⊥ l,i is the derivative along the normal passing through S i to the opposite edge. For computing the derivatives along the edges, we use formula [START_REF] Carlson | Third order accuracy of the 4-point hexagonal net grid[END_REF], with h∂---→ SiSj f (S i ) = hf h + (S i ) and h ∈ {±r i , i = 1, 2, 3}. For computing the normal derivatives at m i , please see Remark 1.5. The basis functions are defined by

B l = Σ l Λ l , (17)
with Σ l the matrix containing the basis functions' coefficients in the sub-triangle K l . Furthermore,

B l = (φ l,i , φ l,j , φ l,k , φ l,ik , φ l,ij , φ l,ji , φ l,jk , φ l,kj , φ l,ki , φ ⊥ l,i , φ ⊥ l,j , φ ⊥ l,k ) T , and Λ l = (λ 3 i , λ 3 j , λ 3 k , λ 2 i λ k , λ 2 i λ j , λ 2 j λ i , λ 2 j λ k , λ 2 k λ j , λ 2 k λ i , λ i λ j λ k ) T .
The matrices Σ l are defined with the eccentricity e i of each edge opposite to the vertex S i of the triangle T:

e i = l 2 k -l 2 j l 2 i
where j and k are the indices of the remaining two vertices of T -S j and S k .

For an equilateral triangle, the eccentricity is null which simplifies greatly Σ l :

Σ l =                     0 0 0 9 2 9 2
0 0 0 0 0

1 2 1 0 -3 2 0 3 3 0 0 3 1 2 0 1 0 -3 2 0 0 3 3 3 -1 12 0 0 5 4 1 2
0 0 0 0 0

1 12 0 0 1 2 5 4 0 0 0 0 0 -7 12 0 0 -1 2 -1 4
1 0 0 0 1 2 3 0 0 -1 4 1 4 0 1 0 0 1 2 2 3 0 0 1 4 -1 4 0 0 1 0 1 2 -7 12 0 0 1 4 -1 4 0 0 1 0 1 2 2 3 0 0 1 4 -1 4 0 0 1 0 1 2 -7 12 0 0 1 4 -1 4 0 0 1 0 1 2 1 4 0 0 -1 4 -1 2 0 0 0 1 1                     . ( 18
)
The advantage of HCT-c and HCT-r is that they do not require more points than the lattices of the hexagonal mesh. This is not the case for the Z10 approach.

Remark 1.5. The normal derivatives at the middle of the edges are computed using the derivatives following the same normal vector but at the vertices of the triangle, which are themselves computed using the derivatives at the vertices but with respect to the edges (i.e. ∂---→ SiSj f (S i )). For example, let ∂ ν f (m) be the normal derivative at the middle of the edge [S i S j ]. We first compute all ∂ ν f (S i + κ --→ S i S j ), where κ = -2, . . . , 3 using a linear combination of ∂ h f (S i ) with h ∈ {r 1 , r 2 }. To obtain the normal derivatives at the middle of the edges, we simply use a Lagrange interpolation of degree 5, such as:

∂ ν f (m) = 3 κ=-2 a κ ∂ ν f (S i + κ --→ S i S j ), where a -2 = a 3 = 3 256 , a -1 = a 2 = -25 256
, and a 0 = a 1 = 75 128 .

The Ganev-Dimitrov element

The Ganev-Dimitrov element reproduces polynomials of total degree ≤ 4 and uses 15 degrees of freedom which are the values of the function at the vertices and at the middle of the edges, plus the value of the derivatives at the vertices in the direction of the other two vertices. The computational cost for this element is four times higher than the HCT-r interpolation because of the computations needed at the middle of the edges: there are on average 3 times more edges than vertices. As a matter of fact, the vertices and the middle of the edges form another hexagonal mesh twice as fine as the original mesh. The reason why we tested such a computationally expensive element is to observe whether or not the gain in precision is interesting compared to the extra computing time allocated. The local interpolation operator is:

Π T (f ) = 3 i=1 [f (S i ) φ i + ∂---→ SiSj f (S i ) φ ij + f (m i ) φ ⊥,0 i -n i ∂ νi f (m i )φ ⊥,1 i ]. ( 19
)
Contrary to ( 16), here we use the values of normal derivatives of the basis function at the middle point of the edge opposite to S i , φ ⊥,1

i , but also the value at these middle points φ ⊥,1 i . The basis function are defined by:

B = ΣΛ, (20)
with:

B = (φ 1 , φ 2 , φ 3 , φ 1,3 , φ 1,2 , φ 2,1 , φ 2,3 , φ 3,2 , φ 3,1 , φ ⊥,0 1 , φ ⊥,0 2 , φ ⊥,0 3 , φ ⊥,1 1 , φ ⊥,1 2 , φ ⊥,1 3 ) T , Λ = (λ 4 1 , λ 4 2 , λ 4 3 , λ 3 1 λ 3 , λ 3 1 λ 2 , λ 3 2 λ 1 , λ 3 2 λ 3 , λ 3 3 λ 2 , λ 3 3 λ 1 , λ 2 2 λ 2 3 , λ 2 3 λ 2 1 , λ 2 1 λ 2 2 , λ 2 1 λ 2 λ 3 , λ 1 λ 2 2 λ 3 , λ 1 λ 2 λ 2 3 ), and, Σ =                          
1 0 0 4 4 0 0 0 0 0 -5 -5 -4 0 0 0 1 0 0 0 4 4 0 0 -5 0 -5 0 -4 0 0 0 1 0 0 0 0 4 4 -5 -5 0 0 0 -4 0 0 0 1 0 0 0 0 0 0 1 0 -0.5 -0.5 0.5 0 0 0 0 1 0 0 0 0 0 0 -1 -0.5 -0.5 0.5 0 0 0 0 0 1 0 0 0 0 0 -1 -0.5 -0.5 0.5 0 0 0 0 0 0 1 0 0 -1 0 0 -0.5 -0.5 0.5 0 0 0 0 0 0 0 1 0 -1 0 0 -0.5 0.5 -0.5 0 0 0 0 0 0 0 0 1 0 -1 0 0.5 -0.5 -0.5 0 0 0 0 0 0 0 0 0 16 0 0 -16 16 16 0 0 0 0 0 0 0 0 0 0 16 0 16 -16 16 0 0 0 0 0 0 0 0 0 0 0 16 16 16 -16 0 0 0 0 0 0 0 0 0 0 0 0 -4 4 4 0 0 0 0 0 0 0 0 0 0 0 0 4 -4 4 0 0 0 0 0 0 0 0 0 0 0 0 4 4 -4

                          . ( 21
)
Remark 1.6. Here the normal derivatives are computed in the same matter as for the HCT elements. Nevertheless, this could probably be improved for the Ganev-Dimitrov elements as we have an additional information: the value of the function at the middle of the edges.

The Mitchell finite element [29] MI(p)

The Mitchell elements use 12 degrees of freedom which are the values of the function at the vertices, the value of the derivatives at the vertices in the direction of the other two vertices, and additionally, it uses mixed derivatives. We will first define the complete Mitchell which uses 15 degrees of freedom, by adding the value of the fourth derivative at the middle of the edges, and then we will set these degree of freedom to zero, as we do not want to compute numerically fourth order derivatives. Given a function f : R 2 → R and a triangle S 1 , S 2 , S 3 , with vertex S i = (x i , y i ) ∈ R 2 , we will define the degrees of freedom through the following auxiliary function g(λ 1 , λ 2 , λ 3 ) = f (λ 1 S 1 + λ 2 S 2 + λ 3 S 3 ), which takes as input the barycentric coordinates.

The first 3 degrees of freedom are the values at the vertices of the triangle:

g 0 = g(1, 0, 0), g 1 = g(0, 1, 0), g 2 = g(0, 0, 1).
Note that g 1 and g 2 are obtained cyclically from g 0 . The following degrees of freedom are the derivatives at the vertices

g 3 = h∂---→ S1S2 f (S 1 ) = lim t→0 + g(1 -t, t, 0) -g(1, 0, 0) t , g 4 = h∂---→ S2S1 f (S 2 ) = lim t→0 + g(t, 1 -t, 0) -g(0, 1, 0) t and g 5 = h (4) (1/2), h(t) = g(t, 1 -t, 0) = 0.
We get the next 6 degrees of freedom cyclically. Mixed derivatives are obtained in the following way

g 12 = h 2 ∂---→ S1S3 ∂---→ S1S2 f = lim s,t→0 + g(1 -t -s, t, s) -g(1 -t, t, 0) -g(1 -s, 0, s) + g(1, 0, 0) st ,
and cyclically g 13 and g 14 . Numerically, they are obtained by applying the derivation procedure (8) to previously computed derivatives (i.e. ∂---→ SiSj f (S i )). Then given g i , i = 0, . . . , 14, there exists a unique polynomial of degree ≤ 4 that has these degrees of freedom, and this defines the complete Mitchell finite element. The basis functions can be explicitly computed as follows. We define

b 0 (λ 1 , λ 2 , λ 3 ) = λ 2 1 (3 -2λ 1 + 6λ 2 λ 3 ), b 1 (λ 1 , λ 2 , λ 3 ) = λ 2 1 λ 2 (1 + 2λ 3 ), b 2 (λ 1 , λ 2 , λ 3 ) = b 1 (λ 2 , λ 1 , λ 3 ), b 3 (λ 1 , λ 2 , λ 3 ) = λ 2 1 λ 2 λ 3 , b 4 (λ 1 , λ 2 , λ 3 ) = λ 2 1 λ 2 2 / 24, σ(λ 1 , λ 2 , λ 3 ) = (λ 2 , λ 3 , λ 1 ) and Π T (f ) = 2 j=0 (g j b 0 + g 3+3j b 1 + g 4+3j b 2 + g 5+3j b 4 + g 12+j b 3 ) • σ j .
This reconstruction is exact for all polynomials of degree ≤ 4. We restrict then this space, setting g 5 = g 8 = g 11 = 0. Thus Σ T = {g i , i = 0 . . . 14} \ {g 5 , g 8 , g 11 }. This corresponds to the space of polynomials of degree ≤ 4 which are of degree ≤ 3 on the edges. As in the previous sections, the derivatives are computed numerically from the values at the points, following the derivation procedure [START_REF] Carlson | Third order accuracy of the 4-point hexagonal net grid[END_REF]. Mixed derivatives are obtained by applying the derivation procedure (8) to previously computed derivatives (i.e. ∂---→ SiSj f (S i )). Note that this interpolation procedure is here the analog on triangles of Hermite interpolation on squares as in [START_REF] Hamiaz | The semi-lagrangian method on curvilinear grids[END_REF]. It shares also the same finite element space as the Box-splines of degree 2. All this motivates the introduction of this element (which came to our knowledge during the writing of this paper).

In the sequel, we will refer to MI(p), for the Mitchell finite elements using a finite difference reconstruction of order p ∈ N * for the computation of the derivatives and mixed derivatives, using the procedure defined in section 1.3.1 (e.g. MI4, MI6, MI17). We introduce this notation, as unlike for the preceding elements, we will make a study with respect to p (which was set to 6 for the other elements).

The Poisson finite-difference solver

When computing the origins of the characteristics with the semi-Lagrangian method for the Vlasov-Poisson or guiding-center models we need to compute the solution of the Poisson equation -∆φ = ρ, φ being the potential and ρ the density. We impose here null Dirichlet boundary conditions. In order to solve this equation, we use a simple finite difference scheme. Since the mesh here is hexagonal, a seven point stencil is used as shown in Figure 2. It is composed of the six vertices of a hexagon plus its center. To compute φ 0 , the value of φ at the center 0, the remaining vertices of the hexagon are used. This particular stencil has the property to give a fourth order scheme at little cost [START_REF] Carlson | Third order accuracy of the 4-point hexagonal net grid[END_REF]. Here is the previously described scheme:

-(φ 1 + φ 2 + φ 3 + φ 4 + φ 5 + φ 6 -6φ 0 ) = 3h 2 4 ρ 0 + h 2 24 (ρ 1 + ρ 2 + ρ 3 + ρ 4 + ρ 5 + ρ 6 ).
Compared to the second order scheme on the same stencil, we notice the only difference to be the second term of the equality:

-(φ 1 + φ 2 + φ 3 + φ 4 + φ 5 + φ 6 -6φ 0 ) = h 2 ρ 0 .
Considering the gain of two order of precision at such little cost, we have used this fourth order scheme to compute φ. Remark 2.1. One difficulty that arises here is to define an indexing that allows the resolution of a "computational-friendly" linear system, i.e. a sparse matrix with the non-null terms close to the diagonal to minimize filling in a Cholesky decomposition. This is done by assigning a number following one hexagonal direction, row after row, similarly to how one proceeds on a Cartesian mesh. Here, however, the difference is that the rows are of variable width resulting in a banded matrix. Therefore the matrix here is not constituted of 7 diagonals which makes the Poisson computation longer than on a Cartesian mesh. The width of the band is directly proportional to the number of cells in the hexagonal domain.

The Backward Semi-Lagrangian Scheme

When solving a Vlasov equation, one usually thinks of Lagrangian methods such as PIC [START_REF] Birdsall | Plasma physics via computer simulation[END_REF]. However these schemes are prone to numerical noise and converge slowly in 1/ √ N as the number N of particles increases, typical of a Monte Carlo integration. Another option to solve the Vlasov equation, are Eulerian methods like Finite Difference, Finite Element or Finite Volume methods [START_REF] Filbet | Comparison of eulerian vlasov solvers[END_REF][START_REF] Zaki | A finite element code for the simulation of one-dimensional vlasov plasmas. i. theory[END_REF][START_REF] Banks | A new class of nonlinear finite-volume methods for vlasov simulation[END_REF]. The downside of this type of method is that there is a numerical limit on the time step, the CFL condition.

With the intent of overcoming the pitfalls of these methods, the Semi-Lagrangian method was introduced, first in numerical weather prediction (see [START_REF] Kalnay | Atmospheric Modeling, Data Assimilation and Predictability[END_REF] and articles cited within it), and then for plasma simulations [START_REF] Sonnendrücker | The semi-lagrangian method for the numerical resolution of the vlasov equation[END_REF][START_REF] Cheng | The integration of the vlasov equation in configuration space[END_REF] and is used also for gyrokinetic simulations of plasma turbulence [START_REF] Grandgirard | A drift-kinetic semi-lagrangian 4d code for ion turbulence simulation[END_REF][START_REF] Kwon | Development of semi-lagrangian gyrokinetic code for full-f turbulence simulation in general tokamak geometry[END_REF]. This scheme consists in fixing a Eulerian grid in phase-space and following the trajectory of the equation's characteristics in time to compute the evolution of the distribution function. The advantages of this scheme are the possibility of taking large time steps and its stability. However it is still quite costly in high dimensions (5 or 6D phase space) where the PIC method still largely dominates. Lastly, we can point out that there are many types of Semi-Lagrangian solvers (e.g. depending on the trajectories: Backward or Forward; depending on degrees of freedom on which it is based: grid points, cell average, ...). We have chosen here to use the classical Backward Semi-Lagrangian (BSL) method.

Our model

We consider here a 2D linear or non linear advection equation, with a divergence free advection field A, which can be written in general form

∂ρ ∂t + A • ∇ x ρ(x, t) = 0, ( 22
)
where A is divergence free (i.e. ∇ • A = 0) and the density ρ is known at the initial time (i.e. ρ(x, 0) = ρ 0 (x) is known). The advection field A will either be given and known for all times or it will depend on an electric potential computed from the solution of a Poisson equation, i.e.

A = -∂φ ∂y ∂φ ∂x

, with -∆φ = ρ.

When the advection coefficient is obtained from a Poisson equation, the model is known as the guiding-center model. We will apply the backward Semi-Lagrangian scheme to solve both the advection in a given field and the guiding-center model.

Computing the origin of the characteristics

We consider the model ( 22) on a 2D hexagonal domain, discretized with the hexagonal mesh. The points of the lattice are denoted x = (x, y). The distribution function ρ(x, t) is known on all grid points at the initial time t = 0. Let A x and A y be respectively the first and second components of A. We proceed to apply the BSL method to the Vlasov equation ( 22): First, we need to compute the origin of the characteristics ending at the grid points. These are defined for a given time s ∈ R by

   dX dt = A X(s) = x ⇐⇒            dX 1 dt = A x dX 2 dt = A y X 1 (s) = x, X 2 (s) = y. ( 23
)
The solutions (X 1 , X 2 ) of ( 23) are called the characteristics associated to the Vlasov equation. Now denoting by t n = n∆t, for a given time step ∆t, and X n = X(t n ) for any n, and setting s = t n+1 . The origin, at time t n , X n of the characteristics ending at the grid point X n+1 = x can then be computed by any ODE solver, typically a Runge-Kutta solver if A is known for all times. In the case of the guiding-center model, we use a second order scheme, which is the implicit Adams-Moulton scheme of order two [START_REF] Filbet | High order time discretization for backward semi-lagrangian methods[END_REF], to compute the origin of the characteristics,

X n+1 -X n ∆t = 1 2 (A n+1 + A n ).
Where

A n = A(t n , X n ).
The difficulty here is that A(t n+1 , X n+1 ), depends on ρ n+1 and is unknown, thus an approximation *

A of A at time t n+1 is made thanks to previous computations:

* A = 2 A(t n , X n+1 ) -A(t n-1 , X n+1 ).
The unknown X n is found by solving:

     X n+1 -X n ∆t = 1 2 ( * A + A n ), X(s) = x.
Remark 3.1. Since we need A(t n-1 ), the first step is done using the implicit Euler time scheme.

Updating the distribution function

We know that the density ρ is conserved along these characteristics and thus we can write for any time t: ρ(X(t), t) = ρ(X(s), s) = ρ(x, s).

So in our case, knowing the origin X n of the characteristics, the new value of ρ at t n+1 is given by

ρ n+1 (x) = ρ n+1 (X n+1 ) = ρ n (X n ), ( 25
)
where ρ n is the distribution function at time step t n . The distribution function ρ n is only known on the mesh points, and the origins of the characteristics X n are in general not on a mesh point (see Figure 4). Therefore, we need an interpolation method to compute ρ n at the characteristic's origin, i.e. to approximate ρ n (X n ) needed in the equation [START_REF] Kwon | Development of semi-lagrangian gyrokinetic code for full-f turbulence simulation in general tokamak geometry[END_REF] to get the new value ρ n+1 (x) at the grid points, using the known data on the mesh points at its vicinity. This interpolation method will be either: the quasi-interpolation method using Box-splines, where the Box-splines coefficients, defined in (4), are computed knowing s[k i ] = ρ n (x i ) (see Section 1.2); or the Hermite Finite-Elements interpolation defined in Section 1.3, where the values at the nodes are used and the derivatives or the values at the middle of the edges are computed by finite difference from the values at the nodes.

Localizing the characteristics' origins

One of the advantages of the hexagonal mesh is that it is a uniform mesh. Indeed, even if the mesh is not Cartesian, localizing the characteristics' origin is computationally very efficient, unlike the case of unstructured meshes where iterations are generally required. The procedure is as follows. Let (X 1 , X 2 ) the Cartesian coordinates of the characteristics' origin, obtained by solving [START_REF] Hamiaz | The semi-lagrangian method on curvilinear grids[END_REF]. Then to obtain the hexagonal coordinates (k 1 , k 2 ) of the lowest point of the rhomboid encapsulating the point, we simply need to solve the system x = Rk, where R is the matrix whose columns are the unit vectors given in [START_REF] Abiteboul | Solving the vlasov equation in complex geometries[END_REF], and take the integer value. Denoting by (r ij ) the coefficients of the matrix R, we get

       k 1 = r 22 X 1 -r 12 X 2 r 11 r 22 -r 12 r 21 = 1 √ 3 (X 1 + X 2 ) , k 2 = -r 21 X 1 + r 11 X 2 r 11 r 22 -r 12 r 21 = X 1 + X 2 . ( 26
)
After obtaining (k 1 , k 2 ), we know the rhomboid (composed by two triangular cells) containing the characteristics' origin. To determine the exact cell on which the origin is located, we only need to verify if the abscissa of the point is greater than the abscissa of the mesh point at (k 1 , k 2 ) or not. In the first case the point belongs to the cell on the right, else to the cell on the left.

General algorithm

Below, we summarize the full algorithm to compute the distribution function ρ n+1 solution of the guiding-center model [START_REF] Guscaglia | Finite element methods for the stokes system based on a zienkiewicz type n-simplex[END_REF].

Initialization At time t = 0, we suppose that ρ(x, 0) is given and we evaluate it at the grid points. We denote this data ρ 0 h (meaning ρ discretized, at the time t n = 0, at the initial step). Time Loop Increment a given time step ∆t, such that: t n+1 = t n + ∆t.

• Solve the Poisson equation to compute advection field A n ;

• Compute the characteristics' origins, X n , using an ODE solver for [START_REF] Hamiaz | The semi-lagrangian method on curvilinear grids[END_REF], Runge-Kutta or Adams-Moulton as described above;

• Interpolate (using either Box-splines or one of the Hermite Finite Elements) the distribution function ρ n on X n to compute ρ n+1 h ;

• Update the known values:

ρ n = ρ n+1 .
Remark Boundary conditions will need to be used between the first and the second step of the time loop (i.e. before the interpolation step) for characteristics that leave the computational domain. In this paper we focus only on null Dirichlet boundary conditions.

Numerical results

In this section we present the numerical simulations we performed to test our methods. With the aim of studying the convergence, the dissipation, and the efficiency of the schemes, we first study the circular advection test case. To study the accuracy of the results, we compare them to the analytical solution, which is known. Then we proceed to the guiding-center simulation. As there is no analytical solution for this test case, we study quantities of the system that we know should be conserved.

Circular advection

We focus here on the circular advection test case. The model is defined by:

∂ t f (x, y, t) + y ∂ x f (x, y, t) -x ∂ y f (x, y, t) = 0. ( 27
)
Since this equation is not coupled to a Poisson model, we can study in detail the differences between the interpolation methods previously presented. Additionally, we can find an analytical solution with the method of characteristics thanks to which we can study the convergence of our schemes. Here, we take a Gaussian pulse as initial distribution function:

f 0 (x, y) = exp - 1 2 (x -x c ) 2 σ 2 x + (y -y c ) 2 σ 2 y . ( 28
)
On a hexagonal mesh centered at the origin of radius 8, we take

σ x = σ y = 1 2 √ 2 .
Let us set here x c = 2 and y c = 2. The distance from the pulse to the limit of the domain makes the boundary effects insignificant, thus we can take a null Dirichlet boundary condition, meaning that we consider that there is no inflow to the domain. To study the convergence in space we took N c = 20, 40, 60, ..., 160. We recall that N c is the number of cells on the radius L. With the maximum time of evaluation, t max , at 6π, we chose to keep a constant CFL (here there is not an actual CFL, but we decide to keep the ratio ∆t/∆x constant). In Figure 5, we plotted the L 2 and L ∞ norms for different space discretization. We can see that for coarse meshes, all the methods are globally the same, with a slightly better accuracy for elements Z10 and Ganev-Dimitrov. But as the mesh gets finer, we can quickly see that the splines converge quicker to better results. Only the Ganev-Dimitrov elements are more accurate.

In Figure 6 (left), we represent the error in L ∞ norm versus the number of cells, N c . The figure is similar to the previous one (see Figure 5), but we added the Mitchell elements for comparisons and changed the CFL to CFL = 0.25. The Ganev-Dimitrov element (GaDi) leads again to the best result and the second is Z10. Both have more degrees of freedom (4 times more for GaDi and 3 times more for Z10) for a fixed number of cells N c . Others reconstructions have the same number of degrees of freedom. Furthermore, we notice that the P 1 element (classical linear interpolation) leads to very poor accuracy, which fully justifies the use of higher order methods.

In Figure 6 (right), we selected the order 3 schemes (although we added also Z9, which is not exactly accurate to the order three), and added the third order slope for comparison. On the one hand, we observe a super-convergence property for the splines; this may be due to the test-case (i.e. due to its symmetry) or to the quasi-interpolation that is, in some cases, more accurate than a classical interpolation [START_REF] Condat | Efficient reconstruction of hexagonally sampled data using three-directional box-splines[END_REF]. On the other hand, the Mitchell elements behave favorably: MI4 is as accurate as HCTC and MI6, MI17 are the most accurate, before the splines take over due to this super-convergence, whereas MI6 and MI17 remain third order. We expect higher order convergence when the CFL goes to zero (see [START_REF] Hamiaz | The semi-lagrangian method on curvilinear grids[END_REF]). The Z9 method is the least accurate, as expected. Note also that with respect to Figure 5, the error is bigger for all the methods; being that the number of interpolations done is multiplied by 8, as we go from CFL = 2 to 0.25. All previous figures study the convergence with respect to the hexagonal step or the number of cells N c , thus in Figure 7 we wanted to use the total number of points instead. We can see that the most noticeable difference is that the scheme using Ganev-Dimitrov elements is the worst at the beginning whereas in Figure 6 it is the best one since the beginning. This is probably due to the approximation of the normal derivatives, which requires additional degrees of freedom at the middle of the edges, as we noted in Remark 1.6. Furthermore, after a given mesh refinement, the Box-spline interpolation seems to be the most accurate scheme. In Figure 8, we can see that the performance converges quite quickly, for all the methods. It is also pretty obvious that, even if the splines are more accurate, the cost is higher than most of the Hermite Finite Element methods. The computations where made on a Intel Core with a i7-3630QM CPU with 8 processors at 2.40GHz. In Figure 9 (left), we represent, with respect to N c , the number of points treated per µ-second, also called the efficiency, on the machine irma-hpc2, whose characteristics are: HP Proliant DL 585 G6, 4 processors AMD Opteron 6 Cores at 2.8GHz, 128 Go RAM, 2.4 To of disk. All the algorithms are not at the same level of optimization: HCTC has not been optimized; P1, Z9 and splines have been optimized using different strategies. We think that splines could be further optimized. P1 which is the least costly method is also clearly the fastest method. We obtain around 1GFlops (when the efficiency is about 30) which is a correct value with respect to [START_REF] Mehrenberger | Vlasov on gpu (vog project)[END_REF] for example. It is followed by the Z9 element, and the Mitchell elements which need also the computation of the mixed derivatives. Finally splines and HCTC have almost the same cost; note that previously, on Figure 8, splines were slower; HCTC were also faster, but the runs were not done on the same computer. We then represent on Figure 9 (right), the error in L ∞ norm with respect to time. If GaDi is the most accurate on the finest meshes, the Mitchell elements are among the best in the current implementation. Going higher in degrees leads to time overhead for the computation of the derivatives, this yields that MI17 is not the most competitive, whereas MI6, because it is less accurate, is better, as it is faster. Note as well that Z9 behaves really well. Box-splines are not as competitive, but the situation could change, by optimizing even further the code. We should notice that, for the Mitchell elements, the storage of the different derivatives is multiplied by 13, which leads to costly memory access. Whereas, for Box-splines, this is not the case (coefficients are computed when needed, making this method less costly memory-wise, but slower). Firstly, in Figure 10 (left), we study the performance with classical interpolation schemes on polar mesh (we use cubic splines and Hermite of degree p; see e.g. [START_REF] Hamiaz | The semi-lagrangian method on curvilinear grids[END_REF]). The polar mesh is obtained by discretizing the annulus Ω = {(r cos(θ), r sin(θ)), r ∈ [r min , r max ] × [0, 2π]}, with r min = 10 -5 , r max = 8 using a grid of N c cells in r and 2N c cells in θ direction. In Figure 10 (left), we compare again the number of points advected per µ-second. We see that in this polar setting, the classical cubic splines method has a very good efficiency and the current Hermite interpolation H(p) has a decreased efficiency and this gets worse the degree p increases. On the hexagonal grid, only the P1 interpolation has a better efficiency and Z9 can approach similar efficiency to the splines in polar geometry, while being better than the polar Hermite methods. On the other hand, the efficiency is quite comparable for the Hermite and Mitchell methods for a given degree p. Box-splines are still the less efficient for the moment.

Secondly, we compare the error in L ∞ norm with respect to the CPU time, see Figure 10 (right). Cubic splines in polar geometry outperforms the other methods. Choosing another testcase more favorable to the waste of points due to the polar geometry could change the situation. On the other hand, the second better method with respect to this diagnostic is MI6 , which is encouraging (it is even better than polar splines for very low resolution and it has a more stable and foreseeable behavior). We should again warn the reader that all these studies have been done for a given implementation and that the situation could change if all the methods are fully optimized.

We then look Figure 11 where the L ∞ error is plotted with respect to the number points, making again the comparison between polar and hexagonal mesh. For low number of points, Boxsplines are the least accurate but the situation changes when the number of grid points increases. Note that Box-splines interpolation is not the complete analog of cubic splines interpolation on a polar grid, as it is quasi-interpolant. We can distinguish that MI17 is the most accurate, before Box-splines become better. On the other hand, the accuracy remains quite similar between polar and hexagonal geometry, especially for the Hermite/Mitchell methods. Note also that the cubic splines method is almost the same as H6 (this is often the case, see [START_REF] Steiner | Numerical computation of the gyroaverage operator and coupling with the Vlasov gyrokinetic equations[END_REF]). From this model, and all results extracted from the simulations, we deduce already several conclusions. Firstly, from the comparisons with the P 1 element, we are reminded of the importance of using high order methods. We can also say that the Z9 element should be avoided for this model, and probably, for any other model. Overall, the methods that stood out for the hexagonal mesh, are the Mitchell Element and the Box-splines. We should notice that both schemes can be computed for an arbitrary order, thus yielding more accurate results at a higher cost. When comparing the two methods in terms of accuracy, see Figures 6 and7, splines give slightly better results. While, when we look into performance, notably Figures 9 (left) and 10 (left), Mitchell elements are faster, even when comparing M I17 and Box-splines. The most interesting results come from Figures 9 (right) and 10 (right), where the Mitchell elements are undoubtedly the best interpolation methods for circular advection with our parameters. Still, we strongly believe that for really accurate results, i.e. for finer meshes, the heavy costs of Box-splines would be payed off by their accuracy. In the future, an optimization of the algorithm that computes the Box-splines should be considered. Nevertheless, at the current state of our schemes, our preference is towards the Mitchell elements, such as M I17, and, in second place, towards Box-splines.

Guiding-center model -Diocotron instability test case

We consider here a guiding-center approximation of the 2D Vlasov-Poisson system. This also corresponds to the reduced gyrokinetic model obtained [START_REF] Filbet | Mixed semi-lagrangian/finite difference methods for plasma simulations[END_REF] when all quantities are homogeneous in the direction parallel to the magnetic field. Here the magnetic field is set to B = 0 0 1

T .

Then the model reads

   ∂ρ ∂t + E ⊥ • ∇ x ρ(x, t) = 0 (29a) -∆φ = ∇ • E = ρ(x, t) (29b
)

with E = (E x , E y ) = -∇φ and E ⊥ = (-E y , E x ).
By neglecting the effect of boundary conditions (here, we took null Dirichlet), the guidingcenter model verifies the following properties:

1. Positivity of density ρ 0 ≤ ρ(x, y, t). This model, is commonly used in 2D simulations to study the particle density, as it describes highly magnetized plasmas in the poloidal plane of a tokamak.

We chose here to study the diocotron instability [START_REF] Crouseilles | A new fully two-dimensional conservative semi-lagrangian method: applications on polar grids, from diocotron instability to itg turbulence[END_REF]. The initial density is given by:

ρ 0 (x ⊥ ) =    (1 + ε cos( θ)) exp (-4(r -6.5) 2 ), if r -≤ x 2 + y 2 ≤ r + , with θ = atan2(y, x).
0, otherwise.

(30) Here r + (respectively r -) is the maximum (resp. minimum) radius of an annulus where ρ 0 is not null. θ is the radian angle given by (x, y). As for the parameters values, we take ε = 0.001, r -= 5, r + = 8, = 6, dt = 0.1 and the hexagonal step is 14 160 with a radius of 14 and and a hexagonal parameter N c = 160. In this part, we will not test the Z10 approach as it requires a special resolution of the Poisson equation that has not been implemented. Indeed, computing the values of the field at the center of the triangles can't be combined with the resolution at the vertices. Moreover to even the computational time of each method we chose to take N c = 80 for the Ganev Dimitrov element as it results in the computations on a mesh with N c = 160 (see Section 1.3.4). We notice that 6 vortices is the main mode, even after changing the parameters. In fact, if we take = 6, with ε small enough, we still see the mode 6 appear. With ε big enough, i.e. Figure 13: Potential at time zero for the guiding-center simulation at least 0.1, the modes different from 6 can be visible for a time but they are not stable, thus we see the fusion or the apparition of vortices until the sixth mode takes over. For instance, as illustrated by Figure 12, we can see the ninth mode turning into the sixth mode by fusion of vortices. This instability can be explained with Figure 13. The influence of the geometry is clear as the potential is not round, but already deformed as a hexagon. This phenomenon is clearly caused by the boundary conditions (null Dirichlet) and the shape of the geometry; If other boundary conditions were imposed, we should be able to see results similar to what we can get in a polar mesh (where any given mode can be captured [START_REF] Crouseilles | A new fully two-dimensional conservative semi-lagrangian method: applications on polar grids, from diocotron instability to itg turbulence[END_REF]).

Remark 4.1. When running the guiding-center model with test-case [START_REF] Nguyen | A comparative study of several classical, discrete differential and isogeometric methods for solving poisson's equation on the disk[END_REF], see Figure 12, no obvious differences are visible between the different interpolation methods, which makes the diagnostics all the more important to compare the results computed. specially when comparing the Z9 elements to the splines; the other Hermite elements have globally a worse positivity conservation. We notice that if Box-splines conserve better the mass, the Z9 approach conserves better the L 1 norm. Also we note that Box-splines and the HCTC element give very near results whichever the diagnostic considered.

In Figures 17,18, and 19, we consider the same diagnostics as before, but we include a comparison with the Mitchell elements. Note that the Ganev Dimitrov has not the same degrees of freedom as the other methods. We define, for this method, the total number of degrees of freedom N GaDi tot = 2N c , whereas, for all the other methods N tot = N c . On Figure 17, we consider the case where N c = 128, ∆t = 2 -3 . We then refine the mesh by 2, on Figure 18 and the time step by 2 on Figure 19, in order to see some effects of changing time and space resolutions.

We notice that the Ganev-Dimitrov FE interpolation has good accuracy when we use the same number of cells (e.g. for N c = 128, N tot = 128 and N GaDi tot = 256) but this is not fair as the complexity is not the same. So, we should compare, as done previously, using for example N tot = N GaDi tot = 128 or N tot = N GaDi tot = 256. In that case, the method is no more competitive as the quantities are badly conserved in comparison to the other methods. The situation may be different if the grid is finer, but this then becomes more difficult to solve, as the Poisson solver would take much longer or take too much memory space, at least in the current implementation.

Concerning the mass, the Box-splines method outperforms the other methods (probably due to their property of partition of unity). HCTC has near same accuracy, but it is more oscillating. Z9 is clearly less accurate. We notice that the Mitchell elements are better than Z9 and that increasing p leads to better results. With MI17, we are not at the same level of accuracy as the one for the Box-splines, but we approach it. Note that the results are different on Figure 18, probably because the time step is too big with respect to the resolution in space. This can affect the convergence in time: the solution is more complex and a bigger time step leads to worse mass conservation. Furthermore, this may also explain the fact that MI4 has better mass conservation than MI17. For energy conservation, results are quite similar, mainly, only Z9 behaves better. Energy (as mass) is better conserved, when the time step gets smaller. L 1 norm conservation is improved with the Mitchell elements, in particular with MI17 which is better than Z9; MI6 is better than HCTC and Box-splines, but not than Z9.

We see that L 2 norm is better conserved when increasing the degree p of the Mitchell elements. HCTC, Z9 and MI4 are almost at the same level: slightly better than Box-splines. Note that MI17 is almost at the level of GaDi which requires 4 times more point (see Figure 17). As previously mentioned, this is probably caused by the reconstruction of the normal derivatives (see Remark 1.6). Concerning the L ∞ norm, the best results are obtained for Mitchell's elements, in particular MI17. Box-splines and HCTC do not behave as good as the latter, and Z9 behaves also well. Finally, for the minimum density, HCTC is one of the worst method. Box-splines seem to behave the best. H17 or Z9 come close next, when the time step is small enough, otherwise huge negative values appear for all the methods, except Box-splines.

Conclusion

In this paper we tested Semi-Lagrangian schemes adapted to a hexagonal mesh. The strategies differentiated in the type of interpolation used: on the one hand, we developed an interpolation method based on Box-splines -spline basis specific to the hexagonal mesh-and on the other hand, we introduced an interpolation method based on Hermite Finite Elements. Furthermore, we presented a Poisson solver based on finite differences on this mesh. The first simulations were made on the circular advection test case. This allowed to compare the order of the methods, as well as their efficiency. The splines approach seemed to be more accurate for finer grids, but at higher cost. Whereas, the Mitchell element, was the best method when studying the efficiency. Next, we simulated the 2D guiding-center model. The two methods yield comparable results. The interpolation using Hermite elements is always faster than the Box-splines quasiinterpolation. Box-splines seem to have a better mass and positivity conservation, whereas the Mitchell elements conserved the norms far more accurately and the Z9 seem unbeatable regarding the conservation of energy. Nevertheless, none of the methods stood out from the others. Further research is needed to compare the two methods in more complex problems. Although overall, we believe the splines interpolation method could be optimized and therefore be undoubtedly favored for the circular advection and the guiding-center model when working with fine meshes.

) Remark 1 . 1 .

 11 We notice that uniform B-Splines are Box-splines where the generating matrix is Ξ = h (e 1 , e 2 ), where e 1 = (0, 1) T , e 2 = (1, 0) T and h ∈ R + is the step of the uniform mesh. For a B-spline of degree d, the multiplicity of e 1 and e 2 are both d + 1.

Figure 3 :

 3 Figure 3: On the left: Box-splines χ [r1,r2] , χ [r1,r2,r3] (shifted for better visualization), and χ 2 [r1,r2,r3] (also translated). On the right: 2d projections of the Box-splines onto the x plane.

Algorithm 1 :

 1 Quasi-interpolation with Box-splines Data: Domain denoted Ω. The function f is known at the mesh points x i , sample(i ) = f (x i ). Result: Approximate value of the function f at arbitrary points x (result) forall the x do Initialize result = 0; Initialize coeffs = BoxSplineCoeff(sample, degree); Initialize K = number of points in the vicinity of x i (depends on the type of pre-filter used); Compute hexagonal coordinates: k = [ u v ] where (u v) T = R -1 x ; for k = 0 to K -1 do /* Treat the points on the enveloping rhomboid of radius = degree */ Compute hexagonal coordinates of point in vicinity of x using k -

Figure 4 :

 4 Figure 4: Semi-Lagrangian step: Tracing back characteristics.

Figure 5 :

 5 Figure 5: Errors with L 2 norm (left) and L ∞ norm (right) with CFL = 2

Figure 6 :

 6 Figure 6: Order of convergence (inclusion of Mitchell elements; CFL = 0.25) of all the methods (left) and only for schemes of order 3 (right).

Figure 7 :

 7 Figure 7: Order of convergence (CFL = 0.25) with respect to the number of degrees of freedom

Figure 8 :

 8 Figure 8: Comparison of performances for the circular advection test case (CFL = 2.)

Figure 9 :

 9 Figure 9: Comparison of performances for the circular advection test case (inclusion of Mitchell elements; CFL = 0.25; on machine irmahpc2)

Figure 10 :

 10 Figure 10: Comparison of performances for the circular advection test case and comparison with classical interpolation on polar mesh; CFL = 0.25; on machine irmahpc2

Figure 11 :

 11 Figure 11: Order of convergence (inclusion of Mitchell elements) and comparison with classical interpolation on polar mesh; CFL = 0.25

3 .

 3 L p norm conservation, for 1 ≤ p ≤ ∞ d dt ||ρ|| L p (D)

Figure 12 :

 12 Figure 12: Time evolution of the guiding-center model with ε = 0.1, at time = 1, 16, 38, 73, 109

Figure 14 :

 14 Figure 14: Time evolution of the relative error of mass and energy

Figure 15 :

 15 Figure 15: Relative error of L 1 and L 2 norms

Figure 16 :

 16 Figure 16: Time evolution of the density's minimum

17: Guiding-center case. Evolution of mass, energy, L 1 , L 2 , L ∞ and minimum density

18: Guiding-center case. Evolution of mass, energy, L 1 , L 2 , L ∞ and minimum density

19: Guiding-center case. Evolution of mass, energy, L 1 , L 2 , L ∞ and minimum density

* This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.