N Besse

F Filbet
email: filbet@math.u-strasbg.fr

M Gutnic
email: gutnic@math.u-strasbg.fr

I Paun
email: ipaun@math.u-strasbg.fr

E Sonnendrücker

An adaptive numerical method for the Vlasov equation based on a multiresolution analysis

Introduction

Plasmas, which are gases of charged particles, and charged particle beams can be described by a distribution function f (t, x, v) dependent on time t, on position x and on velocity v. The function f represents the probability of presence of a particle at position (x, v) in phase space at time t. It satisfies the so-called Vlasov equation

∂f ∂t + v • ∇ x f + F (t, x, v) • ∇ v f = 0. ( 1
)
The force field F (t, x, v) consists of applied and self-consistent electric and magnetic fields:

F = q m (E self + E app + v × (B self + B app )),
where m represents the mass of a particle and q its charge. The self-consistent part of the force field is solution of Maxwell's equations

- 1 c 2 ∂E ∂t + ∇ × B = µ 0 j, ∇ • E = ρ ǫ 0 , ∂B ∂t + ∇ × E = 0, ∇ • B = 0.
The coupling with the Vlasov equation results from the source terms ρ and j such that:

ρ(t, x) = q R d f (t, x, v) dv, j = q R d f (t, x, v)v dv.
We then obtain the nonlinear Vlasov-Maxwell equations. In some cases, when the field are slowly varying the magnetic field becomes negligible and the Maxwell equations can be replaced by the Poisson equation where:

E self (t, x) = -∇ x φ(t, x), -ε 0 ∆ x φ = ρ. (2)
The numerical resolution of the Vlasov equation is usually performed by particle methods (PIC) which consist in approximating the plasma by a finite number of particles. The trajectories of these particles are computed from the characteristic curves given by the Vlasov equation, whereas selfconsistent fields are computed on a mesh of the physical space. This method allows to obtain satisfying results with a few number of particles. However, it is well known that, in some cases, the numerical noise inherent to the particle method becomes too important to have an accurate description of the distribution function in phase space. Moreover, the numerical noise only decreases in √ N , when the number of particles N is increased. To remedy to this problem, methods discretizing the Vlasov equation on a mesh of phase space have been proposed. A review of the main methods for the resolution of the Vlasov equation is given in these proceedings [START_REF] Filbet | Numerical methods for the Vlasov equation[END_REF].

The major drawback of methods using a uniform and fixed mesh is that their numerical cost is high, which makes them rather inefficient when the dimension of phase-space grows. For this reason we are investigating here a method using an adaptive mesh. The adaptive method is overlayed to a classical semi-Lagrangian method which is based on the conservation of the distribution function along characteristics. Indeed, this method uses two steps to update the value of the distribution function at a given mesh point. The first one consists in following the characteristic ending at this mesh point backward in time, and the second one in interpolating its value there from the old values at the surrounding mesh points. Using the conservation of the distribution function along the characteristics this will yield its new value at the given mesh point. This idea was originally introduced by Cheng and Knorr [START_REF] Cheng | The integration of the Vlasov equation in configuration space[END_REF] along with a time splitting technique enabling to compute exactly the origin of the characteristics at each fractional step. In the original method, the interpolation was performed using cubic splines. This method has since been used extensively by plasma physicists (see for example [START_REF] Feix | Eulerian codes for the Vlasov equation[END_REF][START_REF] Ghizzo | A Vlasov code for the numerical simulation of stimulated Raman scattering[END_REF] and the references therein). It has then been generalized to the frame of semi-Lagrangian methods by E. Sonnendrücker et al. [START_REF] Sonnendrücker | The Semi-Lagrangian Method for the Numerical Resolution of Vlasov Equations[END_REF]. This method has also been used to investigate problems linked to the propagation of strongly nonlinear heavy ion beams [START_REF] Sonnendrücker | Simulation of heavy ion beams with a semi-Lagrangian Vlasov solver[END_REF].

In the present work, we have chosen to introduce a phase-space mesh which can be refined or derefined adaptively in time. For this purpose, we use a technique based on multiresolution analysis which is in the same spirit as the methods developed in particular by S. Bertoluzza [START_REF] Bertoluzza | An adaptive collocation method based on interpolating wavelets[END_REF], A. Cohen et al. [START_REF] Cohen | Fully adaptive multiresolution finite volume schemes for conservation laws[END_REF] and M. Griebel and F. Koster [START_REF] Griebel | Adaptive wavelet solvers for the unsteady incompressible Navier-Stokes equations[END_REF]. We represent the distribution function on a wavelet basis at different scales. We can then compress it by eliminating coefficients which are small and accordingly remove the associated mesh points. Another specific feature of our method is that we use an advection in physical and velocity space forward in time to predict the useful grid points for the next time step, rather than restrict ourselves to the neighboring points. This enables us to use a much larger time step, as in the semi-Lagrangian method the time step is not limited by a Courant condition. Once the new mesh is predicted, the semi-Lagrangian methodology is used to compute the new values of the distribution function at the predicted mesh points, using an interpolation based on the wavelet decomposition of the old distribution function. The mesh is then refined again by performing a wavelet transform, and eliminating the points associated to small coefficients.

This paper is organized as follows. In section 2, we recall the tools of multiresolution analysis which will be needed for our method, precizing what kind of wavelets seem to be the most appropriate in our case. Then, we describe in section 3 the algorithm used in our method, first for the non adaptive mesh case and then for the adaptive mesh case. Finally we present a few preliminary numerical results.

Multiresolution analysis

The semi-Lagrangian method consists mainly of two steps, an advection step and an interpolation step. The interpolation part is performed using for example a Lagrange interpolating polynomial on a uniform grid. Thus interpolating wavelets provide a natural way to extend this procedure to an adaptive grid in the way we shall now shortly describe.

For simplicity, we shall restrict our description to the 1D case of the whole real line. It is straightforward to extend it to periodic boundary conditions and it can also be extended to an interval with Dirichlet boundary conditions. The extension to higher dimension is performed using a tensor product of wavelets and will be addressed at the end of the section.

For any value of j ∈ Z, we consider a uniform grid G j of step 2 -j . The grid points are located at x j k = k2 -j . This defines an infinite sequence of grids that we denote by (G j ) j∈Z , and j will be called the level of the grid.

In order to go from one level to the next or the previous, we define a projection operator and a prediction operator. Consider two grid levels G j and G j+1 and discrete values (of a function) denoted by (c j k ) k∈Z and (c j+1 k ) k∈Z . Even though we use the same index k for the grid points in the two cases, there are of course twice as many points in any given interval on G j+1 as on G j . Using the terminology in [START_REF] Cohen | Fully adaptive multiresolution finite volume schemes for conservation laws[END_REF], we then define the projection operator

P j j+1 : G j+1 → G j , c j+1 2k → c j k ,
which is merely a restriction operator, as well as the prediction operator

P j+1 j : G j → G j+1 , such that c j+1 2k = c j k , c j+1 2k+1 = P 2N +1 (x j+1 2k+1 ),
where P 2N +1 stands for the Lagrange interpolation polynomial of odd degree 2N + 1 centered at the point (x j+1 2k+1 ). Using the just defined prediction operator, we can construct on G j a subspace of L 2 (R) that we shall denote by V j , a basis of which being given by (ϕ j k ) k∈Z such that ϕ j k (x j k ′ ) = δ kk ′ where δ kk ′ is the Kronecker symbol. The value of ϕ j k at any point of the real line is then obtained by applying, possibly an infinite number of times, the prediction operator.

In the wavelets terminology the ϕ j k are called scaling functions. We shall also denote by ϕ = ϕ 0 0 . Let us notice that

ϕ j k (x) = ϕ(2 j x -k).
It can be easily verified that the scaling functions satisfy the following properties:

-Compact support: the support of ϕ is included in [-2N -1, 2N + 1].
-Interpolation: by construction ϕ(x) is interpolating in the sense that ϕ(0) = 1 and ϕ(k) = 0 if k = 0. -Polynomial representation: all polynomials of degree less or equal to 2N +1 can be expressed exactly as linear combinations of the ϕ j k . -Change of scale: the ϕ at a given scale can be expressed as a linear combination of the ϕ at the scale immediately below:

ϕ(x) = 2N +1 -2N -1 h l ϕ(2x -l).
Moreover the sequence of spaces (V j ) j∈Z defines a multiresolution analysis of L 2 (R), i.e. it satisfies the following properties:

-. . . ⊂ V -1 ⊂ V 0 ⊂ V 1 ⊂ . . . ⊂ V n ⊂ . . . ⊂ L 2 (R). -∩V j = {0}, ∪V j = L 2 (R). -f ∈ V j ↔ f (2 •)V j+1 .
-∃ ϕ (scaling function) such that {ϕ(x -k)} k∈Z is a basis of V 0 and {ϕ j k = 2 j/2 ϕ(2 j x -k)} k∈Z is a basis of V j .

As V j ⊂ V j+1 , there exists a supplementary of V j in V j+1 that we shall call the detail space and denote by W j :

V j+1 = V j ⊕ W j .
The construction of W j can be made in the following way: an element of V j+1 is characterized by the sequence(c j+1 k ) k∈Z and by construction we have

c j k = c j+1 2k . Thus, if we define d j k = c j+1 2k+1 -P 2N +1 (x j+1 2k+1
), where P 2N +1 is the Lagrange interpolation polynomial by which the value of an element of V j at the point (x j+1 2k+1 ) can be computed, d j k represents exactly the difference between the value in V j+1 and the value predicted in V j . Finally, any element of V j+1 can be characterized by the two sequences (c j k ) k of values in V j and (d j k ) k of details in W j . Moreover this strategy for constructing W j is particularly interesting for adaptive refinement as d j k will be small at places where the prediction from V j is good and large elsewhere, which gives us a natural refinement criterion. Besides, there exists a function ψ, called wavelet such that {ψ j k = 2 j/2 ψ(2 j x -k)} k∈Z is a basis of W j . In practise, for adaptive refinement we set the coarsest level j 0 and the finest level j 1 , j 0 < j 1 , and we decompose the space corresponding to the finest level on all the levels in between:

V j1 = V j0 ⊕ W j0 ⊕ W j0+1 ⊕ • • • ⊕ W j1-1 .
A function f ∈ V j1 can then be decomposed as follows

f (x) = +∞ l=-∞ c j0 l ϕ j0 l (x) + j1-1 j=j0 +∞ l=-∞ d j l ψ j l (x),
where the (c j0 l ) l are the coefficients on the coarse mesh and the (d j l ) l the details at the different level in between. In two dimensions, the prediction operator which defines the multiresolution analysis is constructed by tensor product from the 1D operator. In practise three different cases must be considered (see figure 1 for notations):

c j+1 2k 1 +2,2k 2 +1 c j k 1 ,k 2 c j+1 2k 1 ,2k 2 +1 c j+1 2k 1 +1,2k 2 +1 c j+1 2k 1 +1,2k 2 +2 c j k 1 ,k 2 +1 c j k 1 +1,k 2 +1 c j+1 2k 1 +1,2k 2 c j k 1 +1,k 2
1. Refinement in x (corresponding to points c j+1 2k1+1,2k2 and c j+1 2k1+1,2k2+2 ): we use the 1D prediction operator in x for fixed k 2 . 2. Refinement in v (corresponding to points c j+1 2k1,2k2+1 and c j+1 2k1+2,2k2+1 ): we use the 1D prediction operator in v for fixed k 1 .

Refinement in v (corresponding to point c j+1

2k1+1,2k2+1 ): we first use the 1D prediction operator in v for fixed k 1 to determine the points which are necessary for applying the 1D prediction operator in x for fixed k 2 which we then apply.

The corresponding wavelet bases are respectively of type ψ(x)ϕ(v), ϕ(x)ψ(v) and ψ(x)ψ(v) where ϕ and ψ are respectively the scaling function and the 1D wavelet. We then obtain a 2D wavelet decomposition of the following form:

f (x, v) = k1,k2 c j0 k1,k2 ϕ j0 k1 (x) ϕ j0 k2 (v) + j1-1 j0 d row,j k1,k2 ψ j k1 (x) ϕ j k2 (v) + d col,j k1,k2 ϕ j k1 (x) ψ j k2 (v) + d mid,j k1,k2 ψ j k1 (x) ψ j k2 (v) . ( 3
)
3 The algorithms

We want to numerically solve the Vlasov equation ( 1) given an initial value of the distribution function f 0 . We start by describing the method based on an interpolation using the wavelet decomposition of f in the non adaptive case. Then we overlay an adaptive algorithm to this method.

For those two algorithms, we first pick the resolution levels for the phasespace meshes, from the coarsest j 0 to the finest j 1 . Although these levels could be different in x and v, we consider here for the sake of conciseness and clarity that they are identical.

We also compute our scaling function on a very fine grid so that we can obtain with enough precision its value at any point.

The non adaptive algorithm

We are working in this case on the finest level corresponding to j 1 keeping all the points. Initialization: We decompose the initial condition in the wavelet basis by computing the coefficients c k1,k2 of the decomposition in V j0 for the coarse mesh, and then adding the details d j k1,k2 in the detail spaces W j for all the other levels j = j 0 , . . . , j 1 -1. We then compute the initial electric field.

Time iterations:

-Advection in x: We start by computing for each mesh point the origin of the corresponding characteristic exactly, the displacement being v j ∆t.

As we do not necessarily land on a mesh point, we compute the values of the distribution function at the intermediate time level, denoted by f * , at the origin of the characteristics by interpolation from f n . We use for this the wavelet decomposition (3) applied to f n from which we can compute f n at any point in phase space.

-Computation of the electric field: We compute the charge density by integrating f * with respect to v, then the electric field by solving the Poisson equation (this step vanishes for the linear case of the rotating cylinder where the advection field is exactly known). -Advection in v: We start by computing exactly the origin of the characteristic for each mesh point, the displacement being E(t n , x i )∆t. As we do not necessarily land on a mesh point, we compute the values of the distribution function at the intermediate time level, denoted by f n+1 , at the origin of the characteristics by interpolation from f * . We use for this the wavelet decomposition of f * given by (3) used at the previous step.

The adaptive algorithm

In the initialization phase, we first compute the wavelet decomposition of the initial condition f 0 , and then proceed by compressing it, i.e. eliminating the details which are smaller than a threshold that we impose. We then construct an adaptive mesh which, from all the possible points at all the levels between our coarsest and finest, contains only those of the coarsest and those corresponding to details which are above the threshold. We denote by G this mesh.

-Prediction in x: We predict the positions of points where the details should be important at the next time split step by advancing in x the characteristics originating from the points of the mesh G. For this we use an explicit Euler scheme for the numerical integration of the characteristics.

Then we retain the grid points, at one level finer as the starting point, surrounding the end point the characteristic. -Construction of mesh Ĝ: From the predicted mesh G, we construct the mesh Ĝ where the values of the distribution at the next time step shall be computed. This mesh Ĝ contains exactly the points necessary for computing the wavelet transform of f * at the points of G.

Numerical results

We show here our first results obtained with the adaptive method. We consider first a linear problem, namely the test case of the rotating cylinder introduced by Zalesak [START_REF] Zalesak | Fully multidimensional flux-corrected transport algorithms for fluids[END_REF] to test advection schemes. Then we consider a classical nonlinear Vlasov-Poisson test case, namely the two stream instability.

The slit rotating cylinder

We consider the following initial condition: The advection field is (v, -x), which corresponds to the Vlasov equation with an applied electric field E app (x, t) = -x and without self-consistent field. Figure 2 represents the evolution of the rotating cylinder on a half turn with a coarse mesh of 16 × 16 points and 4 adaptive refinement levels. We notice that the cylinder is well represented and that the mesh points concentrate along the discontinuities.

f (0, x, v) = 1 if √ x 2 + v 2 < 0.

The two-stream instability

We consider two streams symmetric with respect to v = 0 and represented by the initial distribution function

f (0, x, v) = 1 √ 2π v 2 exp(-v 2 /2)(1 + α cos(k 0 x)),
with α = 0.25, k 0 = 0.5, and L = 2 π/k 0 . We use a maximum of N x = 128 points in the x direction, and N v = 128 points in the v direction with v max = 7, and a time step ∆t = 1/8. The solution varies first very slowly and then fine scales are generated. Between times of around t ≃ 20 ω -1 p and t ≃ 40 ω -1 p , the instability increases rapidly and a hole appears in the middle of the computational domain. After t = 45 ω -1 p until the end of the simulation, particles inside the hole are trapped. On figure 3 we show a snapshot of the distribution function at times t = 5 ω -1 p and t = 30 ω -1 p for a coarse mesh of 16 × 16 points and 3 levels of refinement. The adaptive method reproduces well the results obtained in the non adaptive case.

Conclusion

In this paper we have described a new method for the numerical resolution of the Vlasov equation using an adaptive mesh of phase-space. The adaptive algorithm is based on a multiresolution analysis. It performs qualitatively well. However, there is a large overhead due to the handling of the adaptive mesh which has not been optimized yet. The performance of the code needs to be improved before we can recommend this technique for actual computations. We are currently working on optimizing the code and trying different kinds of wavelets, as well as obtaining error estimates for the adaptive method.

Fig. 1 .

 1 Fig. 1. Mesh refinement in 2D.

-

  Advection in x: As in the non adaptive case. -Wavelet transform of f * : We compute the c k and d k coefficients at the points of G from the values of f * at the points of Ĝ. -Compression: We eliminate the points of G where the details d k are lower than the fixed threshold. -Computation of the electric field: As in the non adaptive case. -Prediction in v: As for the prediction in x. -Construction of mesh Ĝ: As previously. This mesh Ĝ contains exactly the points necessary for computing the wavelet transform of f n+1 at the points of G determined in the prediction in v step. -Advection in v: As in the non adaptive case. -Wavelet transform of f n+1 : We compute the c k and d k at the points of G from the values of f n+1 at the points of Ĝ. -Compression: We eliminate the points of G where the details d k are lower than the fixed threshold.

  5 and if x < 0 or |v| > 0.125, 0 else. The computational domain is [-0.5, 0.5] × [-0.5, 0.5].

Fig. 2 .

 2 Fig. 2. Rotating cylinder: evolution for a coarse mesh of 2 4 × 2 4 points and 4 adaptive refinement levels. Snapshots of the cylinder and the corresponding adaptive mesh: (upper) after one time step, (lower) after 1/2 turn.

Fig. 3 .

 3 Fig. 3. Two stream instability for a coarse mesh of 2 4 × 2 4 , and 3 adaptive refinement levels, (left) at time t = 5ω -1 p , (right) at time t = 30ω -1 p .