Targeting realistic geometry in Tokamak code Gysela

In magnetically confined plasmas used in Tokamak, turbulence is responsible for specific transport that limits the performance of this kind of reactors. Gyrokinetic simulations are able to capture ion and electron turbulence that give rise to heat losses, but require also state-of-the-art HPC techniques to handle computation costs. Such simulations are a major tool to establish good operating regime in Tokamak such as ITER, which is currently being built. Some of the key issues to address more realistic gyrokinetic simulations are: efficient and robust numerical schemes, accurate geometric description, good parallelization algorithms. The framework of this work is the Semi-Lagrangian setting for solving the gyrokinetic Vlasov equation and the Gyseka code. In this paper, a new variant for the interpolation method is proposed that can handle the mesh singularity in the poloidal plane at r=0 (polar system is used for the moment in Gysela). A non-uniform meshing of the poloidal plane is proposed instead of uniform one in order to save memory and computations. The interpolation method, the gyroaverage operator, and the Poisson solver are revised in order to cope with non-uniform meshes. A mapping that establish a bijection from polar coordinates to more realistic plasma shape is used to improve realism. Convergence studies are provided to establish the validity and robustness of our new approach.


Introduction
Understanding and control of turbulent transport in thermonuclear plasmas in magnetic confinement devices is a major goal. This aspect of first principle physics plays a key role in achieving the level of performance expected in fusion reactors. In the ITER design 6 , the latter was estimated by extrapolating an empirical law. The simulation and understanding of the turbulent transport in Fusion plasmas remains therefore an ambitious endeavor.
The Fusion energy community has been engaged in high-performance computing (HPC) for a long time. For example, gyrokinetic simulations are time-hungry (thousands up to millions of CPU-hours) and we then need large amount of computational time that are typically provided by advanced computational facilities. Computer simulation is and will continue to be a key tool for investigating several aspects of Fusion energy technology, because right now there is no burning plasma experiments like ITER. Some of the key issues to address realistic simulations of the Tokamak are: efficient and robust numerical schemes, accurate geometric description, good parallelization algorithms.
The gyrokinetic framework considers a computational domain in five dimensions (3D in space describing a torus geometry, 2D in velocity). Time evolution of the system consists in solving Vlasov equation non-linearly coupled to a Poisson equation (electrostatic approximation, quasineutrality is assumed). The code has the originality to be based on a semi-Lagrangian scheme [20] and it is parallelized using an hybrid OpenMP/MPI paradigm [3,16].
Let z = (r, θ, ϕ, v , µ) be a variable describing the 5D phase space. The time evolution of the ionic distribution function of the guiding-centerf ( z) (main unknown) is governed by the gyrokinetic Vlasov equation (simplified version without right-hand side terms): The guiding-center motion described by the previous Vlasov/transport equation is coupled to a field solver (3D quasi neutral solver which is a Poisson-like solver) that computes the electric potential φ(r, θ, ϕ) (adiabatic electron limit): We will not describe this last equation (details can be found in [8,17]). This Poisson-like equation gives the electric field φ that corresponds to the particle distribution at each time step t. The derivates of J 0 φ along the torus dimensions are computed. Then, these quantities act as a feedback in the Vlasov equation, they appear into the term d z dt B * f . The Vlasov solver represents the critical CPU part, i.e. usually more than 90% of computation time. This equation is solved by splitting it into the advection equations (X G = (r, θ)): B * ∂ tf + ∂ ϕ B * dϕ dtf = 0 (φ operator), B * ∂ tf + ∂ v B * dv dtf = 0 (v operator).
Each advection consists in applying a shift operator along one or two dimensions. A Strang splitting procedure is employed to reach second order accuracy in time. The sequence we choose is (v /2,φ/2,X G ,φ/2,v /2), where the factor 1/2 is a shift over a reduced time step dt/2. In this work, we will propose solutions to improve theX G operator in the Vlasov solver, the gyroaverage J 0 that appears in Eq. (2), and the 2D Poisson equation we need to solve that comes from ρ 2 i ∇ 2 ⊥ eφ Ti term in Eq. (2). These three operators are tightly coupled to the geometry in the poloidal plane which is perpendicular (transverse) to the magnetic field direction. Conversely, thev andφ operators are quite independent from the poloidal geometry because they act in other dimensions thanX G . The paper is organized as follows: in the first section, the original poloidal geometry and meshing is described, the new non-uniform approach focusing the poloidal plane is explained, and the mapping that handles realistic geometry is given. Then, in the second section, the interpolation method on non-uniform polar mesh is investigated, but also advection and gyroaverage operators on such a mesh. Also, we show the numerical method chosen for the 2D Poisson solver. Finally, numerical results and convergence studies are presented in the third section.

New geometry and mapping
Changing the mesh of the poloidal plane while keeping a polar coordinate system should allow us first, to loosen the meshing in order to reduce the typical concentration of points near the center r = 0 and second, to have the mesh match more closely the magnetic surfaces of the plasma. We should then have an improvement in execution time by reducing the overall number of points as well as an improvement in accuracy thanks to the grid being closer to the typical pattern of simulated phenomena. The non-uniform meshing will also allows us to focus on a specific location of the plane that we want to solve by using more points there and only solving roughly elsewhere.

Original polar mesh
We fix N r , the number of points in the radial direction and N θ , the number of points in the poloidal direction. The original polar mesh, as it is defined in Gysela, is such as r i = r min +i∆r with r min > 0, i ∈ 0, N r − 1 , ∆r = rmax−rmin Nr−1 , and also θ j = j 2π N θ with j ∈ 0, N θ − 1 . It is worth noting that r min and r max act as boundary conditions. For each operator that is applied within the poloidal domain, specific ad-hoc approaches are setup to handle what is happening in the central hole r ∈ [0, r min ]. We will not detail the set of ad-hoc boundary conditions that are described in [6].

New non-uniform polar mesh
The new poloidal grid that we want to use is sketched in Figure 1. The idea is to have, for each different circle labeled by r coordinate, a different number of point in the radial direction θ. For instance, in Figure 1 (p. 6), the first layer (inner circle) has four points, the second to fourth layers have eight points and the remaining layers have sixteen points. This allows either to have a density of grid point which is nearly uniform on the plane, or to model finely a subset of the plane which is better solved with more grid points. This meshing or quite similar approaches have already been used in a set of papers [19,12,18]. However, in these previous works the setting and the equations solved were quite different from what we investigate here. Therefore, we have mainly only retained the meshing strategy while redesigning the operators and tools that apply on the mesh.
We have r min = ∆r 2 and r max = r min + (N r − 1)∆r so that ∆r = rmax−rmin Nr−1 , which leads to ∆r = r max N r − 1 2 , and the radial points are Now, for each one of the r i we choose a number of points along θ: N θ [i] , and a grid spacing: , according to what we want to do. Either to focus on a specific region of the plane or to reduce the overall number of points used on the plane and keep the same accuracy.

Mapping
The previous approach can be combined with a general mapping, the polar mapping being only a special case. We focus here on mappings with analytical formula and whose inverse can also be expressed by a formula (to shorten execution time) which was one of the concluding points of [1]. This is of course the case for the polar mapping, but we can also find other more general cases, that can have relevance for the description of the geometry of a tokamak. We consider here the case of a large aspect ratio Tokamak equilibrium, and the mapping that derives from it, as in [5,2]. For the polar mapping x = r cos(θ), y = r sin(θ), the inverse mapping is given by For the large aspect ratio mapping (see [2]; the formula is similar, only ω is changed into π − ω), we have the formula where δ, E, T stand for Shafranov shift, elongation and triangularity. The P notation corresponds to a relabeling of the surfaces. We refer to [2], for the physical interest of such mapping in the tokamaks plasma community. We take here P = 0, T = 0, ω = θ, together with E(r) = E 0 r and δ(r) = δ 0 r 2 ; this clearly restrict the range of geometries, but enables to get an explicit formula for the inverse. We get The inverse mapping can be explicitly given. Putting we are lead to solve x +δ We find Note that r is well defined as soon asδ 0 x + x 2 +ỹ 2 ≤ 1, and when δ 0 = E 0 = R 0 = 0, we recover the polar mapping 1 . We refer to [1] (there, the inverse mapping is also needed) and [11,10] for some works concerning the semi-Lagrangian method combined with a mapping. In the following, we will take Figure 3 (page 6) shows a non-uniform grid combined with this specific mapping.
1 This is not the case for the other solution of the polynomial of degree 2 in r 2 : 2 Operators in complex geometry

Lagrange interpolation in 2D
Let first consider a uniform mesh to introduce the notations, i.e ∀i ∈ 0, [, but also let us define h and k in 0, N r − 1 and 0, N θ − 1 such as r h ≤ r < r h+1 and θ k ≤ θ < θ k+1 where r h = (h + 1 2 )∆r and θ k = k∆θ. Given an order of interpolation p − 1, the Lagrange interpolation polynomial equals m,n are the Lagrange basis polynomial. Then the basis polynomial L (p) m,n (r, θ) associated to the point (r h+m , θ k+n ) reads We can defined a unique set of (β, γ) such as r = r h +β∆r and θ = θ k +γ∆θ, with (β, γ) Then Eq. (6) can be simplified to If the radial position r goes above r max then the coefficient is computed the same way but a Dirichlet condition is used and f (r, θ) is cast to f (r max , θ). If the radial position is located in the interval r ∈ [0, ∆r 2 [, the interpolation scheme has to be adapted because we are crossing the most inner radius of the grid. Let suppose, we have h + m < 0 in Eq. (5), we have to recast r and θ coordinates at the same time to cross the center at r = 0. The new coordinates of a mesh point with h + m < 0 located at (r h+m , θ s ) are set to (r −h−m−1 , θ s+ N θ 2 [N θ ] ). We basically continue the stencil on the radially opposite side of the grid by performing a π rotation.
Considering a non-uniform mesh as described in Section 1, we need to take into account the cases where the interpolation stencil covers several radii as shown in Fig. 2. In that case, the number of points along θ for each radius may be different and so the index of the nearest point in θ direction may be different. For instance on this Figure, the indexes of the interpolation points on radius r 1 at (1.5 ∆r) are 2, 3, 4, 5 (2∆θ [1] , 3∆θ [1] , 4∆θ [1] , . . . ) and on radius r 2 (2.5 ∆r) these indexes are 5, 6, 7, 8 (5∆θ [2] ,6∆θ [2] ,. . . ). In this way, we always use the closest known points, leading to a good accuracy. To adapt the interpolation calculation to non-uniform meshes, with notation of Eq. (5), one only need to first perform the interpolation in θ before the one along r. It allows to easily take into account that the number N θ [i] depends on radius r i . Indeed, if the first interpolation was along r, points along the θ direction would possibly not be available (there is possibly not the same number of points in θ for each radius) and it would require extra 1D interpolations along θ to fix this problem. Algorithm 1 summarizes how the 2D Lagrangian interpolation is performed for non-uniform meshes.

Gyroaverage operator
The gyroaverage operator is a key element in solving the Vlasov-Poisson system of equations, since it allows for the transformation of the guiding center distribution into the actual particle  ∆θ [2] ∆θ [1] ∆θ [h ] /* Poloidal position on radii m */ for n ∈ l, u do distribution, thus reducing the dimensionality of the system of one. The cyclotronic motion of the particles around the magnetic field lines at a distance below the Larmor radius is neglected without loss of accuracy, since this motion is much faster than the turbulence effects usually investigated; moreover, even modern computational power doesn't allow for such highly costly simulations.
A gyroaverage operator has been constructed on the new grid, and we will here briefly describe the numerical implementation which has been adopted in this context.
The gyroaverage operator depending on the spatial coordinates in the polar plane is defined as follows [21]: where x G is the guiding center radial coordinate: it is related to x, the position of the particle in the real space, through the Larmor radius ρ, i.e. x = x G + ρ, which in turn is defined as: where α ∈ [0, 2π] represents the gyrophase angle and e ⊥1 , e ⊥2 the unit vectors of a Cartesian basis in a plane perpendicular to the magnetic field direction b = B/|B|. The function f and g in equation (7) are defined such that f : (r, θ) ∈ R + × R → f (r, θ) is a polar function and g : ( is a Cartesian function such that g(r cos θ, r sin θ) = f (r, θ) for any pair (r, θ). The two functions represent an arbitrary field quantity respectively defined on a grid with polar and Cartesian coordinates. It can be shown [21] that the gyroaverage operator defined in equation (7) can be expressed as a function of the Bessel function of first order, and thus in the Fourier space the gyroaverage is reduced to a multiplication with a Bessel function. In this context though, another approach has been used in order to compute the gyroaverage operator, based on the 2D Lagrangian interpolation. In summary, this method consists in averaging the value of the function over N points equally distributed on a circle of radius ρ: since these points will unlikely correspond to grid points, an interpolation method is used in order to retrieve the value of the function, according to the interpolation procedure described in Section 2.1. This procedure is clarified in Figure 4: the function value for which we want to compute the gyroaverage is marked by an orange circle •, and the red circumference marks the gyroradius which has been considered. Three triangle green points are chosen to compute the gyroaverage, and since they do not correspond to any grid point, the value of the function must be retrieved with a preliminary interpolation, using the nearest grid points available, shown in figure as blue squares gathered around the triangles.
We can write the rigorous expression of the operator in the following way [21]: where α = ∆α, ∆α = 2π/N . P(f ) is the Lagrange interpolator operator. Radial projection on the border of the domain is used if the points selected for the gyroaverage lie outside the domain for large radius. The requirements on the gyroaverage operator are to be accurate enough in order not to disrupt the data, and to be cheap enough from a computational point of view, since it needs to be applied many times during a simulation. It is expected that the present implementation on the new grid will make the application of the gyroaverage operator cheaper and faster, with a general benefit for the global simulation execution time.

Advection operator
Advection consists in the transport of a scalar or vectorial quantity over a vector field. In our case, the transported quantity is the distribution function. The advection is performed backward (Backward Semi-Lagrangian scheme) which means that considering a grid point at time step t N +1 we perform the advection with a velocity field in the opposite direction to find where the quantity was at time step t N (see Figure 14). As the displaced point at time t N seldom corresponds to another grid point, a Lagrange interpolation is performed. The general equation solved by the advection operator for the given distribution function f at point (x, y) is: where v x and v y are the velocities along their respective dimensions and ∆t is the time step. The right-hand side term is solved as explained above by calling the interpolation operator described in 2.1. In Gysela, velocities are defined using a Taylor expansion as described in [7, p. 402].
Time step N Time step N+1

2D finite differences for Poisson solver in polar coordinates
As said in introduction, in a gyrokinetic code the 5D Vlasov equation is coupled to a 3D quasineutrality equation. In Gysela code this last equation is projected in Fourier space in the θ dimension and solved by 1D finite differences in the radial direction. This numerical treatment is well adapted to concentric circular magnetic configuration but will be no longer applicable to more realistic magnetic configuration. Radial and poloidal directions can indeed no more be split and a 2D treatment of the poloidal (r, θ) cross-section is required. A 2D finite element method is often used in the gyrokinetic codes including D-shape magnetic configurations. For the Poisson solver, we will examine two specific meshes: (i) a non-uniform circular mesh (see Figure 1) and (ii) a uniform mesh based on a large aspect ratio equilibrium and mapping (see Section 1.2). We choose to use finite differences to solve this problem. In this section, we consider the 2D Poisson equation in polar coordinates on a domain Ω, with Dirichlet boundary conditions f (r = r max , θ) = g(θ) on ∂Ω.

2D finite differences for a non-uniform circular mesh
Let us first consider equation (10) on a disk Ω = {(r, θ) : 0 < r < r max with r max ∈ R and 0 ≤ θ ≤ 2π} where Ω is described by a non-uniform circular mesh Ω k . To overcome the singularity problem at r = 0, we use the same centered finite difference method as proposed in Lai's paper [14] 2 . One of the trick consists in solving Equation (10) for r ∈ [r min , r max ] with r min = ∆r/2 and a half-integered grid in radial direction and an integered grid in poloidal direction. In this section, we propose an extension of the method proposed for an uniform circular mesh by Lai to a nonuniform one. The difficulty is to adapt the scheme to allow a different number of poloidal points per radius. This implies the adding of interpolations. The scheme proposed in the following is based on Lagrange interpolation of third order which is a good compromise between accuracy and complexity. Let N = N r − 1 be the number of cells in radial direction and N θ [i] be the number of cells along θ on the circle of radius r i . Let us call γ i the ratio between number of poloidal mesh points for r = r i and the one for circle of radius r = r i−1 , namely . Let us add the two constraints on Ω k : (i) N θ [i] is even and (ii) γ i ≥ 1. Then, Ω k is defined as )∆r for all i = 1, 2, · · · , N + 1 (11) where ∆r = r max /(N + 1/2) and ∆θ . Let us notice that these indexes used are different from the one used in section 1.1.2 (indices starting here at 1 instead of 0 previously). Let the discrete values be denoted by ). Then, the discrete version of Eq. (10) becomes, for i = 1, · · · , N and where the boundary values are given: (i) radially by the Dirichlet condition f N +1,j = g j for all is not automatically a mesh point (see Figure 6). The valuef ( The required approximations are computed by using a Lagrange interpolation of third order. So let us consider k ∈ N the integer such that θ k < θ j [i−1] < θ k+1 , then using (5)-(6) notations, where the Lagrange polynomials L n are defined by is computed by Lagrange interpolation of 3rd order by using the mesh points represented by a square.
Let us define, for all i ∈ [1, N + 1], Then, equation (13) reads Let us notice that due to the choice of r 1 = ∆r/2, (1 − λ 1 ) = 0. At the opposite of what is proposed in Lai's paper, let us order the unknowns f i,j radius by radius, such that the unknown The matrix system associated to the discrete equation system (18) The with κ k i = (k − 1)γ i + 1 and L m the Lagrange polynomials defined by Eq. (16) where for more readability, L m (l) = L m (θ j [i] =l ). Let us notice that the column position of the value 1 of first row of C k i is equal to k. Finally, the right hand side vector b can be expressed as Let us notice, that Poisson equation (10) on a circular uniform mesh (N θ

2D finite differences on a mapped uniform mesh
One difficulty was to extend the Poisson solver [14] to a non uniform mesh, as done in the previous subsection. Another one is to deal with a mapping. So, we focus here on this point, starting with a uniform mesh. The combination of both schemes will be the subject of further work and is not tackled here. We refer to [11] for the use of a Mudpack solver, and [10] for the use of a finite element solver based on B-splines. Such solvers might be adapted, but here we consider a specific treatment for the center; so we develop a stand-alone solution with finite differences. Note that some adaptations have to be done with respect to the previous case [14] and we will propose two examples of solvers with 7 and 9 points (we could not get a 5 points solution, here due to the appearance of mixed terms from the mapping as we will see). We consider here the Poisson equation first on a elliptic domain and then for the large aspect ratio mapping (see Eq. (3) and (4) for the latter). We write x(r, θ) = ar cos(θ), y(r, θ) = br sin(θ).
For an ellipse, we have |J| = abr and For the large aspect ratio mapping, we have |J| = (1 − E 0 )br − 2δ 0 br 2 cos(θ), Writing |J|G = (a ij ), we get the equation We consider the following finite difference scheme with 7 points for j = 1, . . . , M and i = 1, . . . , N . We have here f ij = f (r i , θ j ), |J ij | = J(r i , θ j ) and a pq k = a k (r p , θ q ), p, q ∈ 1 2 Z. The system is modified as follows in order to deal with the boundary conditions • u ij is replaced by u ij where j =j + kM, k ∈ Z and 1 ≤j ≤ M .
Note that here a 7 points stencil is needed. We have to take special care on the boundary condition: u 0,j does not cancel and it is replaced by u 1,N/2−j (we assume here that N is even). For the case of a circle, we get the standard 5 points stencil (the terms a 21 , a 2,2 cancel) and u 0,j cancels. Next, we give also another scheme with a 9 points stencil, that is using contributions of u i−1,j−1 and u i+1,j+1 . With respect to [14], we underline that the following adaptations have been done: • a 5 points stencil is (at least seems) no more possible for second order accuracy because of the mixed terms, and several schemes are possible (see [13] p204, [9] p103 and [22] p335 for similar schemes using a stencil with 7 or 9 points) • due to the mixed terms, the term u 0,j does not cancel, we have to use the value u 1,N/2−j , as done in [15].

Interpolation
The interpolation operator is of utmost importance, it is used as a building block by more complex operators. As such, it is essential that this operator remains accurate enough to keep 3 We have the intermediate steps: −3 a i+1/2,j+1/2 + a i−1/2,j−1/2 − a i−1/2,j+1/2 − a i+1/2,j−1/2 u i,j the simulated physics valid. Performance are not detailed in the paper though it is critical and impact almost every piece of the code. It will be presented in future work where it will be integrated in Gysela and compared to previous schemes. The accuracy of the Lagrange interpolation depends on three parameters: the degree of the Lagrange polynomial, the mesh discretization in the r direction and in the θ direction.
In order to test our implementation, we perform interpolations from the polar mesh to a uniform Cartesian grid of size [−r max : r max ] × [−r max : r max ] with 2 N r points in each direction. Points outside of the polar mesh are discarded. The mesh is initialized using a sine product f (x, y) = sin(5x)×cos(4y) for the following tests. Solution is thus analytically known everywhere on the plane. The following figures give different norms (L 1 , L 2 and L inf ) of the error done when performing the interpolation on the whole Cartesian grid. The results are given for the uniform and non-uniform mesh. The base mesh used in the simulation is: N r = 256, N θ = 256 and Lagrange order is 7 for the uniform mesh. For the non-uniform mesh the Lagrange order is also 7, N r = 256 and the N θ [i] are given as such (2 : 32, 8 : 64, 64 : 128, 182 : 256), which reads: there are 32 points in θ direction on the 2 inner most radii (N θ [0] = N θ [1] = 32), 64 points on the 8 following radii, and so on. This gives 15% less points for the non-uniform mesh than for the uniform one with N θ = 256.
In Figure 7, the error is presented against the degree of the Lagrange interpolation which ranges from 1 to 15. Both the uniform and the non-uniform meshes have the same behaviour. There is a convergence phase where the error decreases steadily before stopping at a plateau. It either reaches hardware precision or the accuracy allowed by the meshing on both dimensions. In this case it is hardware precision for a double-precision floating-point (10 −15 ) which is achieved by L inf norm (maximum value of the error). The non-uniform mesh proves to be less accurate because the set of values of N θ [i] we have chosen is good but not optimal.    2 n−1 × 4, 64 : 2 n−1 × 8, 182 : 2 n−1 × 16). On the convergence study in θ direction we find the same difference between the two meshes as shown in Figure 7. The convergence rate is the same which assesses the correctness of the interpolation operator. On Figure 8b the curves perfectly match because the number of points in θ direction is chosen high enough not to influence accuracy (2048 on uniform mesh and ranged from 256 at the center to 2048 at the outter edge on non-uniform mesh). As the meshing method along r has not been changed, both mesh types gives the same convergence results.
On these simple test functions (sine products), using a more complex set of N θ [i] tailored to each specific function allows to reach a reduction of more than a half of the number of points with even fewer accuracy loss compared to a uniform mesh. Whether such N θ [i] set exists for realistic distribution functions is still unknown.

Gyroaverage operator
Numerical results concerning the verification of the implementation of the gyroaverage operator, described in Section 2.2, will now be presented.
A certain family of functions has been considered, whose analytical gyroaverage is known. More precisely, their gyroaverage can be obtained simply by a multiplication of a Bessel function [21]. Given the function f (r, θ) = C m (zr) exp (imθ), where m ≥ 0 is an integer which defines the index of the Bessel function C m (the symbols J m and Y m are used to identify respectively the Bessel functions of the first and the second kind); z ∈ C and (r, θ) represent the usual polar coordinates. The gyroaverage of the function described in Eq. (24) can be written as in [21]: where ρ is the gyroradius, while (r 0 , θ 0 ) define a specific point in the polar mesh.
We have to consider also the boundary conditions. In our case it corresponds to set an homogeneous Dirichlet condition on r max . We can now list the family of functions we used as test cases. It comes directly from the definition (24), and it's written as: where j m, is the -th zero of J m . The function in Eq. 26 is defined on a disk [0, r max ]×[0, 2π] and verifies the Dirichlet boundary condition for which f 1 (r max , θ) = 0, 0 ≤ θ ≤ 2π. The analytical gyroaverage of (26) evaluated at the point (r 0 , θ 0 ) is [21]: The convergence study results will be shown for the first class of functions presented, described by Eq. (26): the order of the Bessel function has been chosen equal to 3, and the first zero has been considered in the argument. A plot of this function can be seen in Figure 9. The uniform and non uniform grid cases have been addressed, and the convergence tests have been performed in both the r and θ directions as well as in the degree of the underlying Lagrange interpolation. Similar tests have been repeated changing the parameters of the test function, namely the order of the Bessel function and the particular zero of the Bessel function chosen: not all of them are shown here, as the results are very close to the ones presented in this section. Figure 10 presents the tests for the convergence in the θ direction: in particular, the logarithm of the L 2 norm, L 1 norm and L ∞ norm are shown with respect to the number of points in the θ direction in a logarithmic scale. The degree of the Lagrange interpolation was fixed to be equal to 5, while the number of points in the r direction and on the gyroaverage circle was respectively equal to N r = 256 and N gp = 128, in order to avoid spurious errors related to these parameters; the gyroradius was set equal to ρ = 0.1. For the uniform grid (described in Section 1.1.1), the number of points in the θ direction was respectively 32, 64, 128 and 256. For the non uniform case however the grid (described in section 1.  For each of these cases, comparing the uniform and non uniform case, we were thus able to reduce the number of points in the θ direction by 22%. For both the uniform and non uniform polar meshes we obtain satisfactory convergence results, in accordance with the theoretical expectations given by the blue dashed curve in the following figures. The theoretical slope is in fact given by n · log(h) + c, where h is the mesh spacing, c an arbitrary constant and the multiplicative factor n in front of the logarithm is given by p + 1, with p being the degree of the Lagrange polynomials used for the interpolation.

Convergence tests
Focusing on the sole logarithm of the L 2 norm, the same convergence test described for the Figure 10b has been repeated for different values of the degree of the Lagrange interpolation polynomial. All the other parameters, namely N r , N θ , ρ and N gp , have been kept as in the setting described for the convergence tests in the non uniform θ-direction. The results of this scan is shown in Figure 11b, where the degree of the interpolation has been changed in the range of (5,7,9). Figure 11a presents the last convergence test which has been performed, namely for the r direction. The logarithm of the L 2 norm, L 1 norm and L ∞ norm are plotted with respect to an increasing number of points in the r direction in a logarithmic scale. Among the parameters kept fixed during the convergence scan in r, the degree of the Lagrange polynomial was set equal to 5, the number of points uniformly distributed in the θ direction was equal to N θ = 504, while the number of points on the gyroaverage circle was set equal to N gp = 128. The gyroradius was still considered to be equal to ρ = 0.1. The number of points in the θ direction had to be chosen uniformly and large enough, in order to avoid a contribution of the error due to discretization in the r direction. The expected theoretical slope, e.g. 6 log(h) + c, due to the degree 5 of the polynomial interpolation used, has been observed in the numerical results.

Uniform and non uniform case comparison
For a conclusive comparison between the uniform and non uniform case, the error in the L 2 norm has been investigated, given the same number of points in the two grids, consequently differently distributed in the space. In particular, the non uniform sequence grid (10 : 32, 30 : 64, 50 :  with a uniform grid with 25 points for each radial direction. The other parameters have been kept fixed during the test, and a sufficient amount of points in the r direction (N r = 256) and for the gyroaverage discretization (N gp = 128) have again been used in order to avoid possible spurious contributions related to these parameters. Among the other parameters, the degree of the Lagrange interpolation was equal to 5 and the gyroradius ρ fixed to 0.1. The results are shown in Figure 12. In Figure 12a we can see that the error is smaller in the uniform case, as expected, since we are restricting the domain region to radii smaller than 0.4. But, if on the other hand we investigate the outer region of the domain (0.4 < r < 0.85), which usually is the one with more need to be accurately resolved due to the interesting physical structures which develop especially here, we can see in Figure 12b that the error is smaller in the non uniform grid case.
Given the results of the convergence tests shown in Figures 10 and 11, we can assess that the gyroaverage operator works satisfactorily on the new grid, thus providing a workable implementation. Considering also the last results presented in Figure 12, we can conclude that the operator is more accurate on the new non uniform grid in the physically interesting region of the domain, thus proving the usefulness of the method proposed.

Advection operator
Several test cases have been performed for the advection operator. The most important is the one testing the resilience of structure when advecting through the center as it is where the mesh is loosen the most. The same meshes as in 3.1 are used (i.e. 15% less points on non-uniform mesh than on uniform mesh). The function studied is the following: otherwise.  where r min ≤ r 1 < r 2 ≤ r max and 0 ≤ θ 1 < θ 2 ≤ 2π. This ensures C 1 continuity for Lagrangian interpolation. The function is then centered on (r = 0.7, θ = 6π 7 ) (r 1 = 0.6, r 2 = 0.8, θ 1 = 5π 7 , θ 2 = π), with r max = 1.0 and then advected through the center with speed (v x = 2 3 , v y = − 1 3 ) until r = 0.7 is reached on the other side. The opposite advection is finally performed to bring back the structure to its original position. The overall displacement is performed in 40 time steps. Initial, middle and final snapshots of the advected function are presented on Figure 13. Middle and final snapshots of the error are presented on Figure 14a and 14c for the uniform mesh and on Figure 14b and 14d for the nonuniform mesh. These show the difference between the exact value and the interpolated one for the whole plane at time steps 1, 20 and 40. During the advection, the error done when interpolating quickly grows when approaching and going through the center (from almost flat error at time step 1 to 14a). But it does not evolve on the outer radii of the mesh, nor when going back through the center. It is more pronounced for the non-uniform mesh. This is best seen in Figure 15 which gives the L 1 and L 2 norms of the error for the whole plane at each time step for both meshes. The structure first undergoes accuracy loss when entering the center near time step 12 for the uniform mesh and near time step 8 for the non-uniform mesh. This means that one of the N θ [i] in that region was not high enough to reach the corresponding accuracy. Radii closer to the center did not have a N θ [i] high enough either as shows the higher overall error for the non-uniform mesh. Once the structure has gone once through the coarsely solved part of the mesh, the accuracy does not undergo any other drastic reduction anymore. It means we have reached the minimum resolution offered by the mesh. This result is satisfactory.
Once again, the operator works as expected. The accuracy highly depends on the choice of the N θ [i] , especially near the center where we wanted the number of points to be scarce. It could be useful to have a tuning tool, which, given the desired accuracy, the number of radii and the typical variation of the function (size of the smallest structure to solve), would give the optimal N θ [i] , but we do not have it yet.
Results with large aspect ratio mapping Results on the large aspect ratio mapping, for the advection are reported on Figures 16,17,18 and 19a,19b. The results are very similar to the previous polar case. We had to change a little the test case, so that the solution does not go outside the domain. The initial condition is now in polar    coordinates:

Reference case with uniform mesh
First tests have been performed on a uniform mesh where the number of poloidal points N θ has been fixed equal to N θ = 2N (N being the radial cell number). The analytic solution is plotted for the case N = 512 on Figure 20 and the corresponding relative error for the numerical solution is shown in Figure 21. The maximum relative error depending on the mesh discretization is summarized in Table 1 for N varying from 8 to 2048. Table 1 also shows that the numerical scheme is as expected of second order. These results act as reference results for the next tests on non-uniform meshes. Let us notice that the results on a uniform mesh are the same than the one obtained by Lai [14]. As said in section 2.4, the matrix system is constructed in the opposite way compared to the one proposed by Lai in order to be able to extend it to meshes with non-uniform number of poloidal points. Both matrix systems have been also compared in terms of CPU time and are equivalent (results for our matrix system are recorded in Table 1). Figure 20: Analytic solution f (r, θ) = exp(r(cos(θ)+sin(θ))) for the case N = 512 and N θ = 1024.

Non-uniform mesh
In this section, tests are now performed on non-uniform meshes defined with N r points in the radial direction and divided into N d sub-domains Ω d depending on the number of points in the poloidal direction N θ[Ω d ] , the maximum number being equal to N θmax = 2(N r − 1) = 2N . Let ν d denotes the reduction of poloidal points for domain Ω d compared to the maximum number, namely N θ[Ω d ] = N θmax /ν d and ν = (ν 1 , · · · , ν N d ). Each sub-domain Ω d of Ω non−uniform is an annulus defined by For all the numerical results presented in this section, the circular domain Ω non−uniform is divided into three sub-domains with ν = [4, 2, 1] (see Figure 6 for example with N r = 8  Table 2 while results for case 2 are detailed in Table 3, both show that the numerical scheme is still of second order. As expected, case 1 is simpler to implement but twice more expensive. As seen in Table 2 and Table 1, CPU time for case 1 is of the same order than for the uniform mesh case because here only the CPU time related to the precomputation (symbolic analysis+factorization) and to the solving are recorded. So calculation of the right hand side with approximation (Lagrange interpolation for missing points) should also be taken into account in case 1. The maximum relative errors obtained for case 1 and case 2 are of the order of the one obtained for a uniform mesh for the smallest number of N θ points. Therefore, as a summary, we have shown in this section that Poisson equation (10) can be solved successfully on a non-uniform circular mesh with the coupling of 2D finite differences in polar coordinates and Lagrange interpolation of third order. The associated numerical scheme proposed in section 2.4 is of second order. Higher order for Lagrange polynomials has not been tested in this paper but the matrix system (20)-(23) could be easily generalized. In both cases, CPU time for the solve step remains quite reasonable even for the biggest meshes. This is a key point for a future implementation in Gysela code where the precomputation step will be only performed one times at the beginning.

Mapping
For validation of the code, we consider the following solution given in cartesian coordinates. We take here a = 1 and b = 0.5 u(x, y) = exp(x) + exp(y) 1 + xy .
We get the results in Tables 4a and 4b Table 5, we give numerical results for f = 2 exp(x + y) and large aspect ratio mapping. We get an order two for the error as expected.
We notice that the error behaves similarly for the two schemes and for the two geometries (ellipse and large aspect ratio). It seems that the results are better for the 7 points scheme in the case of the ellipse and for the 9 points scheme in the case of the large aspect ratio mapping, so that there is no clear trend on what scheme would be the best in general.

Conclusion
In the context of gyrokinetic simulation of turbulence inside a Tokamak plasma, we have developed a strategy that incorporates an adapted non-uniform meshing. Instead of having a circular geometry with a uniform grid along (r, θ) dimensions to describe the poloidal cross-section, we propose a non-uniform spacing along theta direction aiming at computational and memory savings. Additionally, a peculiar mapping coupled with the non-uniform mesh permits to match more complex realistic geometry such as D-shaped plasma, thus exceeding the former limited configuration of circular Tokamak cross-sections. We expect this mapping to rely on analytical formulas in order to keep a relatively low price in term of computations, which is quite crucial for a full-f global code as Gysela is. Several features that are typically used in the gyrokinetic code Gysela have been recast to handle such a new approach. Interpolation, Advection, Gyroaverage and Poisson operators are revisited, upgraded and analyzed in this paper. These operators have been studied separately. All in all, the convergence studies we provide show that this is a workable approach that reach accuracy comparable to uniform meshing. Adaptivity brings really a potential benefit in term of memory savings, but it also permits to refine grid in a specific annulus (small region along radial  Table 5: Convergence for large aspect ratio mapping and 9 points scheme for f = 2 exp(x + y) (fortran implementation with pardiso) direction) that will be useful for physics studies incorporating kinetic electrons.
Future works will target the addition of this method in Gysela. It will require overhauling many data structures and to combine the different operators we have described. As this paper was not focusing on algorithms performance, a subsequent aim will be to write parallel versions of these algorithms and to optimize computation costs in order to compete with the execution time of the former uniform approach. A significant gain in lowering the overall memory footprint is also expected.