Hubert Baty

Florence Drui

Philippe Helluy

Emmanuel Franck

Christian Klingenberg

Lukas Thanhäuser

A robust and efficient solver based on kinetic schemes for Magnetohydrodynamics (MHD) equations

A robust and efficient solver based on kinetic schemes for Magnetohydrodynamics (MHD) equations.

Introduction

The MagnetoHydroDynamics (MHD) system is a fundamental model used in many fields of physics: astrophysics, plasma physics, geophysics... Indeed, the MHD model is commonly adopted as an excellent framework for collisional plasma environments. The numerical approximation of this system is a difficult task. Compared to other models in fluid mechanics, it contains more conservative unknowns and thus more scales and more different wave speeds. It is also subject to complex phenomena such as occurrence of shock waves, current sheet formation, magnetic reconnection, instabilities and turbulent behaviors. An additional specificity is that the magnetic field has to satisfy a divergence-free condition, which is generally difficult to be verified by numerical solutions. In order to deal with the divergence-free condition, we adopt here a modified version of the MHD equations by a divergence cleaning term proposed in [START_REF] Dedner | Hyperbolic divergence cleaning for the MHD equations[END_REF].

We propose a simple scheme, based on an abstract kinetic interpretation, for computing two-dimensional MHD solutions. The kinetic interpretation is vectorial and has been first proposed by Bouchut in [START_REF] Bouchut | Construction of BGK models with a family of kinetic entropies for a given system of conservation laws[END_REF] and Aregba-Natalini in [START_REF] Aregba-Driollet | Discrete kinetic schemes for multidimensional systems of conservation laws[END_REF]. In this approach, the original system of conservation laws is represented by a system of transport equations coupled through a so-called collision source term. The idea originates from the Boltzmann kinetic theory of gases, but it is purely abstract and has no physical meaning. It leads to natural numerical methods, where the transport step and the collision step are made separately. Other interesting features arise from this representation (see [START_REF] Badwaik | Task-based parallelization of an implicit kinetic scheme[END_REF][START_REF] Coulette | High-order implicit palindromic Discontinuous Galerkin method for kinetic-relaxation approximation[END_REF], for instance).

In this paper, we solve the kinetic representation with a Lattice-Boltzmann Method (LBM), where the transport step is solved exactly on a regular grid. The same approach has already been used in [START_REF] Dubois | Simulation of strong nonlinear waves with vectorial lattice boltzmann schemes[END_REF] for solving compressible Euler equations. An important aspect of the LBM is the choice of the relaxation parameter in the collision step. The choice of the parameter allows adjusting the numerical viscosity of the LBM scheme. We provide an analysis in a simplified one-dimensional framework which shows that it is possible to adjust more precisely the numerical viscosity with a generalized matrix relaxation parameter. We also show that the divergence cleaning effect is improved if the relaxation parameter is chosen differently for the physical variables and the divergence cleaning potential.

In order to capture fine structures, it is necessary to consider very fine meshes. We have programmed the algorithm in a very efficient way in order to address recent GPUs (Graphic Processing Units) or multicore CPUs. We describe the implementation, which relies on OpenCL and PyOpenCL, and the memory optimizations used for reaching high performance. The program allows performing full MHD simulations on grids as fine as 7000×7000 within a few hours.

Finally, we apply the whole approach to several MHD simulations: the classical Orszag-Tang test and a more physical study of ideal MHD instabilities (tilt modes) and associated formation of quasi-singular current sheets.

Mathematical model

The MagnetoHydroDynamic (MHD) system is a model used in many fields of physics. It consists of an extension of the compressible Euler equations for taking into account magnetic effects. A difficulty is that the magnetic field has to satisfy a divergence-free condition. This condition expresses that magnetic charges are not observed in Nature. Standard finite volume methods do not guarantee that the numerical magnetic field is divergence-free. More annoying: the divergence errors generally grow with the simulation time, which leads to physically wrong results. For limiting the divergence errors, we adopt the divergence cleaning method described in [START_REF] Dedner | Hyperbolic divergence cleaning for the MHD equations[END_REF].

MHD equations with divergence cleaning

We consider the MHD equations with Divergence Cleaning [START_REF] Dedner | Hyperbolic divergence cleaning for the MHD equations[END_REF] (called the MHD-DC equations in the following)

∂ t       ρ ρu Q B ψ       + ∇ •       ρu ρu ⊗ u + (p + B•B 2 )I -B ⊗ B (Q + p + B•B 2 )u -(B • u)B u ⊗ B -B ⊗ u + ψI c 2 h B       =       0 0 0 0 0      
.

The velocity and magnetic field are denoted

u = (u 1 , u 2 , u 3 ) T , B = (B 1 , B 2 , B 3 ) T ,
the pressure is given by a perfect-gas law with a constant polytropic exponent γ > 1

p = (γ -1)(Q -ρ u • u 2 - B • B 2 ).
The other variables are the density ρ, the total energy Q, the divergence cleaning potential ψ. The divergence cleaning velocity is a positive parameter c h > 0.

When the magnetic field satisfies the divergence-free condition

∇ • B = 0,
and the potential ψ is a constant, then the MHD-DC equations simply reduce to the usual MHD equations. In other words, the MHD-DC equations are a generalization of the MHD equations, where the magnetic field can have a non-vanishing divergence. It is possible to show that the MHD-DC system, as the MHD system, is hyperbolic [START_REF] Dedner | Hyperbolic divergence cleaning for the MHD equations[END_REF]. However this system is not strictly hyperbolic, which leads to some difficulties, such as non-uniqueness of the Riemann problems in some situations [START_REF] Torrilhon | Non-uniform convergence of finite volume schemes for Riemann problems of ideal magnetohydrodynamics[END_REF].

The interest of the MHD-DC formulation is for numerical approximations. Indeed, standard approximations of the usual MHD equations suffer from drifting errors along time of the divergence constraint. This has been observed by several authors [START_REF] Powell | An approximate Riemann solver for magnetohydrodynamics (that works in more than one space dimension[END_REF][START_REF] Dedner | Hyperbolic divergence cleaning for the MHD equations[END_REF][START_REF] Barth | On the role of involutions in the Discontinuous Galerkin discretization of Maxwell and magnetohydrodynamic systems[END_REF]. A review on this topic is developed for instance in [START_REF] Tóth | The ∇ • B = 0 constraint in shock-capturing magnetohydrodynamics codes[END_REF]. Approximations of the MHD-DC generally have a much better behavior. In this generalized model, the divergence errors propagate at the wave speed c h . The errors are then damped at the boundaries of the computational domain or by the numerical diffusion in the domain. In [START_REF] Dedner | Hyperbolic divergence cleaning for the MHD equations[END_REF] other divergence corrections are proposed. We choose a correction that leads to a conservative and hyperbolic model, in order to be able to apply the general theory of vectorial kinetic representations given in [START_REF] Bouchut | Construction of BGK models with a family of kinetic entropies for a given system of conservation laws[END_REF][START_REF] Aregba-Driollet | Discrete kinetic schemes for multidimensional systems of conservation laws[END_REF]. Theoretically, the parameter c h can take any value. But in practice it is generally chosen larger than all the wave speeds of the MHD system (see [START_REF] Dedner | Hyperbolic divergence cleaning for the MHD equations[END_REF]).

We introduce the conservative variables

w = w(x, t) =       ρ ρu Q B ψ       ∈ R m , m = 9.
and the flux (n is a vector of R 3 )

F(w, n) =       ρu • n ρu • nu + (p + B•B 2 )n -B • nB (Q + p + B•B 2 )u • n -(B • u)B • n u • nB -B • nu + ψn c 2 h B • n       ∈ R m .
In this work, we assume that all the fields do not depend on the x 3 space variable. We are thus computing two-dimensional solutions. If we set

n 1 = (1, 0, 0) T , n 2 = (0, 1, 0) T , F 1 = F(w, n 1 ), F 2 = F(w, n 2 ),
the MHD equations can also be written as a two-dimensional system of nine conservation laws

∂ t w + ∂ 1 F 1 + ∂ 2 F 2 = 0, (1)
where we use the notation ∂ i for the partial derivative ∂/∂x i . If the fields depend only on the x 1 space variable, the system reduces to

∂ t w + ∂ 1 F 1 = 0. (2)
In this case, the mathematical analysis is simplified. We shall perform an analysis of the numerical viscosity of the kinetic method in this simplified framework.

Kinetic representation

The Lattice-Boltzmann Method (LBM) originated from the physical kinetic interpretation of the Navier-Stokes equations [START_REF] Chen | Lattice Boltzmann method for fluid flows[END_REF]. In the physical world, the fluid particles can have arbitrary velocities. The main idea of the LBM is that it is possible to construct abstract kinetic interpretations of the Navier-Stokes equations in which the particles velocities can take a few number of given values. This makes it possible to solve the kinetic model directly in an efficient way. We refer to [START_REF] Succi | The lattice Boltzmann equation: for fluid dynamics and beyond[END_REF] for a history of the LBM.

Initially devised for solving incompressible Navier-Stokes equations, the LBM is based on a single scalar particle distribution function. More recently, it has been extended to other systems of conservation laws. Extensions to MHD can be found in [START_REF] Dellar | Lattice kinetic schemes for magnetohydrodynamics[END_REF][START_REF] Martínez | Lattice Boltzmann magnetohydrodynamics[END_REF]. Dellar, in [START_REF] Dellar | Lattice kinetic schemes for magnetohydrodynamics[END_REF] showed that it is not possible to rely on a single scalar kinetic function for approximating the MHD equations. He proposes to represent the magnetic part of the equations with a vectorial kinetic formulation. Another vectorial kinetic interpretation of the MHD equations was also previously proposed in [START_REF] Croisille | Numerical simulation of the MHD equations by a kinetic-type method[END_REF] for building a numerical flux for a Finite Volume method.

The resulting hybrid kinetic model (scalar for the fluid part, vectorial for the magnetic part) is, however, limited to low-Mach number flows.

Several attempts have been done to extend the LBM to arbitrary Mach flows with a scalar kinetic function. For instance, in [START_REF] Frapolli | Entropic lattice Boltzmann model for compressible flows[END_REF][START_REF] Latt | Efficient supersonic flow simulations using lattice Boltzmann methods based on numerical equilibria[END_REF][START_REF] Frapolli | Theory, analysis, and applications of the entropic lattice Boltzmann model for compressible flows[END_REF] the Maxwell distribution is approximated with more velocities in the lattice. The additional degrees of freedom allow to computing a discrete Maxwellian distribution that matches exactly the moments of the compressible Euler equations. These approaches require quite large lattices. For instance, in [START_REF] Latt | Efficient supersonic flow simulations using lattice Boltzmann methods based on numerical equilibria[END_REF] the authors propose a D3Q39 LBM with 39 kinetic velocities for approximating the five-equation, three-dimensional, compressible Euler equations. In addition, the methods proposed in [START_REF] Frapolli | Entropic lattice Boltzmann model for compressible flows[END_REF][START_REF] Latt | Efficient supersonic flow simulations using lattice Boltzmann methods based on numerical equilibria[END_REF][START_REF] Frapolli | Theory, analysis, and applications of the entropic lattice Boltzmann model for compressible flows[END_REF] become unstable when the temperature exits from a given range. Another method is presented in [START_REF] Renard | Improved compressible hybrid lattice Boltzmann method on standard lattice for subsonic and supersonic flows[END_REF] for improving the stability of the scalar LBM for compressible flow. It relies on the classical D2Q9 or D3Q27 lattices. However, the LBM is supplemented with a finitedifference approximation of the temperature equation. This approach is thus not fully an LBM method.

For addressing transonic and supersonic flows, we think that a fully vectorial kinetic model is preferable. The general theory of lattice vectorial kinetic model is discussed in [START_REF] Bouchut | Construction of BGK models with a family of kinetic entropies for a given system of conservation laws[END_REF][START_REF] Aregba-Driollet | Discrete kinetic schemes for multidimensional systems of conservation laws[END_REF][START_REF] Graille | Approximation of mono-dimensional hyperbolic systems: A lattice Boltzmann scheme as a relaxation method[END_REF][START_REF] Dubois | Simulation of strong nonlinear waves with vectorial lattice boltzmann schemes[END_REF]. It is valid for any hyperbolic system of conservation laws and is no more limited to low-Mach number flows. It relies on a relatively small set of kinetic data. For instance, the LBM proposed by Dubois in [START_REF] Dubois | Simulation of strong nonlinear waves with vectorial lattice boltzmann schemes[END_REF] requires only 30 kinetic transported kinetic quantities compared to the 39 kinetic velocities of [START_REF] Latt | Efficient supersonic flow simulations using lattice Boltzmann methods based on numerical equilibria[END_REF]. In this work Dubois observes that while very robust and promising, his method is quite numerically diffusive. For obtaining accurate results, it is therefore necessary to consider refined meshes.

Let us finally mention that the vectorial kinetic approach is also of interest for low Mach numbers flows and the incompressible limit [START_REF] Bouchut | Second-order entropy satisfying BGK-FVS schemes for incompressible Navier-Stokes equations[END_REF][START_REF] Zhao | Discrete-velocity vector-BGK models based numerical methods for the incompressible Navier-Stokes equations[END_REF].

The vectorial kinetic approach can be used on arbitrary unstructured meshes at any order of approximation [START_REF] Badwaik | Task-based parallelization of an implicit kinetic scheme[END_REF][START_REF] Coulette | Palindromic Discontinuous Galerkin Method[END_REF]. When the lattice velocities are aligned with the mesh, it is possible to adopt the simple exact transport step of the LBM. This is the method that we present here.

We consider a real number λ > 0 and four velocities

v k = (v 1 k , v 2 k ) T , k = 1 . . . 4, defined by v 1 = -λ 0 , v 2 = λ 0 , v 3 = 0 -λ , v 4 = 0 λ .
We have four vectorial distribution functions f k (x, t) ∈ R m , k = 1 . . . 4. The conservative variable w is related to the kinetic data by

w = 4 k=1 f k .
Usually, in the Lattice-Boltzmann philosophy, the kinetic system is a set a transport equations coupled through a relaxation source term and reads (with τ > 0)

∂ t f k + v k • ∇f k = µ k = 1 τ (f eq k -f k ), k = 1 . . . 4. (3)
The equilibrium function f eq k also called the Maxwellian state is chosen as

f eq 1 (w) = 1 4 w - 1 2λ F 1 (w), f eq 2 (w) = 1 4 w + 1 2λ F 1 (w), f eq 3 (w) = 1 4 w - 1 2λ F 2 (w), f eq 4 (w) = 1 4 w + 1 2λ F 2 (w).
One can check that

w = 4 k=1 f eq k , F i = 4 k=1 v i k f eq k . (4)
When the relaxation time τ → 0 then from [START_REF] Barth | On the role of involutions in the Discontinuous Galerkin discretization of Maxwell and magnetohydrodynamic systems[END_REF] we see that f k f eq k , and thus from (4) we recover the MHD equations [START_REF] Aregba-Driollet | Discrete kinetic schemes for multidimensional systems of conservation laws[END_REF].

But in practice, it is quite difficult to solve equation (3) directly because the relaxation source term couples all the kinetic transport equation. It is better to replace the source term in such a way that most of the time it vanishes and thus the coupling is limited. For this, we introduce the Dirac comb Ψ defined by

Ψ (t) = i∈Z δ(t -i∆t),
where δ is the usual Dirac measure and ∆t a positive time step. Let ω be a relaxation parameter ∈ [START_REF] Aregba-Driollet | Discrete kinetic schemes for multidimensional systems of conservation laws[END_REF][START_REF] Badwaik | Task-based parallelization of an implicit kinetic scheme[END_REF]. The source term is then written as

µ k (x, t) = ωΨ (t) f eq k (w(x, t -)) -f k (x, t -) .
In other words, at times t = i∆t, f has jumps in time and the values after the jump

f k (x, t + ) = f k (x, t -) + ω f eq k (w(x, t -)) -f k (x, t -)
. This comes from the definition of the time derivative in the weak sense. The computation of the jumps at times t = i∆t are called the collision steps. The rest of the time, f k is a solution of the free transport equation

∂ t f k + v k • ∇f k = 0. (5)
The coupling indeed occurs only at times t = i∆t.

Let us remark that when ω = 1 the collision step reads

f k (x, t + ) = f eq k (w(x, t -)).
It is a simple projection on the Maxwellian state associated to the conservative data. This scheme is a first order in time approximation of the original equations. An interesting case is the case of over-relaxation ω = 2. Then the scheme is a second order in time approximation of the original equations. See for instance [START_REF] Dellar | An interpretation and derivation of the lattice Boltzmann method using Strang splitting[END_REF]). See also below for a mathematical analysis in the one-dimensional case.

Numerical method

For solving (3) numerically, we first construct a structured grid of the square D =]0, L[×]0, L[. The space step is given by ∆x = L N . and the grid-points are then as follows

x i,j = (i + 1 2 )∆x (j + 1 2 )∆x , i, j ∈ {0, . . . , N -1}.
For a simpler presentation we can assume periodic boundary condition, which amounts to the following equivalences

i + N ≡ i, j + N ≡ j.
We denote by w n i,j and f n k,i,j the approximation of w and f k at the grid points x i,j and time t n = n∆t just after the collision step. The values of the kinetic vectors just before the collision step are denoted f n,- k,i,j . For solving the kinetic system (3) we treat the transport and the relaxation terms separately.

Transport solver

The transport equation ( 5) admits an exact solution

f k (x, t + ∆t) = f k (x -∆tv k , t).
If we assume that the time step satisfies ∆t = ∆x λ , the transport operator then reduces to a simple shift. Before the collision step, we have thus

f n+1,- 1,i,j = f n 1,i+1,j , f n+1,- 2,i,j = f n 2,i-1,j , f n+1,- 3,i,j = f n 3,i,j+1 , f n+1,- 4,i,j = f n 4,i,j-1 . (6)

Relaxation

At the end of the transport step, we can compute the conservative data

w n+1 i,j = 4 k=1 f n+1,- k,i,j . (7)
The usual projection method would be to write

f n+1 k,i,j = f eq k (w n+1 i,j ). (8)
As stated above, it is more accurate to consider an over-relaxation

f n+1 k,i,j = 2f eq k (w n+1 i,j ) -f n+1,- k,i,j .
It is possible to add a small dissipation τ > 0. See [START_REF] Graille | Approximation of mono-dimensional hyperbolic systems: A lattice Boltzmann scheme as a relaxation method[END_REF] or [START_REF] Coulette | Palindromic Discontinuous Galerkin Method[END_REF] for the relaxation scheme in the case τ > 0. This dissipation provides a better stability in numerical cases, when sharp fronts are generated.

The relaxation step finally reads

f n+1 k,i,j = ωf eq k (w n+1 i,j ) -(ω -1)f n+1,- k,i,j , (9)
where ω = 2∆t 2τ +∆t . Typically, one will choose ω = 1.9 in our simulations. We see that the whole algorithm is extremely simple. It is a succession of "shifts" (6) and "collisions" [START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF]. For ω = 2 the scheme is second order but unstable in shocks. For ω = 1 the scheme is very robust, entropy dissipative, but quite diffusive. It would be interesting to construct a rigorous strategy for choosing locally the optimal value of ω. Let us mention that many finite volume schemes have been designed for solving MHD equations, with specific Riemann solvers. We can mention for instance the solver of [START_REF] Bouchut | A multiwave approximate Riemann solver for ideal MHD based on relaxation ii: numerical implementation with 3 and 5 waves[END_REF], based on a relaxation approach, with proven stability and accuracy features. However, the solver presented in [START_REF] Bouchut | A multiwave approximate Riemann solver for ideal MHD based on relaxation ii: numerical implementation with 3 and 5 waves[END_REF] is more complicated to program and less computationally efficient. We can thus compensate the slightly lower accuracy of the kinetic scheme by finer meshes, without increasing too much the computational load.

Boundary conditions

For the moment, we have assumed periodic boundary conditions. But Dirichlet or Neumann conditions are also possible. Neumann conditions are a good candidate for handling variables associated with outgoing waves. Indeed, because of the hyperbolic nature of the equations, it is not possible to impose all the physical data at the boundaries. The number of boundary conditions depends on the wave pattern. We refer to [START_REF] Drui | An analysis of over-relaxation in a kinetic approximation of systems of conservation laws[END_REF], for a specific analysis of the boundary conditions in the over-relaxation scheme.

Analysis of the numerical viscosity in the one-dimensional case

In this section, we state some results about the numerical viscosity of the kinetic relaxation scheme in the one-dimensional case. This one-dimensional analysis will give us simple intuitions for adjusting the relaxation parameter in the two-dimensional case. We consider the one-dimensional MHD-DC system (2). For more simplicity, we will denote F = F 1 and x = x 1 . The equations then read

∂ t w + ∂ x F(w) = 0.
For the analysis, it is possible to replace the scalar relaxation parameter ω by a matrix Ω, for more generality. In the following, we establish the comparison between matrices in the usual way, by the comparison of the associated quadratic form (the resulting order is thus not total).

It is then possible to prove the following result:

Theorem 1 If the relaxation matrix satisfies I < Ω < 2I and if f = f eq at the initial time, then, up to second-order terms in O(∆t 2 ), w is a solution of the following system of conservation laws

∂ t w+∂ x F(w) = λ 2 ∆t∂ x (Ω -1 - 1 2 I)(I - 1 λ 2 F (w) 2 )∂ x w +O(∆t 2 ). ( 10
)
Remark 1 The proof is based on standard Taylor expansions. For the scalar case Ω = ωI, it can be found (for instance) in [START_REF] Courtès | Vectorial kinetic relaxation model with central velocity. application to implicit relaxations schemes[END_REF]. The approach is classical in the analysis of the Lattice Boltzmann Method (LBM). See for instance [START_REF] Dubois | Equivalent partial differential equations of a lattice Boltzmann scheme[END_REF][START_REF] Otomo | Two complementary lattice-Boltzmannbased analyses for nonlinear systems[END_REF].

Remark 2 For Ω = I the proof, based on a Chapman-Enskog expansion, is also given in [START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF].

The above analysis allows recovering formally the so-called sub-characteristic condition. Assuming that I < Ω < 2I, the second-order ("viscous") terms have the good sign, which ensures stability of the model, if the following matrix is positive:

V(λ, w) = I - 1 λ 2 F (w) 2 > 0. (11)
A fully rigorous mathematical proof of stability of the kinetic model is given by Bouchut in [START_REF] Bouchut | Construction of BGK models with a family of kinetic entropies for a given system of conservation laws[END_REF], Section 3.2, pp. 140-142. Bouchut's proof is not based on asymptotic expansions but on fully non-linear entropy estimates. Let us emphasize that the vectorial kinetic construction ensures stability even when shock waves occur and is not limited by a low-Mach assumption.

From [START_REF] Chen | Lattice Boltzmann method for fluid flows[END_REF], we formally observe that the scheme is second-order accurate in time in the over-relaxation case, when

Ω = 2I.
Finally, this analysis provides a way to numerically approximate, up to second order in time, the second-order system of conservation laws

∂ t w + ∂ x F(w) -∂ x (A(w)∂ x w) = 0. ( 12
)
If the diffusion matrix A(w) is small, it is natural to take

Ω = 2 I + 2 λ 2 ∆t A(w)V(λ, w) -1 -1 . (13)
Let us remark that if λ is large enough then the above matrix is well defined.

In order to check practically the accuracy of the approximation, we apply the above analysis for a simplified system of two conservation laws.

We consider the one-dimensional isothermal Euler equations with a diagonal diffusion of ∂ xx (ρ, ρu)

T . Hence the full fluid system we are interested in is given by

∂ t ρ ρu + ∂ x ρu ρu 2 + c 2 ρ = ∂ xx ρ ρu . ( 14
)
This system is a very simplified version of the MHD equation, where the magnetic field is assumed to vanish and the gas is supposed to be isothermal. We have chosen to use this specific non-physical diffusion for our first test since it is exactly the type of diffusion that is apparent in a standard finite volume code using a Lax-Friedrichs flux. The diffusion matrix for ( 14) is simply given by

A(w) = 0 0 . (15)
From ( 13), the inverse of the relaxation matrix is given by

Ω -1 = a 1,1 a 1,2 a 2,1 a 2,2 (16)
with

a 1,1 = λ 2 -c 2 -3u 2 ∆t u 4 + (c 2 -λ 2 ) 2 -2u 2 (c 2 + λ 2 ) + 1 2 , a 1,2 = 2u ∆t u 4 + (c 2 -λ 2 ) 2 -2u 2 (c 2 + λ 2 ) , a 2,1 = - 2u(u -c)(c + u) ∆t u 4 + (c 2 -λ 2 ) 2 -2u 2 (c 2 + λ 2 ) , a 2,2 = u 2 + λ 2 -c 2 ∆t u 4 + (c 2 -λ 2 ) 2 -2u 2 (c 2 + λ 2 ) + 1 2 .
For our first test of the matrix relaxation, we solve ( 14) comparing the Lattice Boltzmann scheme using [START_REF] Courtès | Vectorial kinetic relaxation model with central velocity. application to implicit relaxations schemes[END_REF] as relaxation matrix and a standard explicit centered finite volume scheme for approximating [START_REF] Croisille | Numerical simulation of the MHD equations by a kinetic-type method[END_REF]. In this centered scheme, the time step is taken very small in such way that the stability condition is satisfied and that the time integration error can be neglected. In other words, the equivalent PDE of both schemes is [START_REF] Coulette | High-order implicit palindromic Discontinuous Galerkin method for kinetic-relaxation approximation[END_REF]. Hence, given the parameters in both schemes are set to represent the same diffusion , one should get the same type of diffusion for both schemes.

To test this we take the simple case of a stationary viscous shock. The initial data for this shock tube problem are

w l = (ρ l , ρ l u l ) T w r = ρ l u 2 l c , c u l T ( 17
)
where the sound speed c is set as c = 1 while for the left state we have chosen

w l = 2, 3 2 T . (18)
For the diffusion we use = 0.1. For the Lattice Boltzmann scheme, we set λ = 40. The results on a grid of 128 cells at time T = 0.1 are shown in Figure 1. The test shows that the matrix relaxation lattice Boltzmann scheme provides the correct diffusion and therefore results in the right viscous profile. The conclusion of this section is that the relaxation parameter is related to numerical viscosity. In 1D it can indeed be adjusted to fit the exact viscosity up to second order. For two-dimensional computations the analysis is more complicated, but heuristically we expect a similar behavior.

In the following, for simplicity reasons, we only consider scalar relaxation.

GPU implementation

We have implemented the above LBM algorithm. The LBM is particularly well adapted to parallelism. It is possible to provide very efficient implementations on GPU (Graphic Purpose Units) hardware.

OpenCL

Terminology

Today Graphic Processing Units (GPU) have more and more computation power. The Open Computing Language (OpenCL) is a software environment for simplifying the use of the GPUs for general computing. It is also possible to use OpenCL for driving a heterogeneous set of general multicore processors.

For giving an idea on how OpenCL is used in practice, we use the following terminology. An accelerator is a parallel computing device, such as a GPU or a multicore CPU. The host is the computer into which the accelerator is plugged. A kernel is a (generally small) program that is executed on several of the computing cores of the GPU. For instance, the Nvidia GPU GTX 1660 has 448 computing cores. Thanks to the OpenCL command queue management, it is possible to launch several million kernel instances, which are dispatched on the hundreds of cores of the GPU.

OpenCL means "Open Computing Language". The OpenCL runtime is a library of C functions, called from the host, in order to drive the GPU. The OpenCL runtime, because it is written in C is quite heavy to use in practice: the verbosity is high, the API is not very user-friendly and memory management is cumbersome. For this reason it is advised to use OpenCL wrappers written in a higher-level language such as C++ or Python. We have used the Python OpenCL wrapper written by Andreas Klöckner, PyOpenCL [START_REF] Klöckner | PyCUDA and PyOpenCL: A scripting-based approach to gpu run-time code generation[END_REF], which makes OpenCL initializations and calls much easier and shorter to program.

The OpenCL language is a C-like language for writing the kernels that will be executed on the computing cores.

OpenCL is practically available since September 2009. The specification is managed by the Khronos Group (that also drives the OpenGL specification). Several books describe today OpenCL in detail. We can refer, for instance, to [START_REF] Gaster | Heterogeneous computing with OpenCL: revised OpenCL 1[END_REF].

GPU

Very schematically we can consider that a general computing accelerator is made of global memory (typically 16 GB for the Nvidia Tesla V100) and compute units (typically 80).

Each compute unit is made of processing elements (typically 64), also called processors. A compute unit has a reserved local memory (typically 64 kb) shared by the processors of the unit.

A schematic picture of an abstract OpenCL accelerator is given in Figure 2

Programming rules

The same program can be executed on all the processing elements at the same time with the following rules:

1. All the processing elements have access to the global memory; 2. The processing elements have only access to the local memory of their compute unit; 3. The access to the global memory is relatively slow while the access to the local memory is very fast. 4. If possible, it is advised that the processing elements of the same compute unit access neighbour global memory locations, in order to improve "coalescence" (faster read/write access). 5. The memory transfers between the host memory and the GPU are slow and should be avoided; 6. If several processing elements try a read access at the same memory location (global or local) at the same time, all the reads will be successful; 7. If several processing elements try a write access at the same memory location (global or local) at the same time, only one write will be successful. For some hardware, atomic operations maybe available, but should be avoided for performance reasons.

In order to perform a complex task, a kernel has to be executed many times. Each execution of a kernel is called a work-item. A work-group is a collection of work-items running on the processing elements of a given compute unit. They can access the local memory of their compute unit. Each work-item is identified by a unique global ID p.

For more details on OpenCL, we refer for instance to [START_REF] Gaster | Heterogeneous computing with OpenCL: revised OpenCL 1[END_REF][START_REF] Klöckner | PyCUDA and PyOpenCL: A scripting-based approach to gpu run-time code generation[END_REF].

OpenCL implementation of the LBM algorithm

We have implemented the above algorithm using PyOpenCL.

The different devices used in this paper are listed in Table 1. The AMD processor was used inside a virtual environment, which implies a non-negligible loss of performance. PyOpenCL allows us to select either the CPU or the GPU for the computations. With the OpenCL AMD drivers, when the CPU is selected, it is also possible to choose the number of activated CPU cores through a Linux environment variable. This is useful for estimating (in a crude way) the efficiency of the OpenCL parallelism.

The LBM is very simple. Our implementation is made of two OpenCL kernels, a few C functions and a small Python driver for initializing the memory buffers, launching the OpenCL kernel and plotting the results. The role of the first OpenCL kernel is to compute the initial condition directly into the memory buffer created on the OpenCL accelerator. The role of the second kernel is to perform a time step of the LBM. The time-stepping is driven from the PyOpenCL program.

The most important point is to take care of the organizations of the kinetic data into memory. Indeed ensuring coalescent memory access is essential for performance. In practice all the values f k,i,j for kinetic velocities v k in cells (i, j) are arranged in a single memory buffer fn[], with the following storage

f k,i,j = fn[imem], imem=i+j*Nx+k*Nx*Ny,
where Nx is the number of grid points in the x 1 -direction and Ny is the number of grid points in the x 2 -direction. Then, in the OpenCL kernel, each cell (i, j) is associated to the work-item

p = i+j*Nx.
This numbering ensures that neighboring work-items will access neighboring memory locations during the shift algorithm [START_REF] Bouchut | Second-order entropy satisfying BGK-FVS schemes for incompressible Navier-Stokes equations[END_REF]. For instance, if work-item p access data in global memory, for reading or writing, at location imem, then work-item p = p + 1 access location imem+1.

This property is still true in the relaxation algorithm. First, the kinetic data are copied in processor registers in a coalescent way. The computations of the conservative data [START_REF] Bouchut | A multiwave approximate Riemann solver for ideal MHD based on relaxation ii: numerical implementation with 3 and 5 waves[END_REF], equilibrium (8) and collisions (9) are done in registers, which ensure very fast memory access. Finally, the kinetic data are copied back to the global buffer fn[] in a fully coalescent way.

We remark that thanks to the chosen organization into memory, we do not have to use the local cache memory for accelerating the algorithm.

In order to measure the efficiency of the implementation, we perform a memory bandwidth test for a 512 × 512 grid. One time-step of the method implies the read access in the global memory of the set of fields of the previous time-step. The local computations are done in registers. Then there is another write access to global memory for storing the data of the next time-step. The memory size in Gigabyte of one set of fields is

n GB = Nx×Ny×prec × 4 × m 1024 3 ,
where prec is the number of bytes for storing one floating point number (prec = 4 for single precision and prec = 8 for double precision). We then perform a given number of time iterations n iter and measure the elapsed time t elapsed (in s) in the OpenCL kernels. We perform two kinds of experiments.

In the first experiment, we deactivate the numerical computations and only perform the shift operations. The memory bandwidth (in GB/s) of the shift algorithm is then given by b = 2 × n GB × n iter t elapsed .

In the second experiment, we reactivate the computations and measure how the bandwidth is reduced. This allows to evaluating how the elapsed time is shared between memory transfers and computations. The results are given in Table 2. We observe a good efficiency of the shift algorithm in the shift-only case: the transfer rates are not very far from the maximal bandwidth of the device, at least for the GPU accelerators. From these results we also observe that the LBM algorithm is clearly memory bound. When the single precision computations are activated on the GPU devices (GTX, Quadro, V100), the elapsed time of the shift-and-relaxation test is not very different from the shift-only test. For the double precision computations, we observe that the V100 device outperforms all the other GPUs. 5 Numerical applications to MHD

Smooth vortex (performance test)

The smooth vortex test is a classical test for MHD codes. It is described for instance in [START_REF] Dumbser | A simple robust and accurate a posteriori sub-cell finite volume limiter for the Discontinuous Galerkin method on unstructured meshes[END_REF]. Because this is an exact solution, it allows us to assess the accuracy of the solver. Here we also used this test to evaluate the efficiency of the parallel implementation. The test case is built upon a single vortex, which is a stationary solution of the MHD system, to which a constant drift velocity is added. In the moving frame centered on r O (t) = tu drift , with u drift ∈ R 2 , the analytical solution reads in polar coordinates

ρ(r, θ) = ρ 0 , u(r, θ) = u 0 [u drift + h(r)e θ ], B(r, θ) = b 0 h(r)e θ , p(r, θ) = p 0 + b 2 0 2 (1 -h(r)), with b 0 = ρ 0 u 2 0 .
The results shown below are obtained with γ = 5/3 and the parameter set

ρ 0 = p 0 = 1, u 0 = b 0 = 0.2, u drift = (1, 1) T , h(r) = exp[(1 -r 2 )/2].
The computational domain is the square Ω =] -L/2, L/2[×]L/2, L2[, with L = 20. We compute the solution at time t = 10. The grid contains N points in x 1 and x 2 directions. For this smooth test case, we can take a relaxation parameter ω = 2. We compute the error e N in the L 1 norm at the final time between the exact solution and the numerical solution on the first component of the momentum (the other components of the solution would give similar results):

e N = Ω |(ρu 1 ) num -(ρu 1 ) exact | .
Asymptotically, the order of the scheme is evaluated by

β ln(e N /e 2N ) ln 2 .
The obtained numerical results are summed up in Table 3, where we give the convergence study and a performance evaluation of the implementation. OpenCL permits to run the same code on a multicore CPU or a GPU. We have tested several CPU or GPU hardware in single or double precision. Table 3 confirms the order of accuracy of the scheme in the case ω = 2. In addition we observe a good efficiency of the implementation on several types of GPU.

On CPU there is also a speedup achieved by the OpenCL parallelism, but it is very sensitive to the OpenCL drivers. For instance, with the same hardware (an Intel Xeon two-CPU system) the program runs almost three times faster with the OpenCL Intel drivers than with the open source POCL drivers. Let us finally mention that when it is run on only one core, the code is very slow. This is due to the fact that our implementation is not really optimized for correctly harnessing the CPU cache. Here, with a more clever tiling strategy, the one-core run could probably be accelerated by an order of magnitude (see for instance [START_REF] Helluy | Asynchronous OpenCL/MPI numerical simulations of conservation laws[END_REF]). Table 3 Convergence and performance study. Some tests are done in single precision (float32) and others in double precision (float64). The "efficiency" is a comparison for N=1024 with the slowest device. "CU" means "Compute Units": it is the number of activated cores in a CPU computation or of OpenCL compute units for a GPU computation.

CU precision N=128

Orszag-Tang vortex

The Orszag-Tang test case [START_REF] Orszag | Small-scale structure of two-dimensional magnetohydrodynamic turbulence[END_REF][START_REF] Dahlburg | Evolution of the Orszag-Tang vortex system in a compressible medium. i. initial average subsonic flow[END_REF][START_REF] Picone | Evolution of the Orszag-Tang vortex system in a compressible medium. ii. supersonic flow[END_REF] is often used to test a numerical method for MHD. It consists in a vortex system where turbulent structures and shocks develop. The domain is D =]0, 2π[×]0, 2π[ and the boundary conditions are periodic in x 1 = x and x 2 = y. The initial conditions are given in Table 4.

In Figures 3 and4 we present several snapshots of the evolution of the vortex. In Figure 5 we compare the results obtained with several grid refinements. The grid refinement clearly improves the sharpness of the shock profiles.

Tilt instability

The tilt instability has been studied in [START_REF] Richard | Magnetic reconnection driven by current repulsion[END_REF][START_REF] Strauss | An adaptive finite element method for magnetohydrodynamics[END_REF][START_REF] Lankalapalli | An adaptive finite element method for magnetohydrodynamics[END_REF]. The initial magnetic field configuration for tilt instability is a dipole current structure. It consists of two oppositely directed currents embedded in a uniform magnetic field (at large distance). These islands may be kept in an unstable equilibrium through a strong magnetic field. However, when perturbed, an instability develops that leads the islands to be aligned horizontally and to be expelled away. The tilt instability is illustrated for a simulation on a 1024 × 1024 grid in Figures 6 to 8 for different times and the numerical approach presented above. A detailed description of the different stages of this instability can be found in [START_REF] Richard | Magnetic reconnection driven by current repulsion[END_REF][START_REF] Baty | FINMHD: An adaptive finite-element code for magnetic reconnection and formation of plasmoid chains in magnetohydrodynamics[END_REF]. In [START_REF] Richard | Magnetic reconnection driven by current repulsion[END_REF][START_REF] Keppens | Interacting tilt and kink instabilities in repelling current channels[END_REF], the study of the tilt instability was also made for a compressible MHD system, while in [START_REF] Strauss | An adaptive finite element method for magnetohydrodynamics[END_REF][START_REF] Lankalapalli | An adaptive finite element method for magnetohydrodynamics[END_REF] the equations are incompressible, leading to slightly different results.

In the following, we will consider a squared spatial domain of dimensions D =]-3, 3[×]-3, 3[ and r 2 = x 2 +y 2 . The initial condition for the equilibrium is given by: ρ = 1.0, u = 0, p 0 = 1.0, ψ = 0,

k = 3.8317059702, K = -2 kJ 0 (k) , ϕ = KJ 1 (kr) y r , (19)
p = p 0 if r ≥ 1, p 0 + k 2 2 ϕ 2 else, (20)
B 1 = x 2 -y 2 r 4 -1 if r ≥ 1, K ky 2 r 2 J 0 (kr) + x 2 -y 2 r 3 J 1 (kr) else, (21)
B 2 = 2xy r 4 if r ≥ 1, -K kxy r 2 J 0 (kr) -2xy r 3 J 1 (kr) else, (22)
where J 0 and J 1 are the Bessel functions of the first kind of orders zero and one respectively. A perturbation of this equilibrium initially reads:

u = 2 exp (-r 2 )   -y x 0   , with = 10 -3 .
The asymptotic magnetic field strength for large r is unity, and thus defines our normalization. We also point out that we consider non-dimensioned equations in which the magnetic permeability is set to one. Consequently, our unit time is defined as the Alfven transit time across the unit distance (i.e. the initial characteristic length scale of the dipole structure).

Moreover, for the analysis of the simulation results, we recall that the current density is given by: C = ∇ × B, and that we are interested in the third component of C, i.e. C 3 . We will also look at the kinetic energy of the system, which is computed as the integral over the spatial domain of the kinetic energy density:

E kin = D 1 2 ρ|u| 2 .
We apply Dirichlet boundary conditions in this test case: the initial data in the boundary cells are simply kept constant. Finally, the numerical scheme involves the following parameters:

∆x = L x /N x , c h = 6 λ = 20, ∆t = ∆x/λ,
and ω = 1.9 for all variables except from the kinetic variables associated with ψ where ω = 1 (see section 5.3). In our numerical tests, we have noticed that choosing ω = 2 leads to an unstable numerical solution, due to the lack of dissipation.

Current sheets and current peak

Current sheets develop at the edges of the magnetic islands. The width of these current sheets as well as the maximum intensity of the current density depends on the numerical implementation, as studied in [START_REF] Lankalapalli | An adaptive finite element method for magnetohydrodynamics[END_REF], where specific adaptive resolution (adaptive grid and order of resolution) was performed near strong current gradients. In Figure 9, we show an enlargement of the part of the spatial domain that lies around the current sheet lying right side for simulations with different grid sizes and at time t = 6. The finest current structures are obtained with the finest grid (here the 4096 × 4096 grid). In Figure 10, we present the ratio of the maximum current density C 3,max over the initial current density C 3,0 . This ratio increases with the grid refinement and should be infinite in the case of a totally ideal MHD system ( [START_REF] Lankalapalli | An adaptive finite element method for magnetohydrodynamics[END_REF]). The maximum ratio we can obtain is 14.2 with the 7000 × 7000 grid. In comparison, [START_REF] Lankalapalli | An adaptive finite element method for magnetohydrodynamics[END_REF] obtained a ratio of 5 for the simulation with 32, 768 triangles and first-order method, and a ratio of 41 for adaptive grids near the current sheets. In conclusion, our method seems to have a higher numerical resistivity, which limits the maximum current peak, but the computations are probably faster. Indeed, adaptive methods are known to suffer from overhead due to the algorithm complexity.

Kinetic energy growth rate

During the development of the instability, the kinetic energy of the system increases exponentially. The growth rate of the kinetic energy has been studied in [START_REF] Richard | Magnetic reconnection driven by current repulsion[END_REF] and [START_REF] Lankalapalli | An adaptive finite element method for magnetohydrodynamics[END_REF], where different values have been found. In [START_REF] Richard | Magnetic reconnection driven by current repulsion[END_REF], two simulations with two values of β (the ratio between the hydrodynamic and the magnetic pressures) give growth rate values of 1.44 for low β case and of 1.27 for high β case. In [START_REF] Lankalapalli | An adaptive finite element method for magnetohydrodynamics[END_REF], the study is focused on the convergence of the kinetic energy growth rate according to the grid and the order of the numerical method. The converged value for the growth rate is near 1.3 for an initial perturbation = 1.0 • 10 -3 . The simulations of the tilt instability that we perform are post-processed at regular time intervals. In Figures 11 and12, we show the time evolution of the total kinetic energy E kin (t). One can note that the growth of the kinetic energy happens in two stages: in a first stage, the growth is faster than exponentially, then slows down and after two seconds the second stage begins with exponential growth. When the grid resolution is higher, the first stage tends to vanish and the second stage begins sooner.

In a logarithmic scale, the linear regression of the kinetic energy growth is performed with the Scipy function linregress and the regression line is illustrated in Figures 11 and12. The results of the linear regression are reported in Table 5 and Figure 12. The converged growth rate is close to 1.45 for our simulations. It is higher than the results of [START_REF] Richard | Magnetic reconnection driven by current repulsion[END_REF] and [START_REF] Lankalapalli | An adaptive finite element method for magnetohydrodynamics[END_REF]. In [START_REF] Keppens | Interacting tilt and kink instabilities in repelling current channels[END_REF] for a similar configuration, the growth rate is evaluated to 1.498. This is probably very close to the exact rate because the simulation is conducted with a very accurate adaptive scheme. Our results seem thus to be quite correct on the finest mesh.

In Table 5, we summarize the characteristics of the linear regression for each mesh size.

Divergence cleaning effect

In order to ensure a magnetic field which is close to a divergence-free field, we have presented the divergence cleaning procedure in Section 2.1. Here we compare two strategies for the numerical implementation of the divergence cleaning equation:

1. the kinetic variables associated with the conservative variable ψ are solved numerically like any other kinetic variable. This means that in the relax-ation step, the relaxation coefficient ω is set to ω = 1.9 and the order of resolution is close to 2, 2. during the relaxation step, the kinetic variables associated with the macrovariable ψ are relaxed with a coefficient ω = 1, which comes down to setting f n+1 k,i,j = ωf eq k (w n+1 i,j ) -(ω -1)f n+1,- k,i,j = f eq k (w n+1 i,j ).

Now, the order of resolution is 1 and the waves associated with the perturbations of the divergence-free constraint are better damped than with ω = 1.9. In [START_REF] Dedner | Hyperbolic divergence cleaning for the MHD equations[END_REF], several strategies are proposed for damping the divergence errors. One of the methods is to introduce a viscous damping term. In its spirit our proposal is similar, but the viscous term is introduced here in a numerical way. At time t = 0.5, one can notice the effects of the boundary conditions that have propagated towards the center of the domain. The perturbations of the divergence-free constraint have been more attenuated with the second strategy than with the first one. At time t = 1, the divergence-free constraint is mainly violated at the edges of the magnetic vortices and close to the domain boundaries. Once more, the results are better with the second divergence cleaning strategy. Finally, at time t = 5, when the vortices have begun to align, the divergence-free constraint is strongly perturbed, but results are still better with the second strategy.

Conclusion

In this work we have proposed a fast and robust Lattice-Boltzmann solver for the two-dimensional MHD equations with divergence cleaning. The method is general and can be extended to other systems of conservation laws.

We have provided a preliminary one-dimensional analysis in order to evaluate the numerical viscosity of the numerical scheme. We have shown that in principle, the numerical viscosity can be adjusted to a specific viscosity. An interesting challenge would be to extend the analysis to higher dimensions, but this leads to complex calculations because of cross second-order derivatives. If this analysis is possible, this would lead to a scheme where the physical resistive terms could be approximated properly.

The kinetic representation of the equations is very well adapted to massive parallel computing. The simplicity of the algorithm allows achieving almost optimal efficiency on GPU hardware. This leads to the possibility to conduct computations on very fine uniformly refined meshes.

When conducting the simulations, we observed that the limiting factor was memory rather than computation time. In order to extend the method to even finer meshes or to three-dimensional computations, the memory management has to be improved. An obvious possibility is to distribute the computations

Fig. 1

 1 Fig. 1 Viscous Shock Test with = 0.1: Comparison of the viscous profiles of the finite volume and the lattice Boltzmann scheme.

Fig. 2 A

 2 Fig. 2 A (virtual) GPU with 2 Compute Units and 4 Processing Elements

Table 4 Fig. 3

 43 Fig. 3 Snapshots of ρ for the Orszag-Tang configuration recorded at times t = 0.1 and t = 0.2. Grid size is N x × N y = 1024 × 1024.

Fig. 4

 4 Fig. 4 Snapshots of ρ for the Orszag-Tang configuration recorded at times t = 0.3 and t = 0.4. Grid size is N x × N y = 1024 × 1024.

Fig. 5

 5 Fig. 5 Snapshots of ρ for the Orszag-Tang configuration recorded at time t = 0.5 and t = 0.2. Grid sizes are N x × N y = 256 × 256, 1024 × 1024 and 4096 × 4096.

Fig. 6

 6 Fig. 6 Snapshots of the magnetic current density C 3 recorded at (non-dimensioned) times t = 1 and t = 2. Grid size is N x × N y = 1024 × 1024.

Fig. 7

 7 Fig. 7 Snapshots of the magnetic current density C 3 recorded at times t = 3 and t = 4. Grid size is N x × N y = 1024 × 1024.

Fig. 8

 8 Fig.8Snapshots of the magnetic current density C 3 recorded at times t = 5 and t = 6. Grid size is N x × N y = 1024 × 1024.

Fig. 9

 9 Fig. 9 Snapshots of the magnetic current density recorded at time t = 6 in a spatial zone around the current sheets. Grid sizes are N x × N y = 256 × 256, N x × N y = 1024 × 1024, N x × N y = 4096 × 4096.

Fig. 10

 10 Fig.10Ratio of the peak current and the initial current density for simulation grids ranging from 128 × 128 to 7000 × 7000.

Fig. 11

 11 Fig. 11 Growth of the kinetic energy during the simulation. Measures that have been accounted for the linear regression are identified with green marks and the regression line is in red. Results for grid sizes N x × N y = 512 × 512 and N x × N y = 2048 × 2048.

Fig. 12

 12 Fig.12Growth of the kinetic energy during the simulation. Measures that have been accounted for the linear regression are identified with green marks and the regression line is in red. Results for grid size N x × N y = 7000 × 7000. Growth rates from N x = N y = 128 (16384 degrees of freedom) to N x = N y = 7000 (49 × 10 6 degrees of freedom).

Fig. 13

 13 Fig. 13 Snapshots of the divergence of the magnetic field recorded at times t = 0.5 for both strategies of divergence cleaning. Grid size is N x × N y = 1024 × 1024.

Fig. 14

 14 Fig. 14 Snapshots of the divergence of the magnetic field recorded at times t = 1 for both strategies of divergence cleaning. Grid size is N x × N y = 1024 × 1024.

Fig. 15

 15 Fig. 15 Snapshots of the divergence of the magnetic field recorded at times t = 5 for both strategies of divergence cleaning. Grid size is N x × N y = 1024 × 1024.

Table 1

 1 Characteristics of the OpenCL devices tested in this paper. CU stands for "Compute Units" and "PE" for "Processing Elements".

	short name name	type frequency memory cache CU PE
	V100	Tesla V100-PCIE-16GB	GPU 1.4 GHz 16 GB 48 kB 80 5120
	Quadro Nvidia Quadro P6000	GPU 1.5 GHz 24 GB 32 kB 30 3840
	GTX	Nvidia GTX 1660	GPU 1.6 GHz 6 GB 48 kB 22 1408
	AMD	AMD EPYC 7551 32-core	CPU 2 GHz 48 GB 32 kB 1-24 1-24
	Intel	2 x Intel Xeon CPU E5-2609 v4 CPU 1.7 GHz 64 GB 32 kB 16 16
	Iris 640 Intel Iris Graphics 640	GPU 1.0 GHz 4 GB 64 kB 48 192

Table 2

 2 Bandwidth efficiency of the LBM algorithm. Comparison of the data transfer rates of the shift-only algorithm and of the shift-and-relaxation algorithm. The resulting bandwidth is compared with the maximal memory bandwidth advertised by the vendors of the hardware devices.

	Intel float32	17.58	13.38	60
	Intel float64	19.12	17.48	60
	Iris 640 float32	26.20	24.98	34
	Iris 640 float64	20.08	3.78	34
	GTX float32	147.54	146.94	192
	GTX float64	148.76	49.72	192
	Quadro float32	336.45	329.06	432
	Quadro float64	344.50	127.21	432
	V100 float32	692.31	676.44	900
	V100 float64	705.88	610.17	900

prec. b (GB/s, shift-only) b (GB/s, shift-relax) max. theoretical b (GB/s)

Table 5

 5 ). Characteristics of the linear regression for the kinetic energy growth rate for each mesh size: range of values that have been accounted for and growth rate.

	Mesh size	128 2	256 2	512 2	1024 2	2048 2	4096 2	7000 2
	Regression range (s) [5.6, 6.7] [4.9, 5.9] [4.3, 5.3] [3.9, 4.9] [3.7, 4.7] [3.0, 4.5] [3.0, 4.5]
	Growth rate	0.390	0.904	1.082	1.330	1.403	1.433	1.450

  These two options are compared in Figures 13 to 15 for different times. For obtaining those plots, the divergence of the magnetic field is computed with a simple centered finite difference approximation∇ • B B 1 (x + ∆x, y) -B 1 (x -∆x, y) 2∆x + B 2 (x, y + ∆y) -B 2 (x, y -∆y) 2∆y .

Observatoire de Strasbourg •

CEA Saclay •

IRMA, Université de Strasbourg, Inria Tonus, philippe.helluy@unistra.fr •

University of Würzburg

on several GPUs, using a task-based runtime system (such as in [START_REF] Bramas | Optimization of a discontinuous Galerkin solver with OpenCL and StarPU[END_REF]). This approach could also be mixed with a compression strategy in order to reduce memory occupation and data transfers.

Thanks: this work has been finalized during a one-week stay at the Oberwolfach Research Institute for Mathematics (MFO) in September 2020. The authors are very grateful to the staff of the MFO for their warm welcome and perfect organization despite the constraints caused by the pandemic.

The authors also thank Matthieu Boileau for his helpful contribution in the design of the PyOpenCL GPU code.


Etienne Decroly
email: etienne.decroly@univ-amu.fr

Jean-Michel Claverie

Bruno Canard

The WHO mission report struggles to trace the origins of the SARS-CoV-2 epidemic

.

Identifying the origins of the SARS-CoV-2 pandemic is a key issue for the prevention and response to future epidemics. After months of negotiations with the Chinese authorities, WHO could send 17 international experts for a joint study with Chinese experts to elucidate the origin of the virus that causes COVID-19. The report of the commission entitled "WHO Joint Global Investigation of the Origins of SARS-CoV-2: Chinese Part" was published in March 2021 [START_REF]Virus origin / Origins of the SARS-CoV-2 virus. WHO[END_REF]. Unfortunately, this study did not allow identifying the cause of the pandemic. However, the report proposes four hypotheses: 1) zoonotic origin followed by transmission via an intermediate host, 2) zoonotic origin with direct transmission from bats, 3) arrival in China via frozen food or the cold chain, and 4) laboratory accident. These working hypotheses are ranked from highly probable to very unlikely by the commission, without any rational explanation for this ranking and despite the fact the laboratory virus escapes have been previously documented [START_REF] Normile | SEVERE ACUTE RESPIRATORY SYNDROME: Lab Accidents Prompt Calls for New Containment Program[END_REF][START_REF] Heymann | Dangerous pathogens in the laboratory: From smallpox to today's SARS setbacks and tomorrow's polio-free world[END_REF]. The WHO Director-General Tedros Adhanom Ghebreyesus welcomed the work of the commission (whose objectivity was called into question in two open letters written by an international group of scientists) [START_REF]Calls for Further Inquiries Into Coronavirus Origins[END_REF][START_REF] Some | Inquiry Into the Coronavirus Pandemic's Origins[END_REF], but did not endorse the report conclusions and declared that all hypotheses remained on the table because the investigation had not been sufficiently thorough. Furthermore, it proposed to send a second expert group with broader powers, including biosafety specialists, with the aim of firmly establishing SARS-CoV-2 origin. Therefore, the investigation process must continue to understand how SARS-CoV-2 was transmitted to humans, probably from the virus reservoir discovered in bats from Southeast Asia.

Unfortunately, the controversy over the epidemic origin, which initially was a scientific issue, has become part of a political struggle between China and the United States and its closest allies. More than a year after the epidemic start, one might wonder whether it is still possible to solve this question objectively.

Let's go back to the basics, specifically the available robust scientific evidences that might help to inform the debate. First, the sequences of a growing number of SARS-CoV-2-related viruses have been obtained from bat and pangolin samples [START_REF] Wacharapluesadee | Evidence for SARS-CoV-2 related coronaviruses circulating in bats and pangolins in Southeast Asia[END_REF][START_REF] Lam | Identifying SARS-CoV-2-related coronaviruses in Malayan pangolins[END_REF][START_REF] Zhou | A pneumonia outbreak associated with a new coronavirus of probable bat origin[END_REF][START_REF] Zhou | A Novel Bat Coronavirus Closely Related to SARS-CoV-2 Contains Natural Insertions at the S1/S2 Cleavage Site of the Spike Protein[END_REF]. The analysis of these different viral sequences confirms the high recombination capacity of these coronaviruses, mainly in the genomic region that encodes the "Spike" (S) protein. This protein plays a key role in the recognition of the cell receptor Angiotensin-Converting Enzyme 2 (ACE2) and in the virus entry in the host cells [START_REF] Graham | Reservoirs, and the Modular Spike: Mechanisms of Coronavirus Cross-Species Transmission[END_REF][START_REF] Hu | Discovery of a rich gene pool of bat SARS-related coronaviruses provides new insights into the origin of SARS coronavirus[END_REF][START_REF] Sallard | Tracing the origins of SARS-COV-2 in coronavirus phylogenies: a review[END_REF]. The genomic sequences of these viruses display 89 to 96% of identity, which corresponds to more than 1200 mutations per genome or decades of biological evolution (in the absence of epidemics or recombination). The closest known ancestor (RaTG13) comes from a sample collected in 2013 in a mine in the Yunnan province, located more than 1500 km south from Wuhan, where six miners developed an atypical pneumonia form that led to the death of three of them [START_REF] Rahalkar | Lethal Pneumonia Cases in Mojiang Miners (2012) and the Mineshaft Could Provide Important Clues to the Origin of SARS-CoV-2[END_REF]. However, sequence comparison showed that RaTG13 is not the proximal ancestor of SARS-CoV-2, and the direct progenitor is still unknown. A second irrefutable observation is that the COVID-19 epidemic started in Wuhan City in late 2019. The sequences of virus samples collected at the beginning of the epidemic indicate very low genetic variability among these initial isolates [START_REF] Zhou | A pneumonia outbreak associated with a new coronavirus of probable bat origin[END_REF]. This suggests that the epidemic onset, or at least the spread of a virus with efficient human-to-human transmission, is likely to have occurred between September and October 2019.

The key remaining questions are to understand how the virus that initiated the epidemic spread in Wuhan City, and then to identify the missing links between the viruses circulating in bat populations and the human virus. It is currently believed that crossing the species barrier involves a sequence of events that allow bat viruses to: 1) acquire the ability to use the ACE2 receptors present at the surface of human cells; and 2) increase their inter-human transmissibility by acquiring a proteolytic cleavage site sensitive to furin, an ubiquitously expressed protease [START_REF] Coutard | The spike glycoprotein of the new coronavirus 2019-nCoV contains a furin-like cleavage siteabsent in CoV of the same clade[END_REF]. Usually, and fortunately, these two steps constitute a major bottleneck that limits zoonotic transmission and thus the pandemic risks.

From the beginning, the vast majority of the scientific community rapidly proposed zoonotic hypotheses because most viral emergences, including the latest SARS-CoV and MERS-CoV epidemics, are zoonoses for which the intermediate hosts have been quickly identified [START_REF] Sallard | Tracing the origins of SARS-COV-2 in coronavirus phylogenies: a review[END_REF][START_REF] Rahalkar | Lethal Pneumonia Cases in Mojiang Miners (2012) and the Mineshaft Could Provide Important Clues to the Origin of SARS-CoV-2[END_REF]. Furthermore, serological studies indicate that 0.6% of the population of Yunnan province has antibodies against SARS-like coronaviruses, compared to 0% in Hubei province, where Wuhan City is located [START_REF] Li | Serological Evidence of Bat SARS-Related Coronavirus Infection in Humans, China[END_REF]. This suggests that species transfer is possible, albeit not frequent, in the Yunnan. Fortunately, despite these seroconversions, it seems that the possibility of human-to-human transmission of these viruses is rare, because medical doctors did not observe any viral outbreak. This confirms the existence of a species barrier that limits human-to-human transmission. Therefore, the Chinese authorities tried to find evidence of zoonotic transmission (direct or via an intermediate host) in more than 80,000 samples from wild animals, livestock (35 species), and poultry collected in late 2019 and early 2020 [START_REF]Virus origin / Origins of the SARS-CoV-2 virus. WHO[END_REF]. All these samples were negative for viral RNA. Moreover, among the samples collected in Wuhan wet market, only environmental contamination by the human virus was detected. No intermediate animal host could be identified. These experimental observations do not support the main conclusion of the report according to which SARS-CoV-2 was most likely introduced into the human population via an intermediate host. However, lack of proof is not a proof, and it is still possible that SARS-CoV-2 virus passed through an intermediate species susceptible to the virus (e.g. cat, raccoon dog, mink, pangolin, civet cat) without detection of an epizootic outbreak. It should be noted that 15% of the cats tested in Wuhan after the epidemic had neutralising antibodies against SARS-CoV-2 [START_REF] Zhang | A serological survey of SARS-CoV-2 in cat in Wuhan[END_REF].

The WHO report proposes an alternative hypothesis by stating that when SARS-CoV-2 was circulating widely in the human population in 2020, traces of this virus were found on frozen meat packages. There is little scientific evidence to support this hypothesis, and to our knowledge, there is no documented coronavirus outbreak through frozen meat to date. Yet, this hypothesis could explain how this outbreak started in a megacity of eleven million people where potentially infected animals are rare. It is a politically satisfactory hypothesis for the Chinese authorities, but scientifically weak, without robust evidence. Moreover, it just displaces the problem, without solving the question of the origin.

Finally, the report disqualifies the hypothesis of an accidental origin linked to experimental work on coronaviruses in different laboratories in Wuhan City. It is important to note that the priority of the WHO-China joint study was to investigate a zoonotic origin and not to evaluate all possible sources of the pandemic, and that the mandate negotiated by the WHO does not mention this last hypothesis [START_REF]Virus origin / Origins of the SARS-CoV-2 virus. WHO[END_REF]. This hypothesis is based particularly on the fact that the RaTG13 virus, which is currently the closest known neighbour of SARS-CoV-2, was sampled by a laboratory located close to the area where the first cases of SARS-CoV-2 infection were detected. Research projects carried out at the Wuhan Institute of Virology (and co-funded by the NIH [22]) wanted to understand the mechanism of species barrier crossing used by coronaviruses. These studies involved collecting viruses from bats, their sequencing to estimate their diversity, their culture (if possible) using different cell types to elucidate the mechanism of species barrier crossing, and their propagation in various animal models, including mice that express the human ACE2 receptor [START_REF] Hu | Discovery of a rich gene pool of bat SARS-related coronaviruses provides new insights into the origin of SARS coronavirus[END_REF]22]. Moreover, this institute has developed strategies to construct chimeric viruses by exchanging receptor recognition domains, as documented by many of their publicly available publications and research projects [START_REF] Hu | Discovery of a rich gene pool of bat SARS-related coronaviruses provides new insights into the origin of SARS coronavirus[END_REF][START_REF] Zeng | Bat Severe Acute Respiratory Syndrome-Like Coronavirus WIV1 Encodes an Extra Accessory Protein, ORFX, Involved in Modulation of the Host Immune Response[END_REF][START_REF] Menachery | A SARSlike cluster of circulating bat coronaviruses shows potential for human emergence[END_REF]. The aims of these gain-of-function experiments were to identify the molecular determinants that promote infection of human cells, and to develop vaccine strategies against this potentially pandemic family of viruses. These research programmes are entirely compatible with the hypothesis of a laboratory accident, for example through the contamination (even asymptomatic) of a staff member, or due to a negligence in biohazardous waste management. This hypothesis would also explain the emergence of this virus in a city that concentrates the largest centres for coronavirus studies in the world, and where seven separate laboratories (of biosafety level BSL-2 to 4) handle these viruses. It should also be noted that a growing number of laboratories in the worldcivilian and undoubtedly military-are studying bat viruses, particularly from samples taken in Asia. The hypothesis of a leak from another lab other than the one in Wuhan and its importation, although highly unlikely, cannot be excluded. In the absence of detailed information on the activities of these laboratories, the "accidental leak" is a hypothesis that has not been rigorously evaluated yet.

An important observation, not documented in the WHO report but much discussed in the context of accidental leakage, is the presence of a basic amino acid-rich cleavage site in the S protein at the junction between the S1 and S2 domains [START_REF] Coutard | The spike glycoprotein of the new coronavirus 2019-nCoV contains a furin-like cleavage siteabsent in CoV of the same clade[END_REF][START_REF] Johnson | Loss of furin cleavage site attenuates SARS-CoV-2 pathogenesis[END_REF]. This proteolytic cleavage site is absent in all currently sequenced SARS-like coronavirus strains. It plays a key role in the human-to-human transmission of the virus and in its ability to infect a broad spectrum of cells. Indeed, the cleavage between the S1/S2 domains, called "priming", induces a conformational change in the S protein that is required for the optimal recognition of the human ACE2 receptor [START_REF] Wrobel | SARS-CoV-2 and bat RaTG13 spike glycoprotein structures inform on virus evolution and furin-cleavage effects[END_REF]. The origin of this site, which results from the insertion of a 12nucleotide sequence in the coding sequence for the S protein, is discussed in the literature. The supporters of the zoonotic origin propose that the natural insertion of this site is linked to "template switching" mechanisms frequently observed in the case of co-infection by two coronaviruses harbouring sufficiently homologous sequences [START_REF] Gallaher | A palindromic RNA sequence as a common breakpoint contributor to copy-choice recombination in SARS-COV-2[END_REF]. On the other hand, the supporters of the laboratory accident argue that many experiments involving the insertion of sites rich in basic amino acids at the S1/S2 junction have already been carried out to potentiate infection by coronaviruses [START_REF] Cheng | The S2 subunit of QX-type infectious bronchitis coronavirus spike protein is an essential determinant of neurotropism[END_REF][START_REF] Belouzard | Activation of the SARS coronavirus spike protein via sequential proteolytic cleavage at two distinct sites[END_REF]. They are also surprised by the unexpected presence of two successive arginine codons (CGG), rarely observed in this family of viruses, but classically used in codon re-encoding experiments to promote protein expression in eukaryotic cells. The origin of the furin cleavage site is not clear yet, and more virus sequences are needed to trace SARS-CoV-2 origin, including the SARS-CoV-2 precursor without this cleavage site. Unfortunately, entry to the Yunnan cave is now forbidden, and new sampling is impossible (including by the WHO mission).

Although SARS-CoV-2 origin is still unknown, elucidating the origin of this pandemic remains an essential challenge to adapt measures in order to prevent the emergence of new pandemics. If the zoonotic origin results from the passage through an intermediate host, it will be necessary to take measures to monitor potentially infected animals that might also be candidate vectors of these viruses. This possibility completely shapes our relationship with wild areas and their fauna. For humans, a health risk is the combination of a danger (a microbe that we cannot control) and human exposure (that we can regulate). Ensuring that biodiversity-rich areas are not overly "anthropized" is both necessary for their preservation and to avoid new pandemics. The consumption of bushmeat (and its poaching) must be more closely monitored, and perhaps banned. If the COVID-19 pandemic turns out to be the result of a laboratory accident, then the experimental conditions under which potentially pathogenic viruses are collected and studied must be better regulated, particularly "gain-offunction" experiments [START_REF] Butler | Engineered bat virus stirs debate over risky research[END_REF][START_REF] Casadevall | H5N1 moratorium controversy and debate[END_REF][START_REF] Evans | The ethics of biosafety considerations in gainoffunction research resulting in the creation of potential pandemic pathogens[END_REF]. Although conferring new properties to a virus seems sometimes necessary to understand its pathogenicity, we must ask how the benefit/risk balance of such experiments should be evaluated. Is it reasonable to generate in laboratories viruses that might cause pandemics in humans and do not exist naturally? This ethical debate on "gain-of-function" experiments has been going on since 2011, when American and Dutch researchers wanted to construct potentially pandemic influenza viruses from a virus that was initially poorly adapted to transmission by aerosol [START_REF] Imai | Experimental adaptation of an influenza H5 HA confers respiratory droplet transmission to a reassortant H5 HA/H1N1 virus in ferrets[END_REF][START_REF] Russell | The potential for respiratory droplet-transmissible A/H5N1 influenza virus to evolve in a mammalian host[END_REF]. The main question is whether the expected benefit of these experiments can justify taking the risk, or whether they should be considered too dangerous and therefore banned. This dilemma led the United States to declare a moratorium on this type of experiment from 2014. These trade-offs are complex, and it is necessary to assess the potential risks and benefits of all experiments in order to define limits without amputating virology research. We should also remember that ignorance does not eliminate the danger. Natural biological samples do not necessarily and always end up in high-security laboratories. This might happen only after they are found to contain potentially pandemic viruses. Surprisingly, at a time when voices are being raised about the risk of accidents in the few existing BSL4 laboratories, [START_REF] Van Boeckel | The Nosoi commute: a spatial perspective on the rise of BSL-4 laboratories in cities[END_REF] China has just announced plans to build 23 BSL4 (one per province) and 88 BSL3 laboratories [START_REF]China goes on biosafety lab building spree[END_REF]. While waiting for the (unlikely) establishment of an independent international organisation to monitor their activities (which could be modelled on the International Atomic Energy Agency), an emergency solution could be to equip each of these laboratories with a sort of biological "black box" (e.g. an archive of all air filters and electronic laboratory notebooks) that would allow identifying (or exonerating) them as source of accidents or problems.

In conclusion, it is also important to note that the introduction of the mRNA vaccine technology has made it much easier to simultaneously construct a dangerous virus and its vaccine. This technological breakthrough has led to a sudden change of era. The possibility to concomitantly possess these two entities (a virus and its vaccine) introduces a radically new political, strategic and military dimension to which neither virologists nor society should remain indifferent.


M Soto-Crespo

J M Dudley

Nail Akhmediev

Adrian Ankiewicz

J M S O T O -C R E S P O B,

Rogue wave early warning through spectral measurements?

de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Rogue waves in the ocean are well known to have disastrous consequences [1,2]. They appear both in the deep ocean and in shallow waters [START_REF] Akhmediev | Rogue Waves -Towards a Unifying Concept[END_REF]. In contrast to tsunamis and storms associated with typhoons that can be predicted hours (sometimes days) in advance, the particular danger of oceanic rogue waves is their sudden appearance as "waves that appear from nowhere" only seconds before they hit a ship. The grim reality, however, is that although their existence has now been confirmed by multiple different means, there remains uncertainty about their fundamental origins which hinders systematic attempts to study their characteristics and perhaps predict their appearance [START_REF]Extreme Ocean Waves[END_REF].

Clearly, large waves in the ocean can arise from multiple causes, but there is an increasing consensus that nonlinear wave shaping effects play a central role under many conditions to increase the wave amplitude. This idea has recently received much attention motivated by studies in optical fiber [START_REF] Dudley | [END_REF]6] where the dynamics are described by the same nonlinear Schrödinger equation (NLSE) as proposed to describe the group envelope of deep water waves. As a basis, if we take the model where rogue waves are described by the NLSE, then one of the prototypes of the rogue wave is the Peregrine soliton [7][8][9][START_REF] Ten | Simulation of the ocean waves and appearance of freak waves[END_REF][START_REF] Clamond | [END_REF][12]. Higher-order rational solutions also exist [13][14][15] and they can also be associated with extreme amplitudes and rogue wave characteristics [16,17]; indeed, they possess higher amplitudes than the Peregrine soliton.

Despite the importance of these rational solutions in NLSE theory and in studies of rogue waves, they have surprisingly been the subject of limited attention. Recently, however, optical experiments have studied the excitation conditions of the Peregrine soliton in detail, confirmed its existence experimentally, and explicitly characterised its two-dimensional localisation [18]. Of course, previous experiments in wave tanks have studied related nonlinear wave propagation mechanisms [19][START_REF] Huijsmans | Experiments on extreme wave generation using the soliton on finite background[END_REF][START_REF] Karjanto | [END_REF][START_REF] Osborne | Nonlinear Ocean Waves and the Inverse Scattering Transform[END_REF], but the experiments in optics have enabled the measurement of unambiguous signatures of the defining features of the Peregrine soliton.

A particular aspect of experiments in optics is the ease of measuring optical spectra of pulse envelopes using readily-available instruments such as grating spectrometers. This has led naturally to renewed interest into the frequency-domain characteristics of the Peregrine soliton, and spectral measurements of the soliton pulse envelope have revealed a characteristic triangular form when plotted on a log scale. Significantly, theory predicts that the development of the triangular form of the envelope spectrum is an important and universal feature of parametric nonlinear spectral broadening, and indeed experiments studying the emergence of spontaneous nonlinear breather pulses from modulation instability in optics have also exhibited triangular spectra in agreement with theory [START_REF] Dudley | [END_REF].

These results from optics have naturally raised the question whether such spectral measurements would be possible in the ocean, and also whether such observations could provide important signatures of nonlinear effects favoring rogue wave formation. In fact, analysis of the spectral content of the full spectrum of waves on the surface of the ocean has already been studied in the context of the nonlinear wave propagation, and the Benjamin-Feir Index has been introduced as a measure of how the frequency content of the wave spectrum relates to the frequency band of instability gain [START_REF] Osborne | Nonlinear Ocean Waves and the Inverse Scattering Transform[END_REF]24,25]. This can indeed be applied in a significant forecasting context to identify sea state conditions that may be associated with an elevated probability for rogue wave appearance, providing important information allowing avoidance of particular ocean regions when mapping navigation routes. This approach, however, is based on analysis of the spectral characteristics of the sea-state as a whole. Motivated by the experiments in optics discussed above, however, we focus here on how the spectral characteristics of the envelope of propagating waves could be used to reveal the presence of nonlinearity and rogue waves a short time before the appearance of a particular rogue wave event. That is, we are interested in particular signatures of the ocean wave dynamics that can be used as "early warning" indicators in all already agitated sea. Specifically, by considering real-rime spectral analysis of the wave envelope, the appearance of characteristic triangular spectral decay of the wings could provide valuable advance and immediate warning of nonlinear amplification conditions and the emergence of rogue waves. If such spectral signatures were apparent before visual observation, then this could be extremely important in allowing appropriate evasive action to be taken.

Significantly, there is already important evidence linking such spectral decay and rogue wave formation on a smaller scale with capillary waves [26]. The open ocean, of course, presents additional difficulties, and precisely how spectral measurements could be accomplished in practice through sensing and time-series measurements of ship displacement is beyond the scope of this Letter. However, it is clearly an important question to address, and our goal here is to make the first steps in this direction. In particular, we review the properties of the characteristic triangular spectrum for rational soliton solutions of the NLSE and demonstrate that such characteristics are universal, also being observed for higherorder rogue wave solutions. By showing how these spectral characteristics appear at an early stage of evolution, our analysis supports the idea that analysing the shape of propagating waves in the frequency domain may provide important information about the emergence phase of nonlinear wave shaping mechanisms.

Spectra of rogue waves

There are several models for the description of rogue waves, both in optics and marine studies. In the lowest-order approximation, waves in the deep ocean are described by the nonlinear Schrödinger equation (NLSE). Here we start the analysis with the NLSE in dimensionless form:

iu z + 1 2 u tt +|u| 2 u = 0, ( 1
)
where z is the longitudinal distance, t is the retarded time in the frame moving with the group velocity of the waves and |u(x, t)| is the wave envelope. This notation is valid for both fibre optics and ocean wave theory.

One of the NLSE solutions is the Peregrine soliton, suggested as a prototype of rogue wave as early as 1983 [7]. In our scaling, it is given by:

u(z, t) = 4 1 + 2it 1 + 4t 2 + 4z 2 -1 e iz . ( 2
)
It is a solution that is localised both along z and t directions (see upper left panel in Fig. 1). There is no periodicity and the spectrum must be a continuous function of frequency, rather than having a comb structure.

The Fourier transform of Eq. ( 2) can be easily calculated:

F (ω, z) = 1 √ 2π ∞ -∞ u(z, t)e iωt dt = √ 2π 1 + 2iz √ 1 + 4z 2 exp - |ω| 2 1 + 4z 2 -δ(ω) exp(iz), (3)
where the Dirac delta function, δ(ω), appears due to the finite background level. The modulus of this spectrum is given by:

F 1 (ω, z) = √ 2π exp - |ω| 2 1 + 4z 2 . ( 4
)
The delta function is omitted here, so Eq. ( 4) represents the spectrum of the variable part of the solution. This spectrum is shown in Fig. 1 on a log scale. The width of the spectrum evolves from zero at z →-∞ to a finite value at z = 0 and decreases again at z →∞. Remarkably, at any particular z, the spectrum has a triangular shape and the experiments in optics referred to previously have now confirmed this prediction quantitatively [18].

The spectra of higher-order rogue waves

A second order rogue wave solution was first presented in [13] (see also [16]). This solution can be considered as a nonlinear superposition of two Peregrine solitons. There is a whole hierarchy of these solutions, with increasing order, which can be constructed using various tools. In principle, the hierarchy of higher-order polynomial solutions of the NLSE can be written in terms of Wronskians [14]. Although quite general, solutions in this form still need to be simplified, because the determinants generally contain extra factors in the form of polynomials that have to be separated and cancelled, and this is not a trivial task, even for the second order solutions. On the other hand, for any order j, the most convenient explicit way to write the polynomial solutions is the following [16]:

u j (z, t) = G j (z, t) + izH j (z, t) D j (z, t) + (-1) j e iz , (5)
where G 2 =-

3 16 + 3 2 t 2 + t 4 + 9 2 z 2 + 6t 2 z 2 + 5z 4 = t 2 + z 2 + 3 4 t 2 + 5z 2 + 3 4 - 3 4
,

H 2 =- 15 8 z -3t 2 z + 2t 4 z + z 3 + 4t 2 z 3 + 2z 5 = z z 2 -3t 2 + 2 t 2 + z 2 2 - 15 8
,

D 2 = 3 64 + 9t 2 16 + t 4 4 + t 6 3 + 33 16 z 2 - 3 2 t 2 z 2 + t 4 z 2 + 9 4 z 4 + t 2 z 4 + z 6 3 = 1 3 t 2 + z 2 3 + 1 4 t 2 -3z 2 2 + 3 64 12t 2 + 44z 2 + 1 .
With the usual Fourier transform [FT] definition, Eq. ( 3),w ec a n obtain the spectra of higher-order rogue waves, Eq. ( 5).E a c hs p e ctrum contains a delta-function due to the constant background part of the solution. We can eliminate the delta-function by taking the Fourier transform of the varying part. Hence, we take the spectrum of

ψ(z, t) = u j (z, t) exp[-iz]-(-1) j = G j (z, t) + izH j (z, t) D j (z, t) . (6)
In contrast to u(z, t), the function t) approaches zero for large |z| and |t|.

We define F j (ω, z) as the FT of the j-th rogue wave, i.e. the FT of the function given in Eq. ( 6). It is easy to find that the central value at ω = 0 and z = 0i s : using the definition of [17] for N j . The above integral can be found using zeros in the upper half plane of complex t [15].U s i n gt h i s , we find that

F j (0, 0) = 1 √ 2π ∞ -∞ G j (0, t) D j (0, t) dt = N j √ 2π , (7)
F j (0,0) √ 2π
is equal to ( j + 1)/2 when j is odd, and equal toj/2 when j is even. Thus, for the lowest-order cases,

F j (0,0) √ 2π equals (1, -1, 2) when j = (1, 2 , 3), respectively.
For j higher than one, general analytic expressions are complicated, but we can readily illustrate the results graphically. Fig. 2 shows the spectrum for the 2nd order rogue wave, along with the original profile of |u(z, t)|. The spectrum initially (z =-4) has a shape close to being triangular. It appears to have a dip at z = 0, but the details in the above right panel of Fig. 2 show that it has a more complicated structure at this point. Indeed, the second-order rogue wave itself has a more involved shape, and therefore so does its spectrum. On a log-scale, the spectrum has two dips, due to the fact that F 2 (ω, 0) is zero at two values of ω. Other than that the general structure of the spectrum is roughly triangular.

Explicit results for the amplitudes of higher-order rogue wave solutions have been presented in [17]. The expressions for the third-order solution are complicated and are not repeated here. Instead, we directly plot the spectrum for the 3rd order rogue wave in Fig. 3. It has a structure which is, in a way similar to, but more complicated than, the second-order result.

The results presented in this section may become useful as the study of rogue waves continues to develop. When direct measurements of wave profiles are difficult, measuring the spectral profiles can help in identifying the presence of rogue waves. As rogue waves tend to appear in various branches of physics [START_REF] Akhmediev | Rogue Waves -Towards a Unifying Concept[END_REF]26],w e anticipate that developing a unified approach to their study may become extremely useful.

Conclusion

We have shown that the Peregrine soliton, which can serve as a prototype of rogue waves in the ocean has a triangular spectrum at every stage of its development. The spectra of higher order rogue waves are more complicated but can also be calculated, and we have presented spectra for the three lowest order solutions. Although more structured than the spectra for the Peregrine soli-ton, these higher-order solutions present characteristics that are approximately triangular when away from the zero frequency position.

A particular consequence of the results we present here relates to the insight we provide into the evolution of the spectra of these solutions. This is because, within an oceanographic context, our results suggest the potentially extremely important prospect of identifying spectral signatures of the early emergence of rogue waves resulting from nonlinear wave shaping. By linking the appearance of triangular spectra to the generation of high amplitude localised rogue wave structures, our results may provide a fruitful direction for further research aiming to identify characteristic signatures of the conditions from which rogue waves might emerge.

Fig. 1 .

 1 Fig. 1. (Color online.) (Upper left panel): Fundamental ( j = 1) r o g u ew a v e( P e r e g r i n e soliton). (Bottom right panel): The spectrum of the fundamental rogue wave, plotted on a log scale i.e. log |F 1 (ω, z)|. The always-present delta function at ω = 0i s omitted here.

Fig. 2 .

 2 Fig. 2. (Color online.) (Above left panel): second-order ( j = 2) rogue wave. (Bottom panel): spectrum of second order ( j = 2) r o g u ew a v eo nl o gs c a l e ,l o g|F 2 |,d u r i n g evolution along the z-axis. (Above right panel): same spectrum on the central line, z = 0.

Fig. 3 .

 3 Fig. 3. (Color online.) (Above left panel): third-order ( j = 3) rogue wave [17].( B o ttom panel): spectrum of the third order ( j = 3) rogue wave on a log scale, log |F 3 |, showing evolution along the z-axis. (Above right panel): same spectrum on the central line, z = 0.

Acknowledgements

N.A. is grateful to Alexander von Humboldt Foundation as this work was prepared while he visiting Germany. A.A. acknowledges the support of the Australian Research Council (Discovery Project DP0985394). J.M.S.C. acknowledges support from the Spanish Ministerio de Ciencia e Innovación under contracts FIS2006-03376 and FIS2009-09895. J.M.D. thanks the Institut Universitaire de France and the French Agence Nationale de la Recherche projects MANUREVA ANR-08-SYSC-019 and IMFINI ANR-09-BLAN-0065, for support.