Nicolas Crouseilles
email: nicolas.crouseilles@inria.fr

Sever A Hirstoaga
email: hirstoaga@math.unistra.fr

Xiaofei Zhao

Multiscale Particle-in-Cell methods and comparisons for the long-time two-dimensional Vlasov-Poisson equation with strong magnetic field

Keywords: Vlasov-Poisson equation, Highly oscillatory, Particle-in-Cell, Numerical comparison, Two-scale formulation, Stroboscopic averaging, Asymptotic preserving Runge-Kutta, Exponential time differencing

Multiscale Particle-in-Cell methods and comparisons for the long-time two-dimensional Vlasov-Poisson equation with strong magnetic field

Nicolas Crouseilles, Sever Adrian Hirstoaga, Xiaofei

Introduction

In this work, we consider the long-time Vlasov-Poisson equation with a strong external homogeneous magnetic field in four phase space dimensions as [START_REF] Frénod | Long time behavior of the two-dimensional Vlasov equation with a strong external magnetic field[END_REF][START_REF] Frénod | Long time behavior of an exponential integrator for a Vlasov-Poisson system with strong magnetic field[END_REF]. The distribution function f ε (t, x, v) depending on time t ≥ 0, on space x = (x 1 , x 2 ) T and on velocity v = (v 1 , v 2 ) T , is then the solution of the following kinetic equation

∂ t f ε + v ε • ∇ x f ε + 1 ε E ε + 1 ε v ⊥ • ∇ v f ε = 0, t > 0, x, v ∈ R 2 , (1.1a)
f ε (0, x, v) = f 0 (x, v), x, v ∈ R 2 , (1.1b)
where the self-consistent electric field E ε (t, x) = -∇ x φ ε (t, x) (with φ ε the electric potential) is determined by the Poisson equation

-∆φ ε = ρ ε -1, ρ ε (t, x) := R 2 f ε (t, x, v)dv, t ≥ 0, x ∈ R 2 . (1.2)
We denote v ⊥ hereafter as

v ⊥ = Jv = (v 2 , -v 1 ) T , with J = 0 1 -1 0 .
The real-valued scalar function f ε (t, x, v) is the unknown representing the space-velocity distribution of electrons. The given parameter ε > 0 is inversely proportional to the strength of the magnetic field. The function ρ ε is referred as the density.

When ε = 1, it corresponds the classical regime and the equation (1.1) has been widely studied in the literature. When the parameter ε becomes very small, the solution f ε (t, x, v) of the Vlasov-Poisson equation (1.1) is known to exhibit fast oscillations in time with wavelength O(ε 2 ). It has been proved that f ε weak- * converges to f (t, x, v) (see [START_REF] Frénod | Long time behavior of the two-dimensional Vlasov equation with a strong external magnetic field[END_REF]) which solves

∂ t f + E ⊥ • ∇ x f + 1 2 ∆φ v ⊥ • ∇ v f = 0, t > 0, x, v ∈ R 2 , (1.3a) E = -∇ x φ, -∆φ(t, x) = R 2 f (t, x, v)dv -1, t > 0, x ∈ R 2 , (1.3b)
f (0, x, v) = f 0 (x, v).

(1.3c)

By integrating (1.3a) with respect to v ∈ R 2 , we further get the well-known Guiding Center equation satisfied by ρ, the limit of the density ρ ε as ε → 0

∂ t ρ + E ⊥ • ∇ x ρ = 0, t > 0, x ∈ R 2 , (1.4a)
E = -∇ x φ, -∆φ(t, x) = ρ(t, x) -1, t > 0, x ∈ R 2 , (1.4b) ρ(0, x) = R 2
f (0, x, v)dv.

(1.4c) Due to high oscillations in time, if one wants to do accurate simulation of the problem (1.1) using classical numerical solvers, one needs to consider small time steps (typically smaller than ε 2 ) which leads to prohibitive time computations in the asymptotic regime. A way to avoid this restriction is to use asymptotic non stiff models like (1.3) or (1.4); such model can be simulated using ε-independent time steps since the models do not contain any stiffness. Obviously, these models are only valid in the asymptotic regime ε 1 and can not correctly describe the dynamic of (1.1) for intermediate regime or when ε ≈ 1. As a consequence, multiscale approaches turn out to be very competitive since they are able to handle small or large values for ε without adapting the numerical parameters. Then, the CPU time is independent of ε. Our goal in this work is to compare, both in terms of accuracy and of CPU time, different multiscale methods recently introduced in the literature for numerically solving highly oscillatory problems like (1.1).

Many multiscale techniques have been devoted to design numerical integrators dedicated to the solving of highly oscillatory problems under different context. For kinetic models in high dimensions, the multiscale schemes are usually performed on the characteristics under the Particle-in-Cell discretization [START_REF] Sonnendrücker | Numerical Methods for Vlasov-Maxwell Equations[END_REF][START_REF] Jacobs | Implicit-Explicit time integration of a high-order particle-in-cell method with hyperbolic divergence cleaning[END_REF]. To tackle the Vlasov-Poisson equation (1.1), a class of semi-implicit Runge-Kutta (SIRK) methods has been proposed in [START_REF] Filbet | Asymptotically stable particle-in-cell methods for the Vlasov-Poisson system with a strong external magnetic field[END_REF]. The SIRK methods could have high order of accuracy and superior stability, and they have been shown in [START_REF] Filbet | Asymptotically stable particle-in-cell methods for the Vlasov-Poisson system with a strong external magnetic field[END_REF] as asymptotic preserving schemes for (1.1) which solve consistently as ε → 0 the Guiding Center limit model (1.4). Very recently, a two-scale formulation (TSF) integrator has been proposed in [START_REF] Crouseilles | Uniformly accurate Particle-in-Cell method for the long time twodimensional Vlasov-Poisson equation with strong magnetic field[END_REF] for solving (1.1). The TSF approach separates the fast variable out in the equation and solves it as an extra degree-of-freedom, which benefits to reduce the oscillation. The TSF integrator has been illustrated in [START_REF] Crouseilles | Uniformly accurate Particle-in-Cell method for the long time twodimensional Vlasov-Poisson equation with strong magnetic field[END_REF] to offer uniform second order accuracy in solving (1.1) for all ε ∈ (0, 1], while extra degreeof-freedom potentially increases the computational cost. An exponential time differencing (ETD) scheme has been proposed in [START_REF] Frénod | Long time behavior of an exponential integrator for a Vlasov-Poisson system with strong magnetic field[END_REF] for solving the Vlasov-Poisson equation with very large time steps with respect to the size of oscillations. The results in [START_REF] Frénod | Long time behavior of an exponential integrator for a Vlasov-Poisson system with strong magnetic field[END_REF] show that the ETD scheme is very accurate as 0 < ε 1, in both short and long-time simulations. On the other hand, under the hierarchy of the Heterogeneous Multiscale Methods [START_REF] Abdulle | The heterogeneous multiscale method[END_REF][START_REF] Engquist | Heterogeneous multiscale methods: a review, Commun. Figure 12: Contour plot of quantity χ ε (t, v) in example (4.3)[END_REF], the stroboscopic averaging method (SAM) has been presented to solve highly oscillatory equations under abstract framework in [START_REF] Calvo | A stroboscopic numerical method for highly oscillatory problems[END_REF]. It has been recently implemented for integrating a highly oscillatory Schrödinger equation in [START_REF] Castella | Stroboscopic averaging for the nonlinear Schrödinger equation[END_REF][START_REF] Ph | Solving highly-oscillatory NLS with SAM: Numerical efficiency and long-time behavior[END_REF], and the success there of SAM in the dispersive equation motivates their prospectiveness in solving highly oscillatory kinetic models.

In this work, we are going to perform the above mentioned methods: SAM, ETD, SIRK and TSF to the Vlasov-Poisson equation (1.1) under the framework of Particle-in-Cell method, and carry out systematical comparisons between them. For SIRK and TSF, we shall briefly present and review them. For ETD and SAM, we shall adapt and derive them in details based on different formulations of the characteristics. We shall test the performance of all the methods for a wide range of ε, and investigate their accuracy and efficiency with comparison carried out aiming at identifying the best one among all in different regimes of the problem.

The rest of the paper is organized as follows. First, we recall the framework of PIC method and derive different formulations of characteristics in Section 2 as preparations for presenting the schemes. The mentioned different kinds of multiscale schemes are given/reviewed in Section 3. Section 4 presents the results of numerical tests and comparisons. Conclusions are drawn in Section 5.

Particle-in-Cell and framework

In this section, we are going to briefly review the Particle-in-Cell (PIC) discretization for the Vlasov-Poisson equation (1.1). We shall first present the PIC framework and then focus on the integrations of the characteristics. We shall derive three different formulations of the characteristics in order to perform different kinds of numerical integrators and obtain corresponding PIC methods.

PIC framework

The PIC discretization consists in approximating the unknown distribution f ε (t, x, v) of (1.1) by a sum of Dirac masses centred at (x k (t), v k (t)) with weight ω k > 0 for k = 1, . . . , N p and N p ∈ N as

f ε p (t, x, v) = Np k=1 ω k δ(x -x k (t))δ(v -v k (t)), t ≥ 0, x, v ∈ R 2 . (2.1)
Plugging (2.1) into (1.1), in the sense of distribution, for k = 1, . . . , N p , one gets that each particle k obeys the characteristics equation

ẋk (t) = v k (t) ε , (2.2a)
vk (t) = E ε (t, x k (t)) ε + v ⊥ k (t) ε 2 , t > 0, (2.2b) x k (0) = x k,0 , v k (0) = v k,0 . (2.2c)
The weight w k and initial values of the particles x k,0 , v k,0 for k = 0, . . . , N p are prescribed according to the given initial distribution f 0 (x, v) in (1.1b). To determine the weight, by integrating (2.1) at t = 0 in whole space we require

Np k=1 ω k = R 2 ×R 2 f 0 (x, v)dxdv.
Thus, a simple choice of uniform weight for all particles would be

w k = 1 N p R 2 ×R 2 f 0 (x, v)dxdv, k = 1, . . . , N p .
The initializing x k,0 , v k,0 for k = 0, . . . , N p is a rather well-established and classical sampling issue. It could be done by the Monte Carlo type rejection sampling method as a general approach. The detailed process can be found in standard statistics textbooks. As a consequence, this approach introduces a noise of order 1/ N p (see [START_REF] Sonnendrücker | Numerical Methods for Vlasov-Maxwell Equations[END_REF]). When f 0 is of variable separation form, i.e. f 0 (x, v) = χ 1 (x)χ 2 (v), especially for the cold plasma case as widely considered in the literature where χ 1 is an one-dimensional distribution and χ 2 is a Gaussian, the initial positions of particles could be given by the so-called inversion of cumulative distribution function. This approach is a deterministic way and a detailed description can be found in [START_REF] Birdsall | Plasma Physics via Computer Simulation[END_REF]. We omit the details here for brevity. The characteristics equations for particles (2.2) are coupled to the Poisson equation (1.2) through the electric field E ε . Once the positions {x k (t)} k=1,...,Np of the particles are obtained at time t > 0, one needs to evaluate the approximated density

ρ ε p (t, x) = Np k=1 ω k δ(x -x k (t)) ≈ ρ ε (t, x) x ∈ R 2 , (2.3)
then solve the Poisson equation

-∆φ ε p (t, x) = ρ ε p (t, x) -1, E ε p (t, x) = -∇φ ε p (t, x), x ∈ R 2 , (2.4)
on a mesh grid of x in R 2 for E ε p (t, x) ≈ E ε (t, x), and finally interpolate for values of E ε (t, x k (t)) at each particle position. In practice, the Dirac delta function δ(x) is approximated by some regularized basis S(x). For example in 1D, it is approximated by the B-spline function S m (x) of order m ∈ N (see [START_REF] Sonnendrücker | Numerical Methods for Vlasov-Maxwell Equations[END_REF]):

S 0 (x) :=    1 ∆x , |x| ≤ ∆x 2 , 0, else, , S m (x) := 1 ∆x x+ ∆x 2 x-∆x 2 S m-1 (y)dy, m ≥ 1. (2.5)
The case in two dimensions is done by tensor product. The B-spline function S m is defined locally in space which is preferred from numerical point of view, but globally it is only a C m-1 function for m ≥ 1. A smooth but global basis has been considered in [START_REF] Jacobs | Implicit-Explicit time integration of a high-order particle-in-cell method with hyperbolic divergence cleaning[END_REF] to get high order accuracy. The classical PIC scheme is hence completed by a standard numerical integrator towards (2.2).

Reformulations of the characteristics equations

In order to design some of the multiscale PIC methods, one needs to reformulate the characteristics equations (2.2). For instance, one could filter out the main oscillations in (2.2) as y k (t) = e -tJ/ε 2 v k (t), t ≥ 0, with e sJ = cos(s) sin(s) -sin(s) cos(s) ,

then one gets a first formulation of the characteristic equations (2.2) well-adapted for the ETD approach

ẋk (t) = 1 ε e tJ/ε 2 y k (t), (2.7a
)

ẏk (t) = 1 ε e -tJ/ε 2 E ε (t, x k (t)), t > 0, (2.7b)
x k (0) = x k,0 , y k (0) = v k,0 . (2.7c)
In order to later apply the SAM and TSF methods, we need to reformulate the latter equations into a suitable form. To do so, we first introduce the following new unknown

u k,+ (t) := x k (t) + εJe tJ/ε 2 y k (t), u k,-(t) := -εJy k (t), t ≥ 0, (2.8)
from which we recover the original ones thanks to

x k (t) = u k,+ (t) + e tJ/ε 2 u k,-(t), and v k (t) = e tJ/ε 2 y k (t) = e tJ/ε 2 ε Ju k,-(t).
Then from (2.7) we get

     uk,+ (t) = JE ε t, u k,+ (t) + e tJ/ε 2 u k,-(t) , uk,-(t) = -Je -tJ/ε 2 E ε t, u k,+ (t) + e tJ/ε 2 u k,-(t) , t > 0.
(2.9)

Keeping in mind that with PIC method, the electric field E ε is given by (2.3) and (2.4), and for simplification we denote

∇ x • E ε (t, x) = Np k=1 ω k δ x -u k,+ (t) + e tJ/ε 2 u k,-(t) -1, (2.10)
where here E ε is always interpreted as a conservative field generated by a potential. If we denote u ± (t) := u 1,± (t), . . . , u Np,± (t) T , as a column vector, then we could interpret the t variable in E ε (t, x) as defined by u ± (t) and the fast time variable t/ε 2 through the nonlinear relation (2.10), i.e. there is a function F ε that reads as

F ε (t/ε 2 , u + (t), u -(t), x) := E ε (t, x) .
Hence, (2.9) can be seen as a system for all k and then can be written as

u± (t) = F ± t/ε 2 , u + (t), u -(t) , (2.11)
where F ± := F 1,± , . . . , F Np,± T and

F k,+ (s, u + , u -) = JF ε s, u + , u -, u k,+ + e sJ u k,-, F k,-(s, u + , u -) = -Je -sJ F ε s, u + , u -, u k,+ + e sJ u k,-. (2.12)
Denoting the initial data in (2.7c) as column vectors x 0 = (x 1,0 , . . . , x Np,0 ) T and y 0 = (v 1,0 , . . . , v Np,0 ) T , we summarize the reformulated characteristics as

u± (t) = F ± t/ε 2 , u + (t), u -(t) , 0 < t ≤ T, u + (0) = x 0 + εJv 0 , u -(0) = -εJv 0 .
(2.13)

This ODE is well-posed on [0, T ], where the final time T > 0 of the above initial value problem is independent of ε and the function F ± (s, u + , u -) is clearly 2π-periodic in s. Hence, (2.13) meets the requirement to perform the multiscale methods like SAM or TSF method. After solving (2.13), one uses

x(t) = u + (t) + e tJ/ε 2 u -(t), v(t) = J ε e tJ/ε 2 u -(t),
to recover the original particle variables.

In the following section, we are going to present PIC schemes by performing different numerical integrators on either (2.2) for SIRK, (2.7) for ETD or (2.13) for SAM and TSF. For simplicity of notation and convenience of comparison, we shall fix all the involved local solvers in the numerical integrators as second order approximations to present the schemes.

Numerical methods

In this section, we present and review the numerical methods we intend to compare in the next section.

Semi-implicit Runge-Kutta

A class of semi-implicit Runge-Kutta (SIRK) schemes have been proposed in [START_REF] Filbet | Asymptotically stable particle-in-cell methods for the Vlasov-Poisson system with a strong external magnetic field[END_REF] for solving the characteristics. These schemes are derived based on (2.2). Among the (second order) SIRK schemes considered in [START_REF] Filbet | Asymptotically stable particle-in-cell methods for the Vlasov-Poisson system with a strong external magnetic field[END_REF], there is a second order L-stable scheme which has been shown to outperform the others in accuracy and stability.

Denoting γ = 1 -1/ √ 2, the second order L-stable SIRK method for the approximation of (2.2) reads

v n+1 = v n + (1 -γ)∆t ε F (1) + γ∆t ε F (2) , x n+1 = x n + (1 -γ)∆t ε v (1) + γ∆t ε v n+1 , (3.1)
where 1) ,

t (1) = t n + ∆t 2γ , x (1) = x n + ∆t 2γε v (1) , v (1) = v n + γ∆t ε F (
F (1) = Jv (1) ε + E ε (t n , x n ), F (2) = Jv n+1 ε + E ε (t (1) , x (1) ).
The above SIRK is uniformly linearly stable in integrating (2.2) and the higher order versions are also derived in [START_REF] Filbet | Asymptotically stable particle-in-cell methods for the Vlasov-Poisson system with a strong external magnetic field[END_REF]. The SIRK (3.1) has been shown in [START_REF] Filbet | Asymptotically stable particle-in-cell methods for the Vlasov-Poisson system with a strong external magnetic field[END_REF] to solve consistently the Guiding Center equation (1.4) in the limit regime, i.e. for any fixed ∆t > 0 and n ≥ 0, x n → x n as ε → 0, where x n is the numerical solution of the characteristics of (1.4) as

x n+1 = x n + (1 -γ)∆tJE(t n , x n ) + γ∆tJE t (1) , x n + ∆t 2γ JE(t n , x n ) .
In this sense, SIRK is asymptotic preserving in the approximations of the particle positions x(t). Thus, by analyzing the SIRK on the limit model (2.2), an error bound can be

|x n -x(t n )| ≤ C(∆t 2 + ε q ), 0 ≤ n ≤ T ∆t , (3.2)
where q appears to be 1 from our the numerical investigations of the convergence of the Vlasov-Poisson equation (1.1) as ε → 0.

Exponential time differencing

An exponential time differencing (ETD) PIC method is proposed in [START_REF] Frénod | Long time behavior of an exponential integrator for a Vlasov-Poisson system with strong magnetic field[END_REF] for the Vlasov-Poisson equation with strong magnetic field under short time scaling. Here we extend the ETD scheme to the long-time scaling problem (1.1) by using the characteristics formulation (2.7).

By applying the Duhamel's principle between t n and t n+1 to the filtered characteristics formulation (2.7), we get

x(t n+1 ) = x(t n ) + tn+1 tn 1 ε e tJ/ε 2 y(t)dt, y(t n+1 ) = y(t n ) + tn+1 tn 1 ε e -tJ/ε 2 E ε (t, x(t))dt, n ≥ 0. (3.3)
For a chosen (macro) time step ∆t = H > 0, we denote

H = N • (2πε 2 ) + rε 2 , N ∈ N, 0 ≤ r < 2π and t * n = t n + N • (2πε 2 ). (3.4)
When H ≥ 2πε 2 , i.e. N > 0, we write the integral equations (3.3) as

x(t n+1 ) = x(t * n ) + J n x , y(t n+1 ) = y(t * n ) + J n y , n ≥ 0,
where

J n x := tn+1 t * n 1 ε e tJ/ε 2 y(t)dt, J n y := tn+1 t * n 1 ε e -tJ/ε 2 E ε (t, x(t))dt,
and we approximate x(t * n ) and y(t * n ) as

x(t * n ) ≈ x(t n ) + N I n x , y(t * n ) ≈ y(t n ) + N I n y , n ≥ 0, (3.5)
where

I n x := tn+2πε 2 tn 1 ε e tJ/ε 2 y(t)dt, I n y := tn+2πε 2 tn 1 ε e -tJ/ε 2 E ε (t, x(t))dt.
It is clear from (2.7) and by fundamental theorem of calculus,

I n x = x(t n + 2πε 2 ) -x(t n ), I n y = y(t n + 2πε 2 ) -y(t n ), J n x = x(t n+1 ) -x(t * n ), J n y = y(t n+1 ) -y(t * n ).
To evaluate x(t n + 2πε 2 ) and y(t n + 2πε 2 ), we solve (2.7) from t n to t n + 2πε 2 , i.e.

ẋ(t) = 1 ε e tJ/ε 2 y(t), ẏ(t) = 1 ε e -tJ/ε 2 E ε (t, x(t)), t n ≤ t ≤ t n + 2πε 2 .
By defining xn (s) := x(t n + ε 2 s) and ỹn (s) := y(t n + ε 2 s), we can rescale the previous differential system to get the following one

       ẋn (s) = εe tnJ/ε 2 e sJ ỹn (s), ẏn (s) = εe -tnJ/ε 2 e -sJ E ε (t n + sε 2 , xn (s)), 0 < s ≤ 2π, xn (0) = x(t n ), ỹn (0) = y(t n ), (3.6)
which is clearly smooth and can be discretized by a classical integrator. Thus, (3.5) can be rewritten as

x(t * n ) ≈ x(t n ) + N (x n (2π) -x(t n )), y(t * n ) ≈ y(t n ) + N (ỹ n (2π) -y(t n )), n ≥ 0. (3.7)
At last, we push (2.7) from t * n to t n+1 to get x(t n+1 ) and y(t n+1 ). By similarly defining the rescaled variables x * (s) := x(t * n + ε 2 s) and ỹ * (s) := y(t * n + ε 2 s), that is say the last step is to solve

       ẋ * (s) = εe t * n J/ε 2 e sJ ỹ * (s), ẏ * (s) = εe -t * n J/ε 2 e -sJ E ε (t * n + sε 2 , x * (s)), 0 < s ≤ r, x * (0) = x(t * n ), ỹ * (0) = y(t * n ), (3.8)
which is a smooth problem as well and easy to discretize. Hence, at each time level t n , the main ETD scheme is given by (3.7) associated with a local initial-value problem (3.6), and accomplished by another local initial-value problem (3.8) which aims at offering

x(t n+1 ) = x * (r), y(t n+1 ) = ỹ * (r).
For detailed discretization of the ETD scheme, when N > 0 in (3.4), we choose an integer n 0 > 0 to discretize (3.6) and (3.8) in time and denote the (micro) time step h = 2π n0 , h r = r n0 , then the scheme goes as

x * = x n + N S h x n0 (x n , y n ) -x n , x n+1 = S hr x n0 (x * , y * ), (3.9a
)

y * = y n + N S h y n0 (x n , y n ) -y n , y n+1 = S hr y n0 (x * , y * ), n ≥ 0, (3.9b)
where (S h x ) n0 (x n , y n ) ≈ xn (2π) and (S h y ) n0 (x n , y n ) ≈ ỹn (2π) are the numerical solution of (3.6) with initial data (x n , y n ) by a numerical integrator under time step h. (S hr x ) nr (x * , y * ) and (S hr y ) nr (x * , y * ) are defined similarly for (3.8) with initial data (x * , y * ) discretized by time step h r . Here, we choose the numerical integrator S h

x , S h y as the central difference integrator for practical implementation, i.e. solving u = λu with

u 1 = u 0 + hλu 0 , u n+1 = u n-1 + hλu n , n ≥ 1.
When N = 0 in (3.4), that is the case if H < 2πε 2 , we reduce the ETD scheme to a classical central-difference method towards (2.7) with only the macro time step H. In this convention, the ETD scheme is well-defined for all H and ε.

The central difference integrator is stable on the [-1, 1] interval of the imaginary axis. We remark the reason why we choose this 'mean' local integrator rather other more robust integrators such as A-stable schemes, is because we want to investigate the robustness of the whole approach.

The ETD scheme is fully explicit. For solving (2.7), i.e. a single particle, the computational cost of ETD at each (macro) time level is of order O(n 0 ) where n 0 is the number of time grids in micro step, and the memory cost is O(1). The optimal error bound or rigorous error estimates of the ETD scheme is still not clear to us. Nevertheless, the numerical results indicate it is second order accurate in the classical regime and is first order accurate in the limit regime.

Stroboscopic averaging method

The stroboscopic averaging method (SAM) is introduced systematically in [3] under general framework. By a rescaling t = ε 2 s and w ± (s) = u ± (t), the model (2.13) becomes

ẇ± (s) = ε 2 F ± (s, w + (s), w -(s)), 0 < s ≤ T ε 2 . (3.10)
As shown in [START_REF] Calvo | A stroboscopic numerical method for highly oscillatory problems[END_REF][START_REF] Castella | Stroboscopic averaging for the nonlinear Schrödinger equation[END_REF], through a smooth change of variable Φ ± (s, Ψ + , Ψ -) which is 2π-periodic in s, there exists a flow Ψ ± and an averaged field (3.11) such that for some constant C > 0 independent of ε nor s

G ± satisfying    Ψ± (s) = ε 2 G ± (Ψ + (s), Ψ -(s)), 0 < s ≤ T ε 2 , Ψ ± (0) = w ± (0),
|w ± (s) -Φ ± (s, Ψ + (s), Ψ -(s))| ≤ Ce -C/ε 2 , 0 ≤ s ≤ T ε 2 .
For time at the multiples of period, as known as the stroboscopic times, we have the stroboscopic property,

Φ ± (s, Ψ + , Ψ -) = Ψ ± (s), s/2π ∈ N.
Hence solving the averaged system (3.11) becomes a good choice to approximate (3.10) at the stroboscopic time. As derived in [START_REF] Castella | Stroboscopic averaging for the nonlinear Schrödinger equation[END_REF][START_REF] Ph | Solving highly-oscillatory NLS with SAM: Numerical efficiency and long-time behavior[END_REF], the unknown function G ± is approximated by

G ± (Ψ + , Ψ -) = 1 4πε 2 (S ± (2π, Ψ + , Ψ -) -S ± (-2π, Ψ + , Ψ -)) + O(ε 4 ), (3.12)
or

G ± (Ψ + , Ψ -) = 1 24πε 2 (-S ± (4π, Ψ + , Ψ -) + 8S ± (2π, Ψ + , Ψ -) -8S ± (-2π, Ψ + , Ψ -) + S ± (-4π, Ψ + , Ψ -)) + O(ε 8 ), (3.13)
where for some Ψ ± and P > 0, S ± (P, Ψ + , Ψ -) is the solution of

ẇ± (s) = ε 2 F ± (s, w + (s), w -(s)), 0 < s ≤ P, w ± (0) = Ψ ± , (3.14)
and S ± (-P, Ψ + , Ψ -) is the solution of

ẇ± (s) = ε 2 F ± (s, w + (s), w -(s)), -P ≤ s < 0, w ± (0) = Ψ ± , (3.15)
Then a numerical approximation to (2.13) can be constructed by solving the macro-dynamics which is a reverse rescaling on (3.11):

˙ Ψ ± (t) = G ± ( Ψ + (t), Ψ -(t)), 0 < t ≤ T, Ψ ± (0) = u ± (0), (3.16)
under the help of the micro-dynamics (3.14)-(3.15) for obtaining G ± from (3.12) or (3.13). Assume T /ε 2 ∈ 2πN, choose a macro time step ∆t = H > 0 and a micro step h for discretizing the time domains N 0 = T H , n 0 = 2π h , and denote the numerical solution to

u n ± = u ± (t n ) ≈ Ψ(t n ) for some t n = nH, 0 ≤ n ≤ N 0 .
A central-difference time discretization for both macro part and micro part reads: choose u 0 ± = u ± (0) and

u 1 ± = u 0 ± + HG 0 ± , u n+1 ± = u n-1 ± + 2HG n ± , n = 1, . . . , N 0 , (3.17)
where

G n ± = 1 4πε 2 (S h ± ) n0 (u n + , u n -) -(S -h ± ) n0 (u n + , u n -) , n ≥ 0, (3.18)
or

G n ± = 1 24πε 2 -(S h ± ) 2n0 (u n + , u n -) + 8(S h ± ) n0 (u n + , u n -) -8(S -h ± ) n0 (u n + , u n -) + (S -h ± ) 2n0 (u n + , u n -) , (3.19)
and (S h ± ) n0 (u n + , u n -) is the numerical flow of (3.14) with P = 2π, Ψ ± = u n ± by the central-difference discretization with time step h, and (S -h ± ) n0 (u n + , u n -) is similarly the numerical flow of (3.15) with time step -h.

We call the SAM (3.17) with (3.18) as SAM2, and address the SAM (3.17) with (3.19) as SAM4 for short. As studied in [START_REF] Ph | Solving highly-oscillatory NLS with SAM: Numerical efficiency and long-time behavior[END_REF][START_REF] Calvo | A stroboscopic numerical method for highly oscillatory problems[END_REF], with a given smooth enough function F ± in (3.10), SAM2 for integrating (3.10) gives an error bound

|u ± (t n ) -u n ± | ≤ C(H 2 + ε 2 h 2 + ε 4 ), 0 ≤ n ≤ T H , (3.20)
and SAM4 gives an error bound

|u ± (t n ) -u n ± | ≤ C(H 2 + ε 2 h 2 + ε 8 ), 0 ≤ n ≤ T H . (3.21)
At each (macro) time level, the memory cost of SAM2 and SAM4 is O(1) and the computational cost is O(n 0 ) where n 0 is the number of time grids in micro step. SAM4 is twice more expensive than SAM2.

Two-scale formulation

The two-scale formulation (TSF) method is a general approach to tackle highly oscillatory problems [START_REF] Ph | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF]. It has recently been developed to solve the kinetic models [START_REF] Crouseilles | Uniformly accurate forward semi-Lagrangian methods for highly oscillatory Vlasov-Poisson equations[END_REF][START_REF] Crouseilles | Uniformly accurate Particle-in-Cell method for the long time twodimensional Vlasov-Poisson equation with strong magnetic field[END_REF] on characteristics.

To solve (1.1) under PIC, the two-scale formulation is performed for (2.13) as

∂ t U ± (t, τ ) + 1 ε 2 ∂ τ U ± (t, τ ) = F ± (τ, U ± (t, τ )), t > 0, τ ∈ T, (3.22)
with initial data satisfying

U + (0, 0) = x(0) + εJy(0) =: u 0 + , U -(0, 0) = -εJy(0) =: u 0 -. (3.23)
The augmented equation (3.22) recovers the solution of (2.7) via

x(t) = U + (t, t/ε 2 ) + e tJ/ε 2 U -(t, t/ε 2 ), y(t) = 1 ε JU -(t, t/ε 2 ).
The initial data of the augmented problem (3.22) is not fully determined. By using the Chapman-Enskog expansion, one can decompose the problem (3.22) into a macro part and a micro part, and a well prepared initial data in order to bound the time derivatives of U ± (t, τ ) could then be constructed. To present the initial data, we introduce for some periodic function u(τ ) on T,

Πu = 1 2π 2π 0 u(s)ds and L -1 u(τ ) = (I -Π) τ 0 u(s)ds if Πu = 0.
Denote the operator A := L -1 (I -Π). The initial data is constructed in spirit of a two-scale PIC approximation as follows. One starts with X 0th (τ ) := u 0 + + e τ J u 0 -. Using X 0th (τ ) through (2.10), define

-∇ x • E 0th (τ, x) = Np k=1 ω k δ(x -X 0th k (τ )) -1. (3.24)
With E 0th , one can define the 'well-prepared' initial data U 1st ± (τ ) as

U 1st ± (τ ) := u 0 ± + h 0th ± (τ ) -h 0th ± (0), h 0th ± (τ ) := ε 2 AF 0th ± (τ, u 0 ± ), (3.25)
where

F 0th + (τ, u 0 ± ) := JE 0th τ, u 0 + + e Jτ u 0 -, F 0th -(τ, U 1st ± ) := -Je -Jτ E 0th τ, u 0 + + e Jτ u 0 -.
We refer the readers to [START_REF] Crouseilles | Uniformly accurate Particle-in-Cell method for the long time twodimensional Vlasov-Poisson equation with strong magnetic field[END_REF] for the detailed derivation. For t ≥ 0, the two-scale electric field is always computed from (2.10). As a property of using U 1st ± as the initial data for the two-scale problem (3.22), we can have [START_REF] Ph | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF][START_REF] Crouseilles | Uniformly accurate Particle-in-Cell method for the long time twodimensional Vlasov-Poisson equation with strong magnetic field[END_REF]

∂ k t U ± (t, τ ) = O(1), ∂ k t F ± (t, τ, U ± ) = O(1), k = 1, 2, 0 < ε 1.
As a consequence, we have

d k dt k F ± (t, τ, U ± (t, τ )) = O(1), k = 1, 2, 0 < ε 1. (3.26)
An exponential integrator is then proposed to integrate the two-scale equation (3.22) with the constructed initial data (3.25).

The detailed numerical scheme reads: Denote U n ± (τ ) ≈ U ± (t n , τ ) for n ≥ 0 and let U 0 ± = U (0, τ ). We update the U n ± for n ≥ 1 as

(U ± ) 1 = e -il∆t ε 2 (U ± ) 0 + p (F ± ) 0 + q 1 ∆t (F ± ) * ,1 -(F ± ) 0 , (3.27a)
(U ± ) n+1 = e -i ∆t ε 2 (U ± ) n + p (F ± ) n + q 1 ∆t (F ± ) n -(F ± ) n-1 , n ≥ 1, (3.27b)
where

U n ± (τ ) = Nτ /2-1 =-Nτ /2 (U ± ) n e i τ , F n ± (τ ) = Nτ /2-1 =-Nτ /2 (F ± ) n e i τ , F * ,1 ± (τ ) = Nτ /2-1 =-Nτ /2 (F ± ) * ,1
e i τ , n ≥ 0, and

F n + (τ ) = JE n τ, U n + (τ ) + e Jτ U n -(τ ) , F n -(τ ) := -Je -Jτ E n τ, U n + (τ ) + e Jτ U n -(τ ) , n ≥ 0, F * ,1 + (τ ) = JE * ,1 τ, U * ,1 + (τ ) + e Jτ U * ,1 -(τ ) , F * ,1 -(τ ) := -Je -Jτ E 1, * τ, U * ,1 + (τ ) + e Jτ U * ,1 -(τ ) , with (U ± ) * ,1 = e -i ∆t ε 2 (U ± ) 0 + p (F ± ) 0 , U * ,1 ± (τ ) = Nτ /2-1 =-Nτ /2 (U ± ) * ,1 e i τ ,
and

p := ∆t 0 e -i ε 2 (∆t-s) ds =    iε 2 e -i ∆t ε 2 -1 , = 0, ∆t, = 0, q := ∆t 0 e -i ε 2 (∆t-s) sds =        ε 2 2 ε 2 -ε 2 e -i ∆t ε 2 -i ∆t , = 0, ∆t 2 2 , = 0.
The numerical electric field E n is given by collecting all the two-scaled characteristics {U n ±,k } Np k=1 and solving the Poisson equation:

E n (τ, x) = -∇ x Φ n (τ, x), -∆Φ n (τ, x) = Np k=1 ω k S(x -X n k (τ )) -1, X n k (τ ) = U n +,k (τ ) + e τ J U n -,k (τ ), E * ,1 (τ, x) = -∇ x Φ * (τ, x), -∆Φ * (τ, x) = Np k=1 ω k S(x -X * k (τ )) -1, X * k (τ ) = U * ,1 +,k (τ ) + e τ J U * ,1 -,k (τ ), (3.28)
where S(x) is the regularised basis.

The proposed TSF scheme (3.27)-(3.28) is self-consistent for solving (3.22). That is to say during the computation, the scheme runs with U n ± (τ ) only. It does not require approximations of the original components x(t), y(t) for (2.7) or X(t, τ ), Y (t, τ ). Moreover, the scheme is fully explicit and efficient thanks to fast Fourier transform in τ . It has been shown in [START_REF] Crouseilles | Uniformly accurate Particle-in-Cell method for the long time twodimensional Vlasov-Poisson equation with strong magnetic field[END_REF] that the TSF scheme (3.27)-(3.28) for solving (3.22) with the 'well-prepared' initial data U 1st ± offers uniformly second order accuracy in the time discretization for all 0 ≤ ε ≤ 1. The finite time error bound of the TSF for solving (2.10) formally goes as

|u ± (t n ) -U n ± (t n /ε 2 )| ≤ C(∆t 2 + N -m0 τ ), 0 ≤ n ≤ T ∆t ,
where m 0 > 0 depends on the smoothness of the solution of (3.22). The memory cost is O(N τ ) and computational cost per time step is O(N τ log N τ ). With {U n ±,k } Np k=1 from the scheme (3.27)-(3.28), one can obtain the approximation to the solution of the original characteristics (2.2) as

x k (t n ) ≈ X n k t n ε 2 , v k (t n ) ≈ J ε e J tn ε 2 U n -,k t n ε 2 , n ≥ 1.
Then together with the PIC approximation (2.1), the numerical scheme of a TSF method for solving the Vlasov-Poisson equation (1.1) is completed. We remark that the reason why we choose to solve (3.22) by an exponential integrator rather a finite difference integrator such as the central difference scheme used in [START_REF] Ph | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF][START_REF] Crouseilles | Asymptotic preserving schemes for highly oscillatory Vlasov-Poisson equations[END_REF], is because the error bound of a second order finite difference discretization depends on the norm of the third order time derivative ∂ 3 t U ± , whereas for the exponential integrator, the error depends on the second order time derivative ∂ 2 t U ± (see [START_REF] Crouseilles | Uniformly accurate Particle-in-Cell method for the long time twodimensional Vlasov-Poisson equation with strong magnetic field[END_REF]). However, the constructed initial data (3.25) is not enough to ensure a uniform bound on the third order time derivative term, and one needs to construct a higher order 'well-prepared' initial data which is more involved.

Numerical results

In this section, we shall conduct numerical experiments to test the accuracy, efficiency and long-time behavior of all the mentioned numerical methods for solving the Vlasov-Poisson equation (1.1) with a wide range of ε ∈ (0, 1]. Comparisons between different methods shall be made. Throughout this section, we shall use the following numerical example unless specified.

We consider the Vlasov-Poisson equation (1.1) with the initial data of a Kelvin-Helmholtz instability type [START_REF] Crouseilles | Conservative semi-Lagrangian schemes for Vlasov equations[END_REF][START_REF] Shoucri | A two-level implicit scheme for the numerical solution of the linearized vorticity equation[END_REF] in space and a double-Maxwellian distribution in velocity space:

f 0 (x, v) = 1 4π (1 + sin(x 2 ) + η cos(kx 1 )) e -(v 1 +2) 2 +v 2 2 2 + e -(v 1 -2) 2 +v 2 2 2
.

(4.1)

The computational domain for x is Ω = [0, 2π/k] × [0, 2π] for some k, η > 0. We choose η = 0.05, k = 0.5 and discretize Ω with 64 points in x 1 -direction and 32 points in x 2 -direction. The Poisson equation (1.2) is always solved on Ω with periodic boundary condition by the fast Fourier transform. We fix the basis function S(x) for the PIC method as the fifth order B-spline, i.e. m = 5 in (2.5). The high order spline function is chosen here to reach a high accuracy in spatial approximation so that the temporal error of each method can be distinguished. We take N p = 204800 and generate the initial position and velocity of each particle by means of the rejection sampling. All the numerical methods are programmed in a sequential way in Fortran and are run on an iMac with processor 2.5Ghz Intel Core i5.

Accuracy test

We first test the short time temporal accuracy of the PIC method coupled with SIRK (3.1), ETD (3.9), SAM2 (3.18), SAM4 (3.19) and TSF (3.27). We solve the Vlasov-Poisson equation (1.1) with (4.1) by the numerical methods with the setup in PIC as described before till T = π 2 . We compute the error of each method regarding the density

ρ ε (t, x) = R 2 f ε (t, x, v)dv,
x ∈ Ω, and the error regarding the second moment in v

ρ ε v (t, x) = R 2 |v| 2 f ε (t, x, v)dv, x ∈ Ω,
where a reference solution is obtained numerically via the TSF with very small time step, e.g. ∆t = 10 -3 and N τ = 64. We measure the error in ρ ε and ρ ε v at t = T = π 2 under maximum norm in space Ω. For ETD (3.9), SAM2 (3.18) and SAM4 (3.19) methods, we denote by n 0 the number of grid points (in time) for the micro-solver and define it as n 0 = T /∆t, and hence ∆t = H = O(h). The errors of SIRK, ETD, SAM2 and SAM4 under different ∆t and ε are shown in Figures 1, 3, 4 and5, respectively. Note that only the error in ρ ε is shown in Figure [START_REF] Abdulle | The heterogeneous multiscale method[END_REF] for the SIRK method; indeed, the method is designed so as to capture the slow dynamic through the Guiding Center model, then plotting the error in ρ ε v is not relevant. The results of TSF under N τ = 32 are shown in Figure 6.

In order to understand the asymptotic error bound of the SIRK scheme in the limit regime and then to have an estimate of the q exponent in (3.2), we investigate the convergence rate of the Vlasov-Poisson equation (1.1) to (1.3) as ε → 0. To do so, we solve the Vlasov-Poisson equation (1.1) accurately by TSF, and solve the limit model (1.3) accurately by classical PIC with fourth order Runge-Kutta scheme. The convergence results in the two quantities ρ ε and ρ ε v as ε → 0 are shown in Figure 2 which indicate the rate q is one. Next, in order to the test convergence of the TSF method in τ -direction, we show the error of the TSF method under different N τ in Figure 7.

Based on the convergence results in Figures 1234567, we can draw the following conclusions: Concerning SIRK method, the second order accuracy is well observed in the classical regime ε = O(1), while in the intermediate regime of ε this scheme loses convergence. The reasons are from two sides. One is due to the temporal oscillations, and the other is due to the formulation (2.2) where the approximation error of E ε (computed numerically through Poisson equation) is amplified by the term 1/ε. As ε becomes very small, the SIRK scheme recovers some convergence rate in the density function ρ ε but the error always scales as O(ε) due to the convergence rate of the model as explained before thanks to Figure 2. Indeed, we see the line in the left figure of Figure 2 matches with the slope in the right figure of Figure 1.

Concerning ETD, second order accuracy is also observed for ε = O(1), but order reduction occurs (first order) in the limit regime 0 < ε 1, which is consistent with the observation in [START_REF] Frénod | Long time behavior of an exponential integrator for a Vlasov-Poisson system with strong magnetic field[END_REF]. In the intermediate regime of ε, the method has strong convergence order reduction. There is a heavy bump in the error when the ETD reduces to classical integrator, but it is accurate in both ρ ε and ρ ε v if the time step is small enough. Concerning SAM methods, the errors of SAM2 and SAM4 stuck at O(ε 4 ) and O(ε 8 ) respectively. This is consistent with the fact that these methods only solve the averaged model which is valid when ε is small. Thus, they do not work for ε = O(1), but in the limit regime when O(ε 4 ) and O(ε 8 ) are negligible, the SAM methods can have second order convergence rate in time and offer very accurate approximation. In the intermediate regime of ε, SAM4 is better than SAM2 to reach a higher accuracy. In this numerical test, SAM4 produces accurate results and works well for 0 < ε ≤ 0.1.

Concerning TSF method, a uniform second order accurate in time for all 0 < ε ≤ 1 is obtained for both ρ ε and ρ ε v . The discretization error in τ -direction converges spectrally. In particular, the smaller the ε is, the less grid points are needed in τ . However, when ε = O(1), a substantial number of points in the τ direction is required to make the error in this direction small enough so as the uniform accuracy is recovered.

To compare the accuracy and efficiency of each numerical method in this test, we list of the errors in ρ ε together with computational time of all the numerical methods under some ∆t in Tables 1, 2 and 3, respectively representing the classical regime, the intermediate regime and the limit regime of ε.

From the comparisons in Tables 1-3, the following comments can be done. First, in the classical regime of ε, TSF, SIRK and ETD have similar accuracy, but TSF is much more expensive due to additional variable τ (which has to be sampled with N τ = 64). SIRK and ETD are then the most efficient method in this case. As ε diminishes (see Table 1), we see that the error of SIRK and ETD for a fixed ∆t increase. As expected, SAM methods do not work in ε = O(1) regime.

Second, in the intermediate regime of ε, the two SAM methods and TSF are more accurate than other methods. SAM4 is more accurate than SAM2 as ∆t is small, but SAM4 is twice more expensive than SAM2. When ∆t becomes small, TSF is more efficient than the SAM methods as we see in Table 2. This is due to the total computational cost of SAM methods which scales as O(n 2 0 ) where n 0 = T /∆t, due to micro step at each macro time level. For TSF, the total computational cost is O(n 0 N τ ) where N τ is fixed (small enough from the results in Figure 7) as ∆t reduces. Hence, the computation time of SAM increases quadratically as ∆t decreases, while that of TSF increases linearly. ETD shows second order accuracy before and after the regime ∆t ≈ 2πε 2 but with a significant change in error and computational time. Indeed, ETD has total computational cost O(n 2 0 ) similarly as SAM2 when ∆t > 2πε 2 , but for ∆t < 2πε 2 as the scheme degenerates to classical finite difference integrator the total cost is O(n 0 ). This comes from the fact that in (3.4), N vanishes when ∆t < 2πε 2 ; then, the error can become large since the size of the macro time step ∆t is similar to the period, which is not enough to get accurate results.

Lastly, in the limit regime of ε, ETD converges linearly and it is as costly as SAM2 for ∆t > 2πε 2 . SIRK could produce an accurate approximation until the error is stuck at O(ε). However, it is much more efficient than SAM methods. If one wants error further less than O(ε), TSF becomes the most efficient scheme since the workload due to the additional variable τ is much reduced.

Long-time behavior

In this part, we focus on more quantitative results and consider the long-time behavior of (1.1) with the initial condition (4.1). We look at the time evolution of the total energy defined as

H(t) := 1 2 R 2 ×R 2 |v| 2 f ε (t, x, v)dxdv + 1 2 R 2 |E ε (t, x)| 2 dx, t ≥ 0. (4.2)
The total energy H(t) of the Vlasov-Poisson equation (1.1) is conserved for all ε at the continuous level and is a good diagnostic to compare the numerical methods. We solve the Vlasov-Poisson equation (1.1) numerically till T = 32π with time step ∆t = π/2 7 (H = h = ∆t) and the same numerical parameters as previously and obtain the approximated energy H(t) at different time. We compute the relative energy error

|H(t) -H(0)|/|H(0)|, 0 < t = nπ 2 ≤ T,
for the different numerical methods under different regimes of ε. 1) and on the other side, the SIRK method (3.1) is not designed to capture the oscillations when ε is small. Thus, we show the time history of the relative energy error (in semi-log-scale) obtained with SIRK, ETD and TSF for ε = 1 and show the results obtained with ETD, SAM2, SAM4 and TSF for ε = 0.1, ε = 0.01 and ε = 0.001 in Figure 8.

In the case ε = 1, as SIRK is a L-stable method, it then enjoys a good long-time behavior. It can be observed than both ETD and TSF preserved in a good way the total energy. Note that for ε = 1 and ε = 0.1, ETD reduces to the classical central difference integrator (since we use ∆t = π/2 7 ). When ε becomes smaller (ε = 0.01 and ε = 0.001), ETD shows a slight increasing in the energy error.

Concerning the SAM methods, we observe that the error increases dramatically after some time. This is because the stroboscopic approximation (3.18) or (3.19) breaks (2.13). The perturbed problem (3.16) introduces non-pure imaginary eigenvalues λ, while the central difference integrator for u(t) = λu(t), t > 0, is not stable when λ is outside the imaginary axis. One solution to improve the long-time performance of SAM methods is to use schemes with larger stability region as the macro integrator for (3.16). For example, one can use the second order Runge-Kutta method (RK2) with time step small enough, but the two-stage method would obviously double the computational time of the whole SAM scheme. The energy error of the SAM2 scheme (3.18) with RK2 as the macro-solver (same ∆t = π/2 7 = H = O(h)) is shown in Figure 9 under different ε. From the results, we see the long-time behavior of the SAM method is significantly improved although there is still a small drift in the energy error. As being applied to study the Schrödinger equation in [START_REF] Ph | Solving highly-oscillatory NLS with SAM: Numerical efficiency and long-time behavior[END_REF], the macro solver with the best long-time behavior is suggested as the implicit midpoint rule which is A-stable and symplectic [START_REF] Hairer | Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF]. However, the implicit midpoint rule will make the integrator fully implicit, which would lead to the expensive inversion of a nonlinear system of size N p in our context. Finally, the TSF method (with N τ = 32) behaves well in the long-time run except for the case ε = 0.1. In the intermediate regime of ε, the TSF method has a significant increase in the energy which has been reported in [START_REF] Crouseilles | Uniformly accurate Particle-in-Cell method for the long time twodimensional Vlasov-Poisson equation with strong magnetic field[END_REF]. This is not due to the stability of the integrator and could be overcome by a refresh procedure. We can choose a t 0 > 0 to stop the algorithm and obtain the numerical solution of (2.2) at t 0 . Since (2.2) is autonomous, we then treat the approximated solution at t 0 as the initial data for (1.1) and restart the TSF method for solving (1.1) from t = 0. That is to say we reformulate (2.2) into (2.13), construct of the initial data (3.25) for the two-scale problem (3.22) and run (3.27) from t = 0 till t = t 0 to repeat the procedure. We plot the energy error without this refresh procedure in Figure 10. From the result in Figure 10, we observe that the energy error of TSF when ε = 0.1 starts to increase at around t = 6, so we choose the refresh time t 0 = 6. The relative energy error of the TSF method with a refresh every t 0 = 6 is shown in Figure 8. Clearly from the results, the TSF with refresh totally overcome the increase of the error and significantly improves its long-time behavior in the problematic intermediate regime.

At last, we compute the dynamics of the following initial data (which is a semi-Gaussian beam in v, as in [START_REF] Crouseilles | Uniformly accurate forward semi-Lagrangian methods for highly oscillatory Vlasov-Poisson equations[END_REF] 3) and we used the same numerical parameters as previously. We then solve the equation by the TSF method (with N τ = 32) for ε = 0.01 and look at the time history of ρ ε (t, x) = The results are shown in Figures 11 and12. First, we remark that the dynamics for ρ ε is similar to the one obtained with an isotropic initial distribution (see [START_REF] Frénod | Long time behavior of an exponential integrator for a Vlasov-Poisson system with strong magnetic field[END_REF][START_REF] Crouseilles | Uniformly accurate Particle-in-Cell method for the long time twodimensional Vlasov-Poisson equation with strong magnetic field[END_REF]). However, the dynamics of χ ε is strongly related to the velocity distribution of f 0 . In this case, the anisotropic character of f 0 induces some effects in the velocity direction which cannot be captured by the sole Guiding Center model. Note that in this regime, ETD and SAM methods produce very similar results.

Conclusion

We considered the long-time Vlasov-Poisson equation with a strong homogeneous magnetic field. The equation contains a small parameter ε which endows fast oscillations in time to the solution. We studied numerically several multiscale Particle-In-Cell (PIC) methods for solving this problem. The methods we considered include asymptotic preserving Runge-Kutta schemes, exponential time differencing schemes, stroboscopic averaging methods and uniformly accurate two-scale formulation method. Systematical comparisons in accuracy, efficiency and long-time behavior between the methods were performed for a wide range of 0 < ε ≤ 1 aiming at identifying the advantages and drawbacks of the different methods.

Figure 1 :

 1 Figure 1: Temporal error of SIRK-PIC at T = π/2 with respect to ∆t = H = O(h) and ε: maximum error in the density ρ ε .

Figure 2 :

 2 Figure 2: Convergence from the Vlasov-Poisson equation (1.1) to the limit model (1.3): maximum error with respect to at T = π/2 in the density ρ ε (left) and in the kinetic energy density ρ ε v (right).

3 ∆ t=T/2 4 ∆ t=T/2 5 ∆ t=T/2 6 ∆ t=T/ 2 7Figure 3 :

 345623 Figure 3: Temporal error of ETD-PIC at T = π/2 with respect to ∆t = H = O(h) and ε: maximum error in the density ρ ε (first row) and in the kinetic energy density ρ ε v (second row).

3 ∆ t=T/2 4 ∆ t=T/2 5 ∆ t=T/2 6 ∆

 3456 t=T/27

Figure 4 :

 4 Figure 4: Temporal error of SAM2-PIC at T = π/2 with respect to ∆t = H = O(h) and ε: maximum error in the density ρ ε (first row) and in the kinetic energy density ρ ε v (second row).

3 ∆ t=T/2 4 ∆ t=T/2 5 ∆ t=T/2 6 ∆ t=T/ 2 7Figure 5 :

 345625 Figure 5: Temporal error of SAM4-PIC at T = π/2 with respect to ∆t = H = O(h) and ε: maximum error in the density ρ ε (first row) and in the kinetic energy density ρ ε v (second row).

∆ t=T/2 3 ∆ t=T/2 4 ∆ t=T/2 5 ∆ t=T/2 6 ∆ 3 ∆ t=T/2 4 ∆ t=T/2 5 ∆ t=T/2 6 ∆

 34563456 t=T/27

Figure 6 :

 6 Figure 6: Temporal error of TSF-PIC at T = π/2 with respect to ∆t and ε: maximum error in the density ρ ε (first row) and in the kinetic energy density ρ ε v (second row).

Figure 7 :

 7 Figure 7: Discretization error of TSF-PIC in τ -direction with respect to Nτ at T = π/2: maximum error in the density ρ ε (left) and in the kinetic energy density ρ ε v (right).

  ) in the Vlasov-Poisson equation (1.1):f 0 (x, v) sin(x 2 ) + η cos(kx 1 )) e -(v 1 ) 2 2σ 2 1 [-r0,r0] (v 2 ),(4.3)where the function1 [-r0,r0] (v 2 ) = 1, if v 2 ∈ [-r 0 , r 0 ], otherwise 1 [-r0,r0] (v 2 ) = 0.We choose η = 0.05, kx 1 = 0.5, σ = 0.2, r 0 = 1.4 in (4.

Figure 10 :

 10 Figure 10: Relative energy error of TSF with refresh under ε = 0.1.

Figure 11 :

 11 Figure 11: Contour plot of quantity ρ ε (t, x) -1 in example (4.3) at different t with ε = 0.01.

Table 1 :

 1 Comparison in classical regime of ε: temporal errors in ρ ε and computational time (seconds) of different PIC methods for the Vlasov-Poisson (1.1) with ∆t = H = O(h) and Nτ = 64.

	ε = 1	∆t = T /2 3	∆t/2	∆t/4	∆t/8	∆t/16
	SAM2	2.72E-1	1.90E-1	1.94E-1	1.95E-1	1.95E-1
	time (cpu)	11.3	45.6	182	740	2930
	SAM4	2.79E-1	1.8E-1	1.84E-1	1.85E-1	1.85E-1
	time (cpu)	23.2	93.6	376	1470	5840
	TSF	2.20E-3	5.82E-4	1.48E-4	3.56E-5	8.00E-6
	time (cpu)	129	264	516	962	1918
	SIRK	1.10E-2	2.9E-3	7.57E-4	1.85E-4	4.74E-5
	time (cpu)	1.4	2.8	5.6	11.2	22.4
	ETD	1.70E-2	4.20E-3	1.10E-3	2.80E-4	7.50E-5
	time (cpu)	0.8	1.2	3.3	6.1	11.8
	ε = 0.5	∆t = T /2 3	∆t/2	∆t/4	∆t/8	∆t/16
	SAM2	9.43E-2	1.48E-2	1.44E-2	1.45E-2	1.45E-2
	time (cpu)	11.3	45.6	182	740	2930
	SAM4	9.61E-2	4.70E-3	2.70E-3	2.50E-3	2.50E-3
	time (cpu)	23.2	93.6	376	1470	5840
	TSF	2.80E-3	4.20E-4	5.60E-5	1.15E-5	2.28E-6
	time (cpu)	129	264	516	962	1918
	SIRK	2.66E-1	7.22E-2	1.72E-2	4.10E-3	1.00E-3
	time (cpu)	1.4	2.8	5.6	11.2	22.4
	ETD	1.23E-1	1.60E-2	3.60E-3	9.08E-4	2.29E-4
	time (cpu)	0.8	1.2	3.3	6.1	11.8

On the one side, the SAM methods (3.18) and

(3.19)

are not valid in classical regime ε = O(

Table 2 :

 2 Comparison in intermediate regime of ε: temporal errors in ρ ε and computational time (seconds) of different PIC methods for the Vlasov-Poisson (1.1) with ∆t = H = O(h) and Nτ = 16.

	ε = 0.1	∆t = T /2 3	∆t/2	∆t/4	∆t/8	∆t/16
	SAM2	5.60E-3	1.10E-3	2.26E-4	1.12E-4	1.14E-4
	time (cpu)	11.3	45.6	182	740	2930
	SAM4	5.70E-3	1.10E-3	2.76E-4	6.50E-5	1.90E-5
	time (cpu)	23.2	93.6	376	1470	5840
	TSF	1.20E-2	3.40E-3	9.77E-4	3.12E-4	5.73E-5
	time (cpu)	39.8	72.4	142	266	532
	SIRK	1.97E-1	1.95E-1	1.94E-1	1.92E-1	2.15E-1
	time (cpu)	1.4	2.8	5.6	11.2	22.4
	ETD	8.58E-2	2.29E-2	1.48	9.37E-1	1.85E-1
	time (cpu)	12	47.6	3.3	6.1	11.7
	ε = 0.05	∆t = T /2 3	∆t/2	∆t/4	∆t/8	∆t/16
	SAM2	5.40E-3	1.30E-3	3.20E-4	8.63E-5	3.34E-5
	time (cpu)	11.3	45.6	182	740	2930
	SAM4	5.40E-3	1.40E-3	3.27E-4	9.01E-5	3.46E-5
	time (cpu)	23.2	93.6	376	1470	5840
	TSF	1.18E-2	3.00E-3	7.55E-4	2.30E-4	7.31E-5
	time (cpu)	39.8	72.4	142	266	532
	SIRK	9.83E-2	1.00E-1	1.02E-1	1.02E-1	1.02E-1
	time (cpu)	1.4	2.8	5.6	11.2	22.4
	ETD	1.07E-1	2.79E-2	8.70E-3	1.70E-3	1.09
	time (cpu)	12	47.6	188	727	11.8

Table 3 :

 3 Comparison in limit regime of ε: temporal errors in ρ ε and computational time (seconds) of different PIC methods for the Vlasov-Poisson (1.1) with ∆t = H = O(h) and Nτ = 4.

	ε = 0.001	∆t = T /2 3	∆t/2	∆t/4	∆t/8	∆t/16
	SAM2	6.00E-3	1.50E-3	3.77E-4	9.16E-5	2.27E-5
	time (cpu)	11.3	45.6	182	740	2930
	SAM4	6.00E-3	1.50E-3	3.77E-4	9.16E-5	2.27E-5
	time (cpu)	23.2	93.6	376	1470	5840
	TSF	1.15E-2	3.00E-3	7.60E-4	1.82E-4	3.66E-5
	time (cpu)	12.4	22.6	41.7	81.6	162
	SIRK	2.60E-3	2.00E-3	1.90E-3	1.90E-3	1.90E-3
	time (cpu)	1.4	2.8	5.6	11.2	22.4
	ETD	1.33E-1	3.29E-2	1.59E-2	8.00E-3	4.00E-3
	time (cpu)	13.7	48.5	193	738	2910
	ε = 0.0001	∆t = T /2 3	∆t/2	∆t/4	∆t/8	∆t/16
	SAM2	5.00E-3	1.50E-3	3.76E-4	9.14E-5	2.25E-5
	time (cpu)	11.3	45.6	182	740	2930
	SAM4	6.00E-3	1.50E-3	3.76E-4	9.14E-5	2.25E-5
	time (cpu)	23.2	93.6	376	1470	5840
	TSF	1.16E-2	3.00E-3	7.61E-4	1.82E-4	3.66E-5
	time (cpu)	12.4	22.6	41.7	81.6	162
	SIRK	2.00E-3	5.18E-4	2.35E-4	2.08E-4	2.01E-4
	time (cpu)	1.4	2.8	5.6	11.2	22.4
	ETD	1.33E-1	3.18E-2	1.58E-2	8.30E-3	4.10E-3
	time (cpu)	13.7	48.5	193	738	2910

R 2 f ε (t, x, v)dv and χ ε (t, v) := R 2 f ε (t, x, v)dx,

Acknowledgements

This work was partly supported by the French ANR project MOONRISE ANR-14-CE23-0007-01. N. Crouseilles and S. Hirstoaga are supported by the Enabling Research EUROFusion project CfP-WP14-ER-01/IPP-03. X. Zhao is supported by the IPL FRATRES.


Van Long Tran
email: vanlong.tran@telecom-sudparis.eu

Eric Renault
email: eric.renault@telecom-sudparis.eu

Viet Hai Ha

Analysis and evaluation of the performance of CAPE

Keywords: CAPE, OpenMP, MPI, high performance computing, parallel programming

MPI (Message Passing Interface) and OpenMP are two tools broadly used to develop parallel programs. On the one hand, MPI has the advantage of high performance while being difficult to use. On the other hand, OpenMP is very easy to use but is restricted to shared-memory architectures. CAPE is an approach based on checkpoints to allow the execution of OpenMP programs on distributed-memory architectures. This paper aims at presenting both an in-depth analysis and an evaluation of the performance of CAPE by comparing the execution model of both CAPE and MPI. Some suggestions are also provided to improve the use of CAPE.

I. INTRODUCTION

In order to minimize programmers' difficulties when developing parallel applications, a parallel programming tool at a higher level should be as easy-to-use as possible. MPI [START_REF]Message passing interface forum[END_REF] (which stands for Message Passing Interface) and OpenMP [START_REF] Openmp | Openmp application program interface version 4[END_REF] are two widely-used tools that meet this requirement. MPI is a tool for high-performance computing on distributedmemory environments, while OpenMP has been developed for shared-memory architectures. If MPI is quite difficult to use, especially for non programmers, OpenMP is very easy to use, requesting the programmer to tag the pieces of the code to execute in parallel.

Some efforts have been made to port OpenMP on distributed-memory architectures. However, no solution has successfully met both requirements: 1) to be fully compliant with the OpenMP standard and 2) high performance. Most prominent approaches include the use of an SSI [START_REF] Morin | Kerrighed: a single system image cluster operating system for high performance computing[END_REF], SCASH [START_REF] Sato | Cluster-enabled openmp: An openmp compiler for the scash software distributed shared memory system[END_REF], the use of the RC model [START_REF] Karlsson | A fully compliant openmp implementation on software distributed shared memory[END_REF], performing a source-to-source translation to a tool like MPI [START_REF] Basumallik | Towards automatic translation of openmp to mpi[END_REF] [START_REF] Dorta | Implementing openmp for clusters on top of mpi[END_REF] or Global Array [START_REF] Huang | Towards a more efficient implementation of openmp for clusters via translation to global arrays[END_REF], or Cluster OpenMP [START_REF] Hoeflinger | Extending openmp to clusters[END_REF]. Among all these solutions, using SSI (for Single System Image) is the most straightforward approach. An SSI includes a DSM (for Distributed Shared Memory) to provide an abstracted shared-memory view over a physical distributedmemory architecture. The main advantage of this approach is its ability to easily provide a fully-compliant version of OpenMP. Thanks to their shared-memory nature, OpenMP programs can easily be compiled and run as processes on different computers in the SSI. However, as the shared memory is accessed through the network, the synchronization between the memories involves an important delay access time which makes this approach hardly scalable. Some experiments [START_REF] Morin | Kerrighed: a single system image cluster operating system for high performance computing[END_REF] show that the larger the number of threads, the lower the performance. As a result, in order to reduce the execution time overhead involved by the use of an SSI, other approaches have been proposed, like SCASH that maps only the shared variables of the processes onto a shared-memory area attached to each process, the other variables being stored in a private memory, and the RC model that uses the relaxed consistency memory model. However, these approaches have difficulties to identify the shared variables automatically. As a result, no fully-compliant implementation of OpenMP based on these approaches has been released so far. Some other approaches aim at performing a source-to-source translation of the OpenMP code to an MPI code. This approach allows the generation of high-performance codes on distributed-memory architectures. However, not all OpenMP directives and constructs can be implemented. Cluster OpenMP, proposed by Intel, also requires the use of additional indicators of its own (ie. not included in the OpenMP standard). As a result, this one cannot be considered as a fully-compliant implementation of the OpenMP standard.

To execute programs including OpenMP directives on distributed-memory architectures, CAPE (which stands for Checkpointing-Aided Parallel Execution) adopts a completely different approach based on the use of checkpointing techniques. This allows to distribute the work of the parallel block programs to the different processes of the system, and to handle the exchange of shared data automatically. Two versions of CAPE [START_REF] Renault | Distributed implementation of openmp based on checkpointing aided parallel execution[END_REF] [START_REF] Ha | Improving performance of cape using discontinuous incremental checkpointing[END_REF] have been released that proves the effectiveness of the approach and the high-performance.

This paper focuses on the in-depth analysis of the CAPE architecture and an evaluation of its performance. A comparison with MPI is proposed and some remarks are also provided for a better use of CAPE.

The next section presents an overview of some related work and serves as an introduction to MPI, OpenMP and checkpointing techniques. Section III introduces CAPE using the discontinuous incremental checkpointer. Section IV compares the operational model of both CAPE and MPI. Section V presents some experimental results to support the analysis in Sec. IV. The last section outlines some conclusions and directions for future research.

II. RELATED WORK

A. MPI

MPI [START_REF]Message passing interface forum[END_REF] is an Application Programming Interface (API) developed during the 90s. This interface provides essential point-to-point communications, collective operations, synchronization, virtual topologies, and other communication facilities for a set of processes in a language-independent way, with a language-specific syntax, plus a small set of language-specific features. MPI uses the SPMD (Single Program Multiple Data) programming model where a single program is run independently on each node. A starter program is used to launch all processes, with one or more processes on each node and all processes first synchronize before executing any instructions. At the end of the execution, before exiting, processes all synchronize again. Any exchange of information can occur in between these two synchronizations. A typical example is the master-slave paradigm where the master distributes the computations among the slaves and each slave returns its result to master after the job complete.

Although it is capable of providing high performance and running on both shared-and distributed-memory architectures, MPI programs require programmers to explicitly divide the program into blocks, one for the master process and the others for slave processes. Moreover, some tasks, like sending and receiving data or the synchronization of processes, must be explicitly specified in the program. This makes it difficult for programmers for two main reasons. First, it requires the programmer to organize the program into parallel structures which are bit more complex structures than in sequential programs. This is even more difficult for the parallelization of a sequential program, as it severely disrupt the original structure of the program. Second, some MPI functions are difficult to understand and use. For these reasons, and despite the fact that MPI achieves high performance and provides a good abstraction of communications, MPI is still considered as assembly for parallel programming.

B. OpenMP

OpenMP [START_REF] Openmp | Openmp application program interface version 4[END_REF] provides a higher level of abstraction than MPI for the development of parallel programs. It consists in a set of directives, functions and environment variables to easily support the transformation of a C, C++ or Fortran sequential program into a parallel program. A programmer can start writing a program on the basis of one of the supported languages, compile it, test it, and then gradually introduce parallelism by the mean of OpenMP directives. For C, C++ or Fortran compilers supporting OpenMP, directives are taken into consideration through the execution of specific codes. OpenMP offers directives for the parallelization of most classical parallel operations including loops, parallel sections, etc.

OpenMP follows the fork-join execution model with threads. When the program starts, it runs only one thread (the master thread) and executes the code sequentially. When it reaches an OpenMP parallel directive, the master thread spawns a team work including itself and a set of secondary threads (the fork phase), and the work allocated to each thread in this team starts. After secondary threads complete their allocated work, results are updated in the memory space of the master thread and all secondary threads are suspended (the join phase). Thus, after the join phase, the program executes only one thread as per the original. Fork-join phases can be carried out several times in a program and can be nested.

OpenMP is based on the use of threads. Thus, it can only run on shared-memory architectures as all threads use the same memory area. However, in order to speed up the computations, OpenMP uses the relaxed-consistency memory model. With this memory model, threads can use a local memory to improve memory accesses. The synchronization of the memory space of the threads is performed automatically both at the beginning and at the end of each parallel construct, or can be performed manually by specifying flush directives.

C. Checkpointing

Checkpointing aims at saving the state of a running program so that its execution can be resumed from that point later in time [START_REF] Plank | An overview of checkpointing in uniprocessor and distributed systems, focusing on implementation and performance[END_REF]. Checkpoints are used in many applications, for example to enable fault-tolerance or for backing up systems.

Checkpointing techniques can be divided into two groups. Complete checkpointing is a technique for which all the memory space and the associated information of a process is saved in each checkpoint. Its major drawbacks are the redundancy of data when several checkpoints are taken consecutively without many changes in the process memory space, and the large size of the checkpoints. Incremental checkpointing aims at saving only the part of the memory space and the associated information that have been updated since the initialization of the program or the last checkpoint. In order to identify the updated memory area, a checkpointed program is usually monitored by another one that performs the saving regularly.

CAPE requires taking checkpoints of each part of the continuing process. As a result, in order to serve CAPE in an optimal way, we developed a checkpointing technique called DICKPT [START_REF] Ha | Discontinuous incremental: A new approach towards extremely lightweight checkpoints[END_REF] (which stands for Discontinuous Incremental Checkpointing). DICKPT is based on incremental checkpoints with additional image capabilities for discrete segments of the process. The analysis of experimental data has demonstrated the superior performance of this new capability both in terms of the reduction of the size of checkpoints and the reduction of the time to generate them.

III. CAPE BASE ON DISCONTINUOUS INCREMENTAL CHECKPOINTING

A. Execution model CAPE is the alternative approach to allow the execution of OpenMP programs on distributed-memory systems. CAPE is based on a process as a parallel unit, which is different from the traditional implementations of OpenMP where the parallel unit is a thread. All the important tasks of the forkjoin model are automatically generated by CAPE based on checkpointing techniques, such as task division, reception of results, updating results into the main process, etc. In its first version, CAPE used complete checkpoints so as to prove the concept. However, as the size of complete checkpoints is very large, it takes a lot of traffic on the network to transfer data between processes and involves a high cost for the comparison of the data from the different complete checkpoints to extract the modifications. These factors have significantly reduced the performance and the scalability of our solution. Fortunately, these drawbacks have been overcome in the second version of CAPE based on DICKPT.

Figure 1 describes the execution model of the second version of CAPE using three nodes. At the beginning, the program is initialized on all nodes and the same sequential code block is executed on all nodes. When reaching an OpenMP parallel structure, the master process divides the tasks into several parts and sends them to slave processes using DICKPT. Note that these checkpoints are very small in size, typically very few bytes, as they only contain the results of some very simple instructions to make the difference between the threads, which do not change the memory space that much. At each slave node, after receiving a checkpoint, it is injected into the local memory space and initialized for resuming. Then, the slave process executes the assigned task, extracts the result, and creates a resulting checkpoint. This last checkpoint is sent back to the master process. The master process then combines them altogether, injects the result into its memory space, sends it to all the other slave processes to synchronize the memory space of all processes and prepares for the execution of the next instruction of the program. In order to translate OpenMP programs, CAPE uses a set of templates. After the source-to-source translation of the OpenMP program, the code is regularly compiled by a C/C++ compiler. Figure 2 shows the different steps of the CAPE compilation process for a C/C++ OpenMP program. The current version of CAPE assumes that the parallel statements satisfy the Bernstein's conditions, which mean that all components of the code can be executed independently, no input and output data are common, and there is no communication between nodes during the execution.

Among all OpenMP parallel constructs, parallel for is probably the most important one, as 1) this is the most used construct in user programs and 2) it can act as an intermediate structure to convert other OpenMP parallel constructs including for example the parallel sections one. As a result, instead of building a complete template for each construct, some are first translated into the parallel for construct. CAPE automatically translates the parallel for construct to a set of function calls that contains the fundamental operations of CAPE. The template for the parallel for construct is shown in Figure 3. In order to make the presentation easier to understand, the template is presented in a version where the number of iterations of the for statement is equal to the number of threads. Below is the list of basic functions used by CAPE.

1) start ( ) set the environment to create DICKPT. 2) stop ( ) suspends the environment to create checkpoints. 3) create ( file ) create a checkpoint and save it to file. 4) inject ( file ) updates the current process with the information contained in checkpoint file provided as a parameter. 5) send ( file, node ) sends the content of file to node. 6) wait_for ( file ) waits and merges all the components of the file . 7) last_parallel ( ) returns TRUE when the current parallel block is the last one of the entire program and FALSE otherwise. 8) merge ( f ile 1 , f ile 2 ) applies the list of modifications from f ile 2 to checkpoint f ile 1 . 9) broadcast ( file ) sends file to all slave nodes. 10) receive ( file ) waits for file to be available.

The main steps of CAPE can be explained as follows: On the master side (lines 1 to 12), the first statement of the for loop is executed. At each for loop iteration, the master distributes jobs to each slave by generating and sending it as an incremental checkpoint. Then, it waits for the results from the slave process. After the result is received, it updates the memory. Finally, the master node broadcasts the result to all slave nodes to guarantee the memory consistency if necessary.

On the slave process side (lines 13 to 25), each process receives a DICKPT from the master. Then, the checkpoint is injected into the memory process space which executes the assigned tasks. At last, the result of the execution on the slave node is extracted into an incremental checkpoint and sent back to the master node to update its memory. IV. ANALYSIS AND PERFORMANCE EVALUATION A. MPI and CAPE are almost the same... Similarities between CAPE and MPI are many including the architectural background, organizational models, and execution models. From the architectural point of view, MPI and CAPE both have the ability to execute on multi-processor architectures, and the most important goal for both is to support systems based on a distributed memory such as clusters, grids, and clouds. From the organizational-and programmingmodel point of view, both use processes and the masterslave paradigm, with a single process attached to a node. Execution models for CAPE and MPI contain similar steps. To execute a program that contains parallel and sequential code alternatively, the steps are i) the initialization of the program, and then ii) the execution of both sequential and parallel sections on a set of processes. The initialization is the same for CAPE and MPI. All nodes are initialized at the beginning of the execution and one of the nodes is chosen as the master, the others being considered as slaves.

#

B. ... but not exactly

The main difference between CAPE and MPI is the determination of the data to exchange and the way data in the memory space of a process are updated. MPI was defined in the scope of a third generation language. The user is required to explicitly define what (the address, the length and the data type) and when (send and receive) to transfer the memory areas. In order to do so, MPI is based on a set of functions which requires the user to specify the operations that have to be executed step by step. However, CAPE is based on OpenMP which was defined as a fourth generation language. In this case, the user specifies the expected result and it is the compiler work to automatically translate this objective into effective code. For the specific case of CAPE, everything is automatically done through the use of checkpoints. The result of the execution of a code block is saved in a discontinuous incremental checkpoint initialized before executing the block and saved at the end of the block. The checkpoint is injected into the memory space of the target process using the inject function (as shown above).

The ability for CAPE to automatically distribute the work, to extract and inject the results into the memory space of the process, creates an advantage in the transparency of the communication of data between processes for the programmer. This definitively helps the implementation of a fullycompatible version of OpenMP for distributed-memory architectures. At the same time, the use of checkpoints is also the main drawback of CAPE as it slows down the execution.

C. Qualitative comparison

Table I presents a qualitative comparison of the different steps when executing a parallel code with both MPI and CAPE. As mentioned above, steps for initialization and execution are similar for both CAPE and MPI. As a result, they are not mentioned on the table.

From Table I, let t b be the time for work division and distribution to the slaves, t c be the computation time, t u be the time for result extraction and returning the data to the master, and t s be the time for the memory update on all processes. The master automatically divides the task using DICKPT and sends them to the slaves.

Transparent to the programmer. Not much time wasted on running.

The programmer divides the code into tasks for all processes.

Not straightforward.

Computation. (tc)

Slaves execute the application code under the supervision of DICKPT.

Overhead due to the monitoring of the process by the checkpointer.

Processes execute the application code independently.

High performance

Result extraction and return to the master. (tu)

Results

automatically extracted and sent back using DICKPT.

Transparent for the programmer. Accurate extraction of results. Overhead spent on filtering results.

The programmer writes the code to determine data to send and where to store the results.

High performance. Difficult for the programmer.

Memory update on all processes. (ts)

Memory space of all processes automatically updated by DICKPT.

Transparent for the programmer. Automatic synchronization. Faster than MPI.

The programmer writes the synchronization code.

Difficult for the programmer.

For any given program, let t tt and t ss be the execution time of the sequential part and the parallel part of the code respectively. Let t be the total amount of time to execute the program. It is given by:

t = t tt + t ss (1)
= t tt + (t b + t c + t u + t s ) (2)
Note that executing a program in parallel is meaningful when the execution time of the parallel part of a program is far greater than the execution time for the initialization and the communication between the nodes. It is the same for the time spent in the synchronization of the memory space of the processes. As a result, in the following, the time for initialization and communication as negligible as compared to the time for computation and synchronization, ie. t b = t s = 0. This leads to:

t = t tt + (t c + t u ) (3)
Then, it is fair to consider that the sequential code executed on all the nodes with both MPI and CAPE is the same and thus that the time is equal for both techniques. The time difference between CAPE and MPI becomes:

∆t = ∆(t c + t u ) (4)
Basically, the computation of the parallel code should be the same for MPI and CAPE, and so the time difference in t c as the cost of monitoring and generating DICKPT, in the application process is always monitored by a checkpointer process. Therefore, the execution time of CAPE applications is higher than the one of MPI. According to [START_REF] Gioiosa | Transparent, incremental checkpointing at kernel level: a foundation for fault tolerance for parallel computers[END_REF], the time difference is in the 2% to 20%. The time difference in t u caused by the process update via DICKPT is more complex. However, for the problems have the high t c , so the t u is too small in the total time of the program. Therefore we can ignore that time.

V. EXPERIMENTAL RESULTS

In order to validate our approach, some performance measurements have been conducted on a Desktop Cluster. This testbed is composed of nodes including Intel (R) Core (TM) 2 Duo E8400 CPUs running at 3 GHz and 2 GB RAM, operated by Linux kernel 2.6.35 with the Ubuntu 10.10 flavor, and connected by a standard Ethernet at 100 Mb/s. In order to avoid as much as possible external influences, the entire system was dedicated to the tests during performance measurements.

The program used for tests is a matrix-matrix product for which the size varies from 3,000×3,000 to 12,000×12,000. Matrices are supposed to be dense and no specific algorithm has been implemented to take into account sparse matrices. Each experiment has been performed at least 10 times and a confidence interval of at least 90% has always been achieved for the measures. Data reported here are the means of the 10 measures.

The execution of both the sequential version and the OpenMP version of the program on one of the nodes gives the result provided in Table II. A single core was used for the sequential execution of the program, while the OpenMP program takes benefits of the two cores. One can check that results in Table II are consistent as the execution time for both sequential and OpenMP versions are directly proportional to the cube of the matrix size. Typically, this means that no important cache effects have polluted the performance measurements, probably because almost all data fit into memory. Moreover, the speedup obtained by OpenMP is 1.8 for the first three matrix sizes and 1.65 for the last one, which are expected values. The next part of the experiment consists in comparing CAPE and MPI. Figure 4 and 5 present the execution time in seconds of the matrix-matrix program for various numbers of nodes and matrix sizes. Note that, despite the fact that processors are dual core, a single core was used during the experiments. Two measures are represented each time: the left one is associated with CAPE using incremental checkpoints, and the right one is associated with MPI. Both Figure 4 and 5 contain two groups of graphs. The group above belongs to the master node, while the one below belongs to the slaves. The line numbers for the code refer to Figure 3. Each group shows four graphs and the meaning is as follows:

• Init is the elapsed time between the beginning of the program and the beginning of the parallel for loop in the matrix-matrix product. On Fig. 3, these are all lines before the first one. In the modeling above, this is time t tt . • Before is the time spent to create and send checkpoints (lines 2 to 5) on the master node. On slave nodes, this includes waiting for and receiving the checkpoint, and then updating the slave process using the checkpoint (lines 16 and 17). For the case of MPI, this is the time to send data to slave nodes. In the modeling above, this is time t b . • Compute is the time to generate the last checkpoint on the master node (lines 7 and 8) and the time to do the job (execute the code of matrix-matrix product) on the slaves (lines 18 and 19). In the modeling above, this is time t c . • Update is the time spend in waiting for and receiving all updates from the slave nodes and inject them into the master process memory (lines 9 and 10). On slave nodes, this is the time to generate the incremental checkpoints and send them to the master node (lines 20 to 22). For the case of MPI, this is the time to send results from slave nodes to the master node. In the modelling above, this is time t u . Figure 4 presents the execution time in the detailed phases with different numbers of nodes. The size of matrices are 9,000×9,000. However, similar trends are observed for the other matrix sizes. In the master node, the execution times for the Init, Before and Compute steps are small as compared to step Update (t tt , t b and t c are very small as compared to t u ). This is because step Update waits for the results computed in the slaves. This means that t u includes the computation time of the slaves to do the assigned tasks. Then, one can see that t is similar to t u on the master node. Similarly to the slaves, the execution times of both Before and Init steps are too small and can be ignored. As a result, the total execution time consists in both the computation time and the update time. This clearly conforms to equations 3 and 4.

Graphs in Figure 4 also show that, for the case with three nodes, the execution time for CAPE is always larger than the one for MPI. This is due to the fact that slaves compute and update the result of half of the matrix, and the size of the checkpoint takes a large part of the total memory space of the process. Therefore, the impact of the checkpointer on the total execution time is very important. However, when increasing the number of machines, the size the checkpoints decreases together with ∆t. When increasing the number of nodes to 31, this time is approximately the same for both CAPE and MPI. This shows that the performance of CAPE is poor in systems with a very small number of nodes. However, no one would use a distributed-memory machine to run a program requiring a small number of threads. Any sharedmemory machine would be sufficient in that case. On the other hand, the same graphs also show that CAPE has the ability to achieve high performance when the size of checkpoints is small. This is involved by the fact that slaves only update a small amount of memory space.

Figure 5 presents the execution time for different matrix sizes for a matrix-matrix product. Thirty one computers have been used to perform the parallel computations. However, the results tend to be similar to the experimental system when the numbers of nodes is different (except for the case with 3 nodes). Comments to chart in Figure 4 are also true for this one. The data on the graph shows that the execution time for both MPI and CAPE are proportional to the size of the matrix. As both MPI and CAPE are only using the data transfer mechanism to update the computation results, the execution speed is also quite similar. In a more detailed analysis, data show that the total execution time for CAPE based on DICKPT is only about 10% higher than the execution time for MPI, except for matrices with a size of 3000x3000 which rate is 1.3. Figure 6 shows the speedup of CAPE using incremental checkpoints for various number of nodes and matrix sizes. The dotted line represents the maximum theoretical speedup. The figure clearly shows that the solution provides an efficiency (the ratio of the speedup over the number of nodes) in the range from 75% to 90% which is very good. Also, it highlights that the larger the size of matrices, the higher the speedup.

VI. CONCLUSION AND FUTURE WORK

CAPE, with its ability to automatically divide the task and establish data communication between processes, has demonstrated its ability to install a conforming implementation of OpenMP on distributed-memory systems. Moreover, the comparison of the performance results of MPI and CAPE by means of both a theoretical and an experimental analysis has proven that the efficiency of CAPE programs is very close to the one of MPI programs, and that the difference between both is getting smaller and smaller as the size of the problem and the number of nodes in the system increases.

For the near future, we will keep on developing CAPE to support the OpenMP constructs that have not been ported yet. Besides, we would like to study other models to take benefits of CAPE so that our model can improve the flexibility of the execution of many more sequential programs.

Fig. 1 .

 1 Fig. 1. The execution model of CAPE.

Fig. 2 .

 2 Fig. 2. Translation OpenMP programs with CAPE

Fig. 3 .

 3 Fig. 3. Template for the parallel for with incremental checkpoints.

Fig. 4 .

 4 Fig. 4. Execution time (in seconds) vs. number of nodes.

Fig. 5 .

 5 Fig. 5. Execution time (in seconds) vs. problem size.

Fig. 6 .

 6 Fig. 6. Speedup vs. number of nodes.

TABLE I THE

 I STEPS THAT EXECUTE A PARALLEL CODE BLOCK OF CAPE AND MPI

	Steps	Approach	CAPE Pros and cons	Approach	MPI Pros and cons
	Work	division
	and distribution
	to slaves.
	(t b )

TABLE II EXECUTION

 II TIME (IN SECONDS) ON A SINGLE NODE.

	Size	Sequential	OpenMP
	3000 x 3000	258.9	142.4
	6000 x 6000	1852.7	1048.7
	9000 x 9000	7314.5	3986.2
	12000 x 12000	14990.5	8999.4