Non-uniform numerical schemes for the modelling of turbulence in the 5D GYSELA code

Predicting the performance of a fusion plasma as a function of the plasma power gains, is one of the crucial challenges in fusion plasma physics. With that in mind, turbulence and heat transport must be modelled in a precise theoretical framework which in this case involves the use of 'first principle' non-linear simulations. The 5D (3 spatial coordinates, 2 velocity coordinates) gyrokinetic equations for each species (ions and electrons), coupled with the 3D Maxwell equations form a self-consistent description of the problem. Studies of transport in the core of a tokamak plasma have now reached maturity. Several first principle codes exist which are capable of handling this problem. However, despite their many successes, their prediction capabilities remain limited by the energy content, in particular in the case of discharges with optimised confinement time. In order to push past this limitation, the gyrokinetic models must be extended towards the region at the edge of a tokamak, and, where possible, should treat the transport at the edge and at the core in the same way. The 5D gyrokinetic non-linear GYSELA (V. Grandgirard, Abiteboul, J. Bigot, et al. 2016) code developed by IRFM/CEA is special in that it is global (simulates the entire torus), makes no scale separation approximations ('full F code') and drives turbulence via particle, momentum and heat sources. The additional use of an immersed boundary condition, imitating the extraction of heat by a limiter, makes the GYSELA code one of the rare codes capable of addressing the problem of edge-core turbulence. It is already capable of studying the impact of the edge on the turbulence in the core for electrostatic simulations (i.e. simulations where electrons are considered to be adiabatic). In order to run these simulations, it uses peta-scale high performance calculations (100 million CPU hours/year). The long-term objective for the code is to simulate a turbulent plasma in both the edge and the core with kinetic electrons for the international tokamak ITER which is currently being built at Cadarache (France). We already know that such simulations will not only require tomorrow's exascale resources, but also major numerical changes in the code. This thesis lies within this context and it has a double objective: (i) develop new scalable numerical methods, adapted to the semi-Lagrangian scheme used in the GYSELA code, capable of solving the problem of large fluctuations and temperature variations (1-2 orders of magnitude) at the edge of the plasma, and (ii) take into account more realistic magnetic configurations than the concentric circles currently simulated by the code. Concerning the handling of steep gradients at the edge of the plasma, as the current 5D grids (3 spatial coordinates, 2 velocity coordinates) already represent more than 100 billion points, the proposed solution for adding more points in the edge region is to use a non-equidistant mesh in the radial direction. Even if splines, which have up till now shown themselves to be the best compromise in terms of calculation precision and cost, are retained, adding non-equidistant meshes requires in-depth changes to the compute kernels. From a theoretical perspective, I present a new approach for quadrature using splines, which limits the condition number for the procurement of such quadrature coefficients. One of the disadvantages of splines for their parallelization is their global nature. I present a local spline method where derivatives are transported between patches, and show its stability for semi-Lagrangian advection. From a numerical perspective, bearing in mind that the GYSELA code is a code with more than 50 000 lines, based on a hybrid MPI/OpenMP parallelization, and optimised for more than 100 000 cores, the choice was made to carry out in-depth studies of the semi-Lagrangian method based on non-uniform splines on a model with reduced dimensionality. The problem examined is a Vlasov-Poisson 1D-1V model, used for studies of the plasma sheath, which is a section of the plasma, which presents numerically troublesome steep gradients. The existing VOICE code (which is a mini version of GYSELA), designed to study such problems, has been modified and optimised on a GPU to operate on a non-uniform mesh. These improvements allowed simulations to be carried out which were previously unattainable, and allowed the validation of the semi-Lagrangian method on non-uniform splines. As regards the new more realistic magnetic configurations in the GYSELA code, the choice was made to implement a Culham equilibrium. This equilibrium has the benefit of being based on analytical formulae, while also taking into account the important shaping parameters of an equilibrium plasma: elongation, triangularity, and Shafranov shift. Co-variant and contra-variant transformation matrices were derived and implemented in the code to allow the 5D Vlasov equations to take this geometry into account. The Poisson equation has up to now been solved numerically by projecting each 2D poloidal slice into Fourier space in the periodic poloidal direction, and by using second order finite differences in the radial direction. This solver has been replaced by a 2D finite elements solver based on splines(Zoni and Guclu 2019) which was extracted from the SELALIB library (SeLaLib Development Team 2018). The inclusion of this new magnetic configuration has been successfully numerically validated on the linear benchmarks used for geodesic acoustic mode (GAM) studies. In parallel, a test platform for the 2D Poisson solver was developed in order to numerically compare this spline finite elements solver to two other multi-grid solvers: (i) a solver based on the AMReX library(al. 2019) which uses finite volumes on a uniform cartesian mesh with embedded boundaries, and (ii) a solver developed by the CERFACS which uses finite differences on a logical mesh(Martin J Kuhn, Kruse, and Rude 2022). (author)