Section: Scientific Foundations

Development of simulation tools

Numerical methods, Vlasov equation, unstructured grids, adaptivity, numerical analysis, convergence, Semi-Lagrangian method The development of efficient numerical methods is essential for the simulation of plasmas and beams. Indeed, kinetic models are posed in phase space and thus the number of dimensions is doubled. Our main effort lies in developing methods using a phase-space grid as opposed to particle methods. In order to make such methods efficient, it is essential to consider means for optimizing the number of mesh points. This is done through different adaptive strategies. In order to understand the methods, it is also important to perform their mathematical analysis. Since a few years we are interested also with solvers that uses Particle In Cell method. This new issue allows us to enrich some parts of our research activities previously centered on the Semi-Lagrangian approach.

Introduction

The numerical integration of the Vlasov equation is one of the key challenges of computational plasma physics. Since the early days of this discipline, an intensive work on this subject has produced many different numerical schemes. One of those, namely the Particle-In-Cell (PIC) technique, has been by far the most widely used. Indeed it belongs to the class of Monte Carlo particle methods which are independent of dimension and thus become very efficient when dimension increases which is the case of the Vlasov equation posed in phase space. However these methods converge slowly when the number of particles increases, hence if the complexity of grid based methods can be decreased, they can be the better choice in some situations. This is the reason why one of the main challenges we address is the development and analysis of adaptive grid methods.

Convergence analysis of numerical schemes

Exploring grid based methods for the Vlasov equation, it becomes obvious that they have different stability and accuracy properties. In order to fully understand what are the important features of a given scheme and how to derive schemes with the desired properties, it is essential to perform a thorough mathematical analysis of this scheme, investigating in particular its stability and convergence towards the exact solution.

The semi-Lagrangian method

The semi-Lagrangian method consists in computing a numerical approximation of the solution of the Vlasov equation on a phase space grid by using the property of the equation that the distribution function $f$ is conserved along characteristics. More precisely, for any times $s$ and $t$, we have

where $\left(𝐗\right(s;𝐱,𝐯,t),𝐕(s;𝐱,𝐯,t\left)\right)$ are the characteristics of the Vlasov equation which are solution of the system of ordinary differential equations

with initial conditions $𝐗\left(t\right)=𝐱$, $𝐕\left(t\right)=𝐯$.

From this property, $f^n$ being known one can induce a numerical method for computing the distribution function $f^{n+1}$ at the grid points $({𝐱}_i,{𝐯}_j)$ consisting in the following two steps:

1. For all $i,j$, compute the origin of the characteristic ending at ${𝐱}_i,{𝐯}_j$, i.e. an approximation of $𝐗(t_n;{𝐱}_i,{𝐯}_j,t_{n+1})$, $𝐕(t_n;{𝐱}_i,{𝐯}_j,t_{n+1})$.

2. As

   $f^{n+1}({𝐱}_i,{𝐯}_j)=f^n(𝐗(t_n;{𝐱}_i,{𝐯}_j,t_{n+1}),𝐕(t_n;{𝐱}_i,{𝐯}_j,t_{n+1})),$

   $f^{n+1}$ can be computed by interpolating $f^n$ which is known at the grid points at the points $𝐗(t_n;{𝐱}_i,{𝐯}_j,t_{n+1}),𝐕(t_n;{𝐱}_i,{𝐯}_j,t_{n+1})$.

This method can be simplified by performing a time-splitting separating the advection phases in physical space and velocity space, as in this case the characteristics can be solved explicitly.

Adaptive semi-Lagrangian methods

Uniform meshes are most of the time not efficient to solve a problem in plasma physics or beam physics as the distribution of particles is evolving a lot as well in space as in time during the simulation. In order to get optimal complexity, it is essential to use meshes that are fitted to the actual distribution of particles. If the global distribution is not uniform in space but remains locally mostly the same in time, one possible approach could be to use an unstructured mesh of phase space which allows to put the grid points as desired. Another idea, if the distribution evolves a lot in time is to use a different grid at each time step which is easily feasible with a semi-Lagrangian method. And finally, the most complex and powerful method is to use a fully adaptive mesh which evolves locally according to variations of the distribution function in time. The evolution can be based on a posteriori estimates or on multi-resolution techniques.

Particle-In-Cell codes

The Particle-In-Cell method [56] consists in solving the Vlasov equation using a particle method, i.e. advancing numerically the particle trajectories which are the characteristics of the Vlasov equation, using the equations of motion which are the ordinary differential equations defining the characteristics. The self-fields are computed using a standard method on a structured or unstructured grid of physical space. The coupling between the field solve and the particle advance is done on the one hand by depositing the particle data on the grid to get the charge and current densities for Maxwell's equations and, on the other hand, by interpolating the fields at the particle positions. This coupling is one of the difficult issues and needs to be handled carefully.

Maxwell's equations in singular geometry

The solutions to Maxwell's equations are a priori defined in a function space such that the curl and the divergence are square integrable and that satisfy the electric and magnetic boundary conditions. Those solutions are in fact smoother (all the derivatives are square integrable) when the boundary of the domain is smooth or convex. This is no longer true when the domain exhibits non-convex geometrical singularities (corners, vertices or edges).

Physically, the electromagnetic field tends to infinity in the neighbourhood of the re-entrant singularities, which is a challenge to the usual finite element methods. Nodal elements cannot converge towards the physical solution. Edge elements demand considerable mesh refinement in order to represent those infinities, which is not only time- and memory-consuming, but potentially catastrophic when solving time dependent equations: the CFL condition then imposes a very small time step. Moreover, the fields computed by edge elements are discontinuous, which can create considerable numerical noise when the Maxwell solver is embedded in a plasma (e.g. PIC) code.

In order to overcome this dilemma, a method consists in splitting the solution as the sum of a regular part, computed by nodal elements, and a singular part which we relate to singular solutions of the Laplace operator, thus allowing to calculate a local analytic representation. This makes it possible to compute the solution precisely without having to refine the mesh.

This Singular Complement Method (SCM) had been developed  [49] and implemented  [48] in plane geometry.

An especially interesting case is axisymmetric geometry. This is still a 2D geometry, but more realistic than the plane case; despite its practical interest, it had been subject to much fewer theoretical studies  [54] . The non-density result for regular fields was proven  [58] , the singularities of the electromagnetic field were related to that of modified Laplacians  [45] , and expressions of the singular fields were calculated  [46] . Thus the SCM was extended to this geometry. It was then implemented by F. Assous (now at Bar-Ilan University, Israel) and S. Labrunie in a PIC–finite element Vlasov–Maxwell code  [47] .

As a byproduct, space-time regularity results were obtained for the solution to time-dependent Maxwell's equation in presence of geometrical singularities in the plane and axisymmetric cases  [70] , [46] .