Asma Toumi
email: asma.toumi@cmla.ens-cachan.fr

Guillaume Dufour
email: guillaume.dufour@onera.fr

Ronan Perrussel
email: perrussel@laplace.univ-tlse.fr

Thomas Unfer
email: thomas.unfer@laplace.univ-tlse.fr

Second order accurate asynchronous scheme for modeling linear partial differential equations ✩

We propose an asynchronous method for the explicit integration of multi-scale partial differential equations. This method is restricted by a local CFL (Courant Friedrichs Lewy) condition rather than the traditional global CFL condition. Moreover, contrary to other local time-stepping (LTS) methods, the asynchronous algorithm permits the selection of independent time steps in each mesh element. We derived an asynchronous Runge-Kutta 2 (ARK2) scheme from a standard explicit Runge-Kutta method and we proved that the ARK2 scheme is second order convergent. Comparing with the classical integration, the asynchronous scheme is effective in terms of computation time.

Introduction

Numerical simulation has become a central tool for the modeling of many physical systems (fluid dynamics, plasmas, electromagnetism, etc.). Multi-scale phenomena make the integration of these physical systems difficult in terms of accuracy and computation time. Time-stepping integration techniques used for solving such problems generally fall into two categories: explicit and implicit schemes. In the explicit schemes, all unknown variables are computed at the current time level from quantities already available. Time step is then limited by the most restrictive CFL condition over the whole computation domain. The implicit method allow to overcome this time step constraint but comes with an additional computational cost arising from solving large linear systems. Thus the implicit approach becomes effective in terms of computational cost only if the time steps used are quite large with respect to the CFL condition. However this is not always possible, especially when strongly coupled and multi-scale phenomena such as microwave plasma generation are studied.

To overcome such problems, a number of local time-stepping approaches have been developed. These methods are restricted by a local CFL condition rather than the traditional global CFL condition. In [START_REF] Osher | Numerical approximations to nonlinear conservation laws with locally varying time and space grids[END_REF], Osher and Sanders proposed a local time-stepping scheme for one-dimensional scalar conservation laws. They gave a thorough analysis of a first order spatial discretization with a local forward Euler time-stepping scheme. This scheme allows each element to take either an entire time step or some fixed number of smaller steps. Tang and Warneck proposed in [START_REF] Tang | A class of high resolution schemes for hyperbolic conservation laws and convection-diffusion equations with varying time and space grids[END_REF] a class of high resolution local time step schemes for hyperbolic conservations by projecting the solution increments at each local time step. Savcenco et al. [START_REF] Savcenco | A multirate time stepping strategy for stiff ordinary differential equations[END_REF] constructed a multirate scheme for parabolic problems. This scheme was obtained by adaptation of an implicit Rosenbrock-type scheme. The idea is to compute a prediction for all the cells of the mesh with the trapezoid formula. Then the fine mesh is updated using information from the coarse mesh by quadratic interpolation. Dawson and Kirby [START_REF] Dawson | High resolution schemes for conservation laws with locally varying time steps[END_REF], developed upwind methods for solving conservation laws which allow local time refinement to be coupled with local spatial refinement. In that scheme a limiter is applied which is adapted to the outcome of previous stages. W. Hundsdorfer et al. [START_REF] Hundsdorfer | Monotonicity conditions for multirate and partitioned explicit Runge-Kutta schemes[END_REF], noted that several multirate schemes for conservation laws have one of the following defects: there are schemes that are locally inconsistant, e.g. [START_REF] Osher | Numerical approximations to nonlinear conservation laws with locally varying time and space grids[END_REF][START_REF] Dawson | High resolution schemes for conservation laws with locally varying time steps[END_REF], and schemes that are not mass-conservative, e.g. [START_REF] Tang | A class of high resolution schemes for hyperbolic conservation laws and convection-diffusion equations with varying time and space grids[END_REF][START_REF] Savcenco | A multirate time stepping strategy for stiff ordinary differential equations[END_REF]. They discussed these two defects for one-dimensional conservation laws. In the same context, an error analysis is presented in [START_REF] Hundsdorfer | Error analysis of explicit partitioned Runge-Kutta schemes for conservation laws[END_REF], for explicit partitioned Runge-Kutta methods and multirate methods applied to conservation laws. Hundsdorfer and al. showed that the different multirate methods studied in [START_REF] Hundsdorfer | Monotonicity conditions for multirate and partitioned explicit Runge-Kutta schemes[END_REF] lead to order reduction of the schemes. To guarantee mass conservation, flux-based decompositions are studied. In this case, the accuracy may deteriorate. Another approach developed by Berger and Oliger [START_REF] Berger | Adaptive mesh refinement for hyperbolic partial differential equations[END_REF] involves automatically taking smaller time steps where the mesh is refined. In their approach refined grids are laid over regions of the coarse mesh. Information is then passed between the grids by means of injection and interpolation. Flaherty et al. [START_REF] Flaherty | Adaptive local refinement with octree load-balancing for the parallel solution of three-dimensional conservation laws[END_REF] developed a parallel, adaptive discontinuous Galerkin method with a local forward Euler scheme which relies on interpolating values in time at interfaces between time steps of different sizes. This scheme, however, does not appear to conserve flux along these interfaces. Also, only first order in time methods are discussed. Note that for these different algorithms, local time steps are usually selected to be fractions of the global time step. Recently a new time integration approach has been applied to equations of non-linear elastodynamics [START_REF] Lew | Asynchronous variational integrators[END_REF]. It is based on a discrete spacetime form of Hamilton's variational principle. This algorithm permits the selection of independent time steps in each mesh element. However, this approach is only applicable to Hamiltonian systems. In [START_REF] Omelchenko | Self-adaptive time integration of flux-conservative equations with sources[END_REF] Omelchenko and Karimabadi presented an asynchronous approach to a diffusion-advection-reaction equation in one dimension. Their method is based on discrete-event simulation (DES). They defined the concept of "flux capacitor", a variable in which they stored the flux between two cells until a later event. Mass conservation is guaranteed when the neighboring cell is updated and the flux capacitors are emptied. In order to improve the accuracy of their asynchronous method, Omelchenko and Karimabadi developed in [START_REF] Omelchenko | A time-accurate explicit multi-scale technique for gas dynamics[END_REF] a second order DES algorithm and showed that, at least numerically, the second order is indeed attained on one-dimensional gas dynamics test problems. More recently, they extended in [START_REF] Omelchenko | HYPERS: a unidimensional asynchronous framework for multiscale hybrid simulations[END_REF] their DES algorithm to multiple dimensions in the case of "logically uniform meshes". In [START_REF] Semiletov | CABARET scheme with conservation-flux asynchronous time-stepping for nonlinear aeroacoustics problems[END_REF], V.A. Semiletov and S.A. Karabasov proposed an asynchronous time-stepping algorithm for the Compact Accurately Boundary Adjusting high-REsolution Technique (CABARET), which is an Euler method for non-linear aeroacoustic problems. Numerical tests up to 3D show that the asynchronous algorithm maintains the same second-order convergence rate as the original single-time stepping scheme and that it also decreases the absolute numerical error. However, the question of the performance in terms of computational time remains open. In [START_REF] Markesteijn | Time asynchronous relative dimension in space method for multi-scale problems in fluid dynamics[END_REF] V.A. Semiletov and S.A. Karabasov adapted their asynchronous method for solving fluid dynamics equations. Their approach is based on a transformation of the governing equations in space and time in order to rewrite the initial problem on a uniform Cartesian grid, with adapted conditions at the grid interfaces. As for the numerical implementation of their scheme, the authors combined their new asynchronous method with their CABARET scheme presented in [START_REF] Semiletov | CABARET scheme with conservation-flux asynchronous time-stepping for nonlinear aeroacoustics problems[END_REF] which results in a quasi-second order scheme. The corresponding numerical tests were limited to two-dimensional problems.

In this paper we focus on the asynchronous method proposed by one of the co-author T. Unfer. His method uses local stability criteria and can be used on an arbitrary mesh. In [START_REF] Unfer | An asynchronous scheme with local time-stepping for multi-scale transport problems: application to gas discharges[END_REF], he developed an upwind asynchronous forward Euler scheme for the transport equation in one dimension. We studied the extension of this method to an arbitrary space dimension and we proved through a thorough analysis that the asynchronous scheme is first order convergent. In addition, we noticed that the asynchronous scheme reduces numerical diffusion and CPU time in comparison to the classical first order scheme. To improve the convergence rate of the asynchronous scheme, T. Unfer proposed in [START_REF] Unfer | An asynchronous framework for the simulation of the plasma/flow interaction[END_REF] an extension of the asynchronous integration procedure to higher order. This method was explored in some details from a practical viewpoint. We then explored it from a more theoretical viewpoint and we proved that the asynchronous scheme as it was presented in [START_REF] Unfer | An asynchronous framework for the simulation of the plasma/flow interaction[END_REF] is at most of order one. The derivation of higher order explicit LTS methods is a serious problem and only a few methods are available. In [START_REF] Taube | A high-order discontinuous Galerkin method with time-accurate local time stepping for the Maxwell equations[END_REF], A. Taube et al. proposed an ADER-DG scheme with time-accurate LTS for the Maxwell equations. They found via numerical study that their approach maintains the accuracy of the global time-stepping ADER-DG scheme "with a very small increase in computational effort". Grote et al. [START_REF] Grote | Runge-Kutta based explicit local time-stepping methods for wave propagation[END_REF] derived an explicit LTS scheme from standard explicit Runge-Kutta (RK) methods. They have proved that their scheme preserves the accuracy of the original RK scheme. However, their approach is applicable only in the case of two distinct regions in the mesh: a "coarse" region with the larger elements and a "fine" region with the smaller elements. In the fine region, the time step must be a fraction of the coarse time step. Here, we develop an asynchronous RK2 scheme in an arbitrary mesh.

This paper is outlined as follows. In section 2, we summarize basic ideas of the asynchronous methodology. Then, we consider an arbitrary space discretization of a given linear partial differential equation and we apply the ARK2 method. In section 3, we prove that the asynchronous scheme is second order convergent. Next, numerical experiments that illustrate the convergence rate of the asynchronous scheme and validate the gain in computation time are presented in section 4. Finally, some conclusions and final remarks are given in section 5.

Numerical asynchronous simulation

For stability reasons, any explicit method has to fulfill a time-step restriction. This may reduce the method efficiency, because automatic grid generation tools for unstructured meshes may produce very small cells in geometrically complicated parts of the computational domain. To overcome this problem, we propose an asynchronous scheme with local time steps.

We consider

a polygonal domain in R d (d ≥ 1). Let T = {K i , i = 1, . . . , N} be a partition of the domain in N polyhedral volumes K i . For all i ∈ J1, NK, N (i) = { j, |K i ∩ K j | = 0} is the set of indices of neighboring volumes of K i ; |K i ∩ K j | denotes the (d -1) positive measure of K i ∩ K j . Let N + (i) (resp. N -(i))
be the set of the indices of the neighbors of the element K i who receive (resp. send) the outflow (resp. inflow) fluxes from (resp. to) K i .

We suppose an arbitrary space discretization of a given linear partial differential equation without a source term, for the sake of simplicity. This leads to a system of coupled ordinary differential equations

dy i dt (t) = B i y(t), ∀K i ∈ T , ( 1
)
where y i and y are respectively vectors that contain the degrees of freedom (dofs) in the element K i and in all the elements, and B i a matrix coming from the space discretization.

Asynchronous methodology

The asynchronous methodology is based on two key ideas. First of all, the asynchronous algorithm permits the selection of independent time steps in each mesh element. Then, the sequence of update times of the cells is no longer defined by continuously adding a uniform time step as in the classical case. We need to introduce the sequence of the most urgent refresh time tag (t i ) i≥1 defined by

t 1 = min j ( t j ) := t i , t 2 = min min j =i ( t j ), 2 t i , . . .
where for all i, the time step t i of the element K i is determined from the local stability restriction. Let (t r i ) r be the sequence of update times of a given element K i which is defined by t r i = r t i , ∀r ≥ 0.

For each element, the source term depends only on the current cell whereas the fluxes depend also on neighbors. Hence the second basic idea consists in separating the source term and the fluxes. In other words, we need to split the term B i (1) in three contributions: V i for a local volume contribution, F + i, j for a local flux contribution and F - i, j for a flux contribution from the neighbor. Then (1) becomes

dy i dt (t) = V i y i (t) + j∈N (i) F + i, j y i (t) + F - i, j y j (t), ∀K i ∈ T .
(

) 2
The concept of the asynchronous time integration can be summarized in three main phases. First, at the start-up time all densities and fluxes are initialized. Then, a second phase is carried out by continuously applying the following steps until the global simulation clock is advanced past the simulation finish time; see also Fig. 1: find the most urgent cell to be refreshed by using the CFL condition, compute the values of the densities which are needed to compute the fluxes, compute the next refresh time of the current cell and then update the values of the fluxes. Finally, the third phase is to build the solution at the output time.

The number of updated cells in the second step in Fig. 1 depends on the treated system and on the space discretization method. For example, in the case of an upwind scheme, we need to update only the outflow fluxes and then only the current cell K i and the neighbors K j , ∀ j ∈ N + (i) are involved; see also [START_REF] Toumi | Asynchronous numerical scheme for modeling hyperbolic systems[END_REF]. Whereas for a centered scheme, we need at least to update all the neighbors. In all cases, only few cells are involved.

Note that the critical point for speeding up an asynchronous method in terms of computation time is the searching algorithm for the most urgent time tag to be treated, some alternative options are presented in [START_REF] Unfer | An asynchronous scheme with local time-stepping for multi-scale transport problems: application to gas discharges[END_REF].

ARK2 algorithm

To present the asynchronous Runge-Kutta 2 algorithm, we start with the semi-discretized system written in the following form as explained in section 2

dy i dt (t) = V i y i (t) + j∈N (i) F + i, j y i (t) + F - i, j y j (t), ∀K i ∈ T .
(3) 2. For all the neighbors of the element i, i.e. j ∈ N (i), we update the tentative values

ỹ j (t cur ) = ỹ j (t pre, j ) + (t cur -t pre, j ) ỹ t j (t mid, j ),
t pre, j is the last time where the tentative dofs has been updated and t mid, j is the "mid-point" time for the element j.

3. Affect the dofs and update the last refresh times

y p i = ỹi (t cur ), t pre,i = t cur and t pre, j = t cur , ∀ j ∈ N (i).
4. After all elements K i such that t cur / t i is an integer have been updated, update their slopes

y t p i = V i y p i + j∈N (i) F + i, j y p i + F - i, j y n j j + (p t i -n j t j ) y t n j j ,
where n j = t cur t j . The final simulation time should correspond to a common meeting time of all cells i.e. a multiple of the lowest common multiple (LCM) of all the local time steps. If it is not the case, we use an interpolation at least of order 2.

Properties of the asynchronous scheme

Convergence of the ARK2 scheme

The asynchronous RK2 scheme is derived from the standard explicit Runge-Kutta 2 method. We shall now prove that the ARK2 scheme is second order convergent.

Notations

1. Let us begin with some notations concerning the temporal evolution of the system between two successive update times of an element K i which are multiples of t i ; see also Fig. 2 (left) (a) Let t r i +1 i := (r i + 1) t i = t cur ≤ T where T is the final time of the simulation. Suppose that the element K i has been updated m times between t r i i and t

r i +1 i at the instants: t 1 pre,i , t 2 pre,i , ..., t m pre,i , 0 ≤ m ≤ t i min j∈N (i) ( t j ) + 1. (b) Let δ k i = t k pre,i -t k-1 pre,i , for all 2 ≤ k ≤ m, δ 1 i = t 1 pre,i -t r i i , δ m+1 i = t r i +1 i -t m
pre,i and

m+1 k=1 δ k i = t i . Then, t k mid,i = t k pre,i -δ k i /2 for 1 ≤ k ≤ m and t m+1 mid,i = t r i +1 i -δ m+1 i /2. (c)
Suppose that for all j ∈ N (i) we have the following assumptions: i. between t r i i and t r i +1 i , the element j has been updated z r i j times at instants which are multiples of t j then 0 ≤ z r i j ≤ m. For example, in Fig. 2 (left), z r i n = 1 and z r i j = 5 where n and j are two neighbors of i. ii. r j t j ≤ r i t i and at the instant t k pre,i , r k j t j denote the last update time of the element j which is multiple of t j just before t k pre,i . Note that there may be situations where r k+1 j t j = r k j t j . In Fig. 2 for instance, for the neighbor n, r 4

j t j = r 3 j t j = r 2 j t j = r 1 j t j = r n t n . iii. t k mid,i, j = t k next,i, j - t k next,i, j -t k pre,i, j 2 = t k next,i, j + t k pre,i, j 2 , with: A. t k pre,i, j = max(r k j t j , r i t i ); B. t k next,i, j = min((r k j + 1) t j , (r i + 1) t i ).
2. The update of an element K i involves its neighbors. The neighbors of K i depend on their neighbors, the neighbors of the neighbors of K i involves their neighbors and so on. Then, we introduce the following notations; see also Fig. 2 (right).

(a) The elements K j , j ∈ N (i), are the neighbors of K i . (b) The elements K {i,1} , {i, 1} ∈ N ( j), are the neighbors of K j and then the neighbors of the neighbors of K i . (c) The elements K {i,2} , {i, 2} ∈ N ({i, 1}), are the neighbors of the neighbors of the neighbors of

K i . (d) The elements K {i,n} , {i, n} ∈ N ({i, n -1}), are the neighbors of the elements K {i,n-1} .
Let h r i {i,t s } be the number of dependencies presented by the element K i between t r i i and t s , where t s is the time of the last synchronization of the solution which can particularly be the initial time t 0 . In other words, the elements K {i,h r i i,ts } use their values computed at the time of the last synchronization of the solution. Note that the study of the dependency will be necessary for the theoretical study of the ARK2 scheme which needs the control of the temporal evolution of the system. However, for the asynchronous algorithm, the update of a given element involves only the neighbors and then no need of the neighbors of the neighbors. This dependency is implicit and it is hidden behind the intermediate steps of the algorithm. 3. Let us introduce some elements for the identification:

(Bx) i = V i x i + j∈N + (i) F + i, j x i + j∈N -(i) F - i, j x j (B P i x) i = V i x i + j∈N + (i) F + i, j x i
where P i is a diagonal projection matrix on the dof of the element i.

Proposition 1. The evolution of any given cell i between the two instants t

r i i = r i t i and t r i +1 i = t r i i + t i is given by the following expression y r i +1 i = y r i i + t i V i y r i i + t i 2 y t r i i + m+1 k=1 δ k i j∈N + (i) F + i, j y r i i + t k mid,i, j -r i t i y t r i i + m+1 k=1 δ k i j∈N -(i) F - i, j y r k j j + t k mid,i, j -r k j t j y t r k j j .
Proof. By definition, we have

y r i +1 i = ỹt m pre,i i + δ m+1 i ỹ t i (t m+1 mid,i ) = ỹt m-1 pre,i i + δ m i ỹ t i (t m mid,i ) + δ m+1 i ỹ t i (t m+1 mid,i ) . . . = y r i i + m+1 k=1 δ k i ỹ t i (t k mid,i )
where for all 1

≤ k ≤ m + 1 ỹ t i (t k mid,i ) = V i y r i i + t i 2 y t r i i + j∈N + (i) F + i, j y r i i + t k mid,i, j -r i t i y t r i i + j∈N -(i) F - i, j y r k j j + t k mid,i, j -r k j t j y t r k j j
To conclude, we just use that

δ k i × t k mid,i, j = (2r i + 1)( t i ) 2 2 , ∀ j ∈ N (i). ( 4
)
In particular, between the two instants t s , any time of the synchronization of the solution, and t s + t i we have ∀ j ∈ N (i)

m+1 k=1 δ k i × t k mid,i, j = ( t i ) 2 2
.

Proof. We first prove that, for a given element j ∈ N (i), between two successive times of update of the element K j which are multiples of t j and which are between t r i i and t r i +1 i , t k mid,i, j is constant. In other words, if r i t i ≤ (r j + 1) t j < (r j + l) t j ≤ (r i + l + 1) t i , l integer, then for all k such that (r j + l) t j < t k pre,i ≤ (r j + l + 1) t j we have

t k pre,i, j = max(r k j t j , r i t i ) = max((r j + l) t j , r i t i ) = (r j + l) t j and t k next,i, j = min((r k j + 1) t j , (r i + 1) t i ) = min((r j + l + 1) t j , (r i + 1) t i ) = (r j + l + 1) t j . As a result, t k mid,i, j = (2(r j + l) + 1) t j 2 . ( 5
)
Now, we prove Equation [START_REF] Grote | Runge-Kutta based explicit local time-stepping methods for wave propagation[END_REF]. Suppose that for all j ∈ N (i), z r i j is the number of times of update of the element K j at times that are multiples of t j between r i t i and (r i + 1) t i . There are three cases:

1. If z r i j = 0 then r j t j ≤ r i t i ≤ (r i + 1) t i ≤ (r j + 1) t j and for all 1 ≤ k ≤ m + 1, t k mid,i, j = (r i + 1) t i + r i t i 2 = (2r i + 1) t i 2 , therefore, m+1 k=1 δ k i × t k mid,i, j = (2r i + 1)( t i ) 2 2 . 2. If z r i j = 1 then we have m+1 k=1 δ k i × t k mid,i, j = ((r j + 1) t j + r i t i ) 2 × ((r j + 1) t j -r i t i ) + ((r i + 1) t i + (r j + 1) t j ) 2 × (r i + 1) t i -(r j + 1) t j = 1 2 (r i + 1) 2 t 2 i -r 2 i t 2 i = (2r i + 1)( t i ) 2 2 3. If z r i j ≥ 2 then r j t j ≤ r i t i < (r j + 1) t j < ... < (r j + z r i j ) t j ≤ (r i + 1) t i < (r j + z r i j + 1) t j . Using (5) we have the following result m+1 k=1 δ k i × t k mid,i, j = ((r j + 1) t j + r i t i ) 2 × ((r j + 1) t j -r i t i ) + z r i j -1 s=1 (2(r j + s) + 1) t j 2 × t j + ((r i + 1) t i + (r j + z r i j ) t j ) 2 × ((r i + 1) t i -(r j + z r i j ) t j ) = ((r j + 1) t j ) 2 2 - (r i t i ) 2 2 + ((r i + 1) t i ) 2 2 - ((r j + z r i j ) t j ) 2 2 + z r i j -1 l=1 (2(r j + l) + 1) × ( t j ) 2 2 .
Note that

z r i j -1 l=1 (2(r j + l) + 1
) is a sum of an arithmetic sequence with a common difference equal to two. Therefore,

z r i j -1 l=1 (2(r j + l) + 1) = z r i j -1 l=0 (2(r j + l) + 1) -(2r j + 1) = z r i j (4r j + 2z r i j ) 2 -(2r j + 1) = z r i j 2 + 2r j (z r i j -1) -1 and then m+1 k=1 δ k i × t k mid,i, j = (2r i + 1)( t i ) 2 2 .
The proof is then complete. 2

Proposition 3. For any element i and at each time t r i i = r i t i , r i an integer, we have the following expression of the slope i y t

r i i = B P i y r i i + j∈N -(i) F - i, j y r j j + δ i, j ⎛ ⎝ B P j y r j j + {i,1}∈N -( j) F - {i,1}, j y r {i,1} {i,1} ⎞ ⎠ + h r i i,ts -1 z=1 ⎛ ⎝ z n=1 δ {i,n-1},{i,n} {i,n}∈N -({i,n-1}) F - {i,n-1},{i,n} ⎞ ⎠ × ⎛ ⎝ B P {i,z} y r {i,z} {i,z} + {i,z+1}∈N -({i,z}) F - {i,z},{i,z+1} y r {i,z+1} {i,z+1} ⎞ ⎠ ,
where for all j ∈ N -(i), δ i, j = r i t ir j t j > 0 and r j t j is the last time, multiple of t j , of updating the element j just before t r i i .

For all

1 ≤ z < h r i i,t s -1,
F + i, j y r i i + j∈N -(i) F - i, j y r j j + δ i, j y t r j j = B P i y r i i + j∈N -(i)
F - i, j y r j j + δ i, j y t r j j , the elements j are the neighbors of the element i and for all j ∈ N -(i), δ i, j = r i t ir j t j > 0 where r j t j = t r j j < t r i i is the last time, multiple of t j , where the element j has been updated. Using the notations illustrated in Fig. 2 .

In the same way, we define for all 3 ≤ z < h r i i,t s -1, the slopes

y t r {i,z} {i,z}
and δ {i,z-1},{i,z} = r {i,z-1} t {i,z-1}r {i,z} t {i,z} > 0. In particular, for z = h r i i,t s -1 the slope is defined by Fig. 3. The temporal evolution of the system between two instants of synchronization t q and t q+1 = t q + t max ; example of five elements. B P {i,h

r i i -1} y r {i,h r i i -1}
{i,h

r i i -1} + {i,h r i i }∈N -({i,h r i i -1}) F - {i,h r i i -1},{i,h r i i } y 0 {i,h r i i } + δ {i,h r i i -1},{i,h r i i } (B y 0 ) {i,h r i i } and then y t r i i = B P i y r i i + j∈N -(i) F - i, j ⎛ ⎝ y r j j + δ i, j ⎛ ⎝ B P j y r j j + {i,1}∈N -( j) F - {i,1}, j y r {i,1} {i,1} ⎞ ⎠ ⎞ ⎠ + h r i i -1 z=1 ⎛ ⎝ z n=1 δ {i,n-1},{i,n} {i,n}∈N -({i,n-1}) F - {i,n-1},{i,n} ⎞ ⎠ × ⎛ ⎝ B P {i,z} y r {i,z} {i,z} + {i,z+1}∈N -({i,z}) F - {i,z},{i,z+1} y r {i,z+1} {i,z+1} ⎞ ⎠ , we supposed that K {i,0} ≡ K j , j ∈ N (i) are the neighbors of K i . 2
To study the convergence rate of the asynchronous scheme, we first need to suppose that the solution is synchronized at the instants t s which are multiples of a fixed time step t max . In the case of meshes where all the elements have a common meeting time, the solution is automatically synchronized at the instants which are multiples of the lowest common multiple (LCM) of all the local time steps and then t max is equal to the LCM. If it is not the case, so that the local time steps are completely independent, we have to suppose that the solution is regularly synchronized at the instants multiples of t max where, in this case, t max is computed in function of the time steps t i . For instance, we can choose t max = max i ( t i ). Note that the hypothesis of regular synchronization of the solution has no influence on the asynchronous aspect; see Fig. 3. In fact, the cells advance independently between t q and t q+1 . It is important to note that the regular synchronization of the solution was imposed only for demonstration purposes. Numerically, the synchronization is only made at the final time of the simulation independently of the mesh; see the ARK2 algorithm. Let (t q ) q≥0 be the sequence defined by t q+1 = t q + t max . The study of the convergence of the ARK2 scheme will be done between two instants of synchronization: t q and t q+1 . Let (t n q ) 1≤n≤N q be the sequence of the most urgent refresh time tags between t q and t q+1 , where N q = i s i and for all i, s i = t max t i

. (t k i,q ) k is the sequence of the temporal evolution of the element i between t q and t q+1 defined by t k i,q = t q + k t i (Fig. 3).

Theorem 1. For all the elements i and at any instant t r i i,q = t q + r i t i such that t q < t r i i,q ≤ t q+1 , the following expression holds true

y t r i i,q i = y t q i + (r i t i )(B y t q ) i + (r i t i ) 2 2 (B 2 y t q ) i (6) + 2r i n=3 ( t i ) n α n,r i (B P i ) n-2 (B 2 y t q ) i + R r i i,asyn , (7)
where:

α n,0 = 0 ∀n ≥ 3, α 3,r i = α 3,r i -1 + (r i -1) 2 +(r i -1) 2 , α 4,r i = α 3,r i -1 + α 4,r i -1 + (r i -1) 2 2 , α n,r i = α n,r i -1 + α n-1,r i -1 + α n-2,r i -1 2 ∀5 ≤ i ≤ 2(r i -1), α 2r i -1,r i = α 2(r i -1),r i -1 + α 2(r i -1)-1,r i -1 2 , α 2r i ,r i = α 2(r i -1),r i -1 2 , and R r i i,asyn ∞ ≤ M r i i n=3 C n,r i ( t i ) n B n ∞ y t q ∞ ( 8
) M r i i = max(2r i , h r i i ) and for all 3 ≤ n ≤ M r i i , C n,r i is a bounded constant.
Proof. The proof of Theorem 1 is very technical and many formulas are long, it is then placed in the appendix. 2

Theorem 2. The asynchronous Runge-Kutta 2 scheme is a second order accurate scheme.

Proof. Let y i (t q+1 ) be the semi-discrete exact solution of the ODE (3) at t q+1 , y t q+1 i,syn is the solution of the classical global time step Runge-Kutta 2 method at t q+1 and y t q+1 i the asynchronous solution. Then,

y i (t q+1 ) -y t q+1 i ∞ ≤ y i (t q+1 ) -y t q+1 i,syn ∞ + y t q+1 i,syn -y t q+1 i ∞ (9)
The classical global time step RK2 scheme is second order convergent then,

y i (t q+1 ) -y t q+1 i,syn ∞ ≤ C 0 ( t max ) 3 ,
where, for t max sufficiently small, C 0 is a constant independent of t max . To complete the proof we thus need to treat the second term on the right of (9). First, using Theorem 1, we get for all i y t s i q,i i

= y t q i + (s i t i )(B y t q ) i + (s i t i ) 2 2 (B 2 y t q ) i + 2s i n=3 ( t i ) n α n,s i (B P i ) n-2 (B 2 y t q ) i + R t s i q,i i,asyn
where for all i, t

s i q,i = t q + s i t i , s i = E t max t i
. As a first step, suppose that t max t i is an integer. As a result,

y t q+1 i ≡ y t s i q,i i = y t q i + ( t max )(B y t q ) i + ( t max ) 2 2 (B 2 y t q ) i + 2s i n=3 ( t i ) n α n,s i (B P i ) n-2 (B 2 y t q ) i + R t q+1 i,asyn
Consequently,

y t q+1 i,syn -y t q+1 i ∞ = E t q+1 i ∞ ,
where

E t q+1 i = 2s i n=3 ( t i ) n α n,s i B n-2 (B 2 y t q ) i - 2s i n=3 ( t i ) n α n,s i (B P i ) n-2 (B 2 y t q ) i -R t q+1 i,asyn = 2s i n=3 ( t i ) n α n,s i j∈N -(i) F i, j n-2 (B 2 y t q ) j -R t q+1 i,asyn , B is linear. Therefore, E t q+1 i ∞ ≤ C 1 ( t max ) 3 , ( 10
)
where for t max sufficiently small, C 1 is a constant independent of t max . Finally, note that for the case where t max t i

is not an integer we just need to synchronize the solution at t q+1 using a second order interpolation as it was explained when the ARK2 algorithm has been presented. 2

3.2. Comparison with the LTS-RK2 scheme presented in [START_REF] Grote | Runge-Kutta based explicit local time-stepping methods for wave propagation[END_REF] In [START_REF] Grote | Runge-Kutta based explicit local time-stepping methods for wave propagation[END_REF], Grote et al. derived an explicit LTS-RK scheme from standard explicit Runge-Kutta (RK) methods in the case of two distinct regions in the mesh: a "coarse" region with the larger elements and a "fine" region with the smaller elements. In the fine region, the time step must be a fraction of the coarse time step. They use the classical RK method of a given order inside the fine region and quadrature formula inside the coarse region. The intermediate values needed at the coarse and fine mesh interface are obtained through a judicious combination of interpolation and Taylor expansion. They have proved that their scheme preserves the accuracy of the original RK scheme. We started from their LTS-RK scheme and we adapted the algorithm presented in [START_REF] Grote | Runge-Kutta based explicit local time-stepping methods for wave propagation[END_REF] to the asynchronous aspect in order to apply the three main phases of the asynchronous integration. In particular, each cell is treated independently. In this paper, we only focus on order two. In the case of the mesh imposed in [START_REF] Grote | Runge-Kutta based explicit local time-stepping methods for wave propagation[END_REF], we proved that the ARK2 degenerates into the LTS-RK2 scheme. More precisely, using the asynchronous algorithm, we get the same result at each global time step t, where t denote the time step dictated by the CFL-condition in the coarser part of the mesh. Whereas, at each time step of size n t/p where t/p is the local time step inside the refined region of the mesh with p and n are integers and 0 < n < p, the result is different. In fact, we have the following result Proposition 4. In the case of two distinct regions in the mesh where the time step of the fine region is a fraction of the time step of the coarse region, the ARK2 scheme degenerates into the LTS-RK2 scheme presented in [START_REF] Grote | Runge-Kutta based explicit local time-stepping methods for wave propagation[END_REF].

Proof. Resuming the notation used in [START_REF] Grote | Runge-Kutta based explicit local time-stepping methods for wave propagation[END_REF]: = c ∪ f where c is the coarse region with the large time step t c and f is the fine region with the small time step t f = t c p , p is an integer. Introduce some elements for the identification

(Bx) i =V i x i + j∈N + (i) F + i, j x i + j∈N -(i) F - i, j x j (B P x) i =V i x i + j∈N + (i) F + i, j x i + j∈N - f (i) F - i, j x j (B(I -P )x) i = j∈N - c (i) F - i, j x j
where N c (i) and N f (i) are the neighbors of the element i which are respectively in the coarse mesh and the fine mesh.

In other words, j ∈ N c (i) then j ∈ N (i) ∩ c and j ∈ N f (i) then j ∈ N (i) ∩ f . P is a diagonal matrix, its diagonal entries are equal to zero or one. It identify the unknowns associated with the locally refined region. (I -P ) identify the unknowns associated with the coarse region. The previous results are proved for any given of mesh, in particular for a mesh with two regions. In fact, using the previous notations we have t max = t c . This is a natural choice because the solution is regularly synchronized at the instants multiples of t c . Then the sequence (t q ) q is defined by t q+1 = t q + t c , for all i ∈ f , (t k i,q ) 0≤k≤p is defined by t k i,q = k t f and for all j ∈ c , (t k j,q ) 0≤k≤1 is defined by t 0 j,q = t q and t 1 j,q = t q+1 . Let us use the previous results of Proposition 1 and Proposition 2 for a fine element, the proof for the coarse elements can be done in the same way. For all i ∈ f we have

y t r i +1 i,q i = y t r i i,q i + t f B P y t r i i,q i + ( t f ) 2 2 ⎛ ⎝ B P y t t r i i,q ⎞ ⎠ i + t f B(I -P )y t q i + (2r i + 1)( t f ) 2 2 B(I -P ) y t t q i , with y t t r i i,q i = B P y t r i i,q i + B(I -P ) y t q + r i t f y t t q i
Comparing with the expression of the slope given in Proposition 3, the dependencies with the neighbors are simplified. In fact, in this particular case, there are only two types of elements the fine cells which are grouped in the matrix P and the coarse cells which are grouped in the matrix (I -P ). All the fine elements has been updated in the same times and it is also the case for the coarse elements. Consequently, the expression of y i (6), i a fine element, in Theorem 1 is written in the following form

y t r i i,q i = y t q i + (r i t f )(B y t q ) i + (r i t f ) 2 2 (B 2 y t q ) i + 2r i m=3 ( t f ) m α m,r i (B P i ) m-2 (B 2 y t q ) i ,
and then we find the same expression given in [START_REF] Grote | Runge-Kutta based explicit local time-stepping methods for wave propagation[END_REF]. 2

Originality of the asynchronous formulation in ARK2

The originality of our approach with respect to the LTS scheme developed in [START_REF] Grote | Runge-Kutta based explicit local time-stepping methods for wave propagation[END_REF], and more generally with respect to any other standard local time-stepping approach, lies within its ability to adapt to any arbitrary mesh and its resilience for being second order convergent scheme independently of the chosen mesh. Indeed, for the standard LTS schemes, local time steps are usually selected to be fractions of the global time step whereas the asynchronous scheme permits the selection of independent time steps in each mesh element.

However, as mentioned in the introduction, the idea of using a different time-step in every element already exists in the literature. In [START_REF] Omelchenko | Self-adaptive time integration of flux-conservative equations with sources[END_REF], Omelchenko and Karimabadi introduced such an asynchronous scheme for drift-diffusion transport. The major difference between their asynchronous approach and the one presented here lies within the way the fluxes are updated. In [START_REF] Omelchenko | Self-adaptive time integration of flux-conservative equations with sources[END_REF], the concept of "flux capacitor", a variable in which the flux between two cells is stored, is used. In our approach, desynchronized fluxes can be used which allows fluxes recomputations in order to ensure the second-order convergence.

In [START_REF] Taube | A high-order discontinuous Galerkin method with time-accurate local time stepping for the Maxwell equations[END_REF], A. Taube et al. proposed an ADER-DG scheme with independent time steps for the Maxwell equations with arbitrary high order of accuracy in space and time, although no proof of the convergence rate is given. Their method consists in replacing the time derivatives with space derivatives. The computation of the flux between the grid cells is then based on the solution of generalized Riemann problems. Our asynchronous method is different, since it is derived from the standard explicit Runge-Kutta 2 method and thus the numerical solution needs to go through intermediate stages in order to be advanced in time.

Numerical tests

We consider the two-dimensional Maxwell equations in transverse electric mode:

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ ∂ ∂t E x = 1 ∂ ∂ y H z in × (0, T ) ∂ ∂t E y = -1 ∂ ∂x H z in × (0, T ) ∂ ∂t H z = 1 μ ∂ ∂ y E x -1 μ ∂ ∂x E y in × (0, T ) (11)
where E = E x , E y , 0

T and H = (0, 0, H z ) T are respectively the electric and magnetic fields. and μ are, respectively, the electric permittivity and magnetic permeability.

is the computational domain and T is the final time of the simulation. For the spatial discretization, we apply the second order P 1 -discontinuous Galerkin scheme with upwind fluxes.

Numerical convergence of the ARK2 scheme

We shall now present numerical experiments which validate the order of convergence of the ARK2 method derived in section 3. The one-dimensional tests assessing the second order convergence of the ARK2 scheme have already been published in [START_REF] Toumi | Asynchronous numerical scheme for modeling hyperbolic systems[END_REF]. In this paper, we consider the two-dimensional Maxwell equations [START_REF] Omelchenko | Self-adaptive time integration of flux-conservative equations with sources[END_REF]. The initial condition:

(E x ) 0 = (E y ) 0 = (H z ) 0 = cos(2π x),
is a rectangular domain of size [0, 1] × [0, 1] and the boundary conditions are supposed periodic.

First, we divide [0, 1] in X -direction, respectively in Y -direction, into tree equal parts. The left and right intervals are discretized with an equidistant mesh of size x, respectively y. The middle interval is discretized with an equidistant mesh of size p x/q, respectively p y/q where p and q are integers; see also Fig. 4 left. We start with a simple case, mesh 1, where the local time steps are fractions of the global time steps with p = 4 and q = 1. Then for a more general case, mesh 2, we suppose that p = 2, q = 3 and then p/q < 1. Finally, mesh 3 is to test the ARK2 algorithm in the case where the time steps are independent, we suppose that the positions x i and y i follow a polynomial low and thus and

x i = 4(1 -α) i -1 N x -1 - 1 2 3 + α i -1 N x -1 - 1 2 + 1 2 , ( 12
)
y i = 4(1 -α) i -1 N y -1 - 1 2 3 + α i -1 N y -1 - 1 2 + 1 2 , ( 13
)
where N x , respectively N y , is the number of cells in X -direction, respectively in Y -direction; see also Fig. 4 right. The parameter α define the slope at x = 1/2 and y = 1/2. For each mesh, each rectangle is splitted into two triangles.

For each mesh, we plot the infinity norm of the error versus various number of cells in X direction in logarithmic scale. Note that for the three treated meshes the number of cells in X direction and the number of cells in Y direction were chosen to be equal. We present a numerical study of the convergence rate of the asynchronous scheme in the case of mesh 1 in Fig. 5 left, in the case of mesh 2 in Fig. 5 right and in the case of mesh 3 in Fig. 6.

Effective gain in computational cost

The asynchronous method has been developed in order to reduce the computation cost in the case of multi-scale problems. To study the performance of our asynchronous method we consider the two-dimensional Maxwell equations [START_REF] Omelchenko | Self-adaptive time integration of flux-conservative equations with sources[END_REF] with null initial conditions and for the boundary conditions we inject by one of the sides of the computational domain , a plane wave with the equation:

u(x, t) = p 1 sin(2π × p 2 × (t -x/c)), (14)
where p 1 = 1 and p 2 = 10 9 are respectively the amplitude and the frequency of the wave and c = 3 × 10 8 m/s is the speed of light. λ = c × T is the wavelength and T = 10 -9 s is the final simulation time. is a rectangular domain of size [0, λ] × [0, λ]. To built a multi-scale mesh, we add a circle inside ; see Fig. 7 left. (the circle represents a two-dimensional projection of a cylindrical electric wire). The circle is of radius r = 0.001 × λ << λ. The multi-scale mesh produced by GMSH contains 8416 triangles with A max = 6.5 × 10 -5 and A min = 2.5 × 10 -9 where A max (resp. A min ) is the area of the The theoretical expected speed-up is an a priori estimate on the time reduction factor which depends only on the mesh. As expected, the real asynchronous speed-up is lower than the theoretical one mainly because of two reason. The first is that the asynchronous scheme access the data in memory in a more random way, yielding a higher number of "cache-miss" within the processor. The second is that the algorithm may require extra computations in order to ensure the second-order accuracy. This computational overhead depends on the implementation and on the test-case, nonetheless experiments conducted in both 2D and 3D on various test-case showed that the ratio between the real and the optimal speed-ups lies within a range from 1.5 to 2 [START_REF] Toumi | Methode numerique asynchrone pour la modelisation de phenomenes multi-echelles[END_REF]. This means that for application with a theoretical speed-up exceeding 3, the asynchronous integration practically provides an effective CPU speed-up.

Strategy for an asynchronous ARK2 parallel computation

The numerical example presented in the previous section shows the efficiency (in terms of computational time) of the asynchronous approach with respect to the global time stepping method. However, it has to be noted that it has been performed in a sequential approach. Conversely, given a sequential algorithm, many strategies for parallel computing have been developed in order to gain computational time. While the focus of this paper is on the theoretical order of accuracy for the ARK2 algorithm, we will present here a first insight on the possible strategies for a parallelization of the ARK2 algorithm, allowing for further computational time saving.

Deriving a parallel approach for the asynchronous integration method has first been studied in the context of a first-order Lax-Wendroff asynchronous time integration algorithm [START_REF] Guivarch | Parallelisation of MACOPA, a multi-physics asynchronous solver[END_REF]. In this section, we will focus on its application to the ARK2 method, following the openMP memory paradigm. Still, the ARK2 method is also compatible with distributed memory strategy introduced in [START_REF] Guivarch | Parallelisation of MACOPA, a multi-physics asynchronous solver[END_REF]. For shared memory, the key idea is to have as many DTS as openMP threads. The master thread is the master of time and is the only allowed to change the simulation time. All threads (master and slaves) have an action Thread safety is ensured since each action belongs to a single thread and is therefore rescheduled within its own DTS. Data safety is also ensured because the algorithm has been split in two parts. The first part consists in computing intermediate variables such as interface fluxes or source terms and generating an "integration list", which is the lists of all cells for which the time derivatives have to be recomputed at the current time. The second part consists in integrating the equations in the dedicated cells and recomputing the time derivatives using the conservation equation (flux balance and source term) according to the integration list. This list is then split within the threads. Hence, date safety is guaranteed provided that a single cell appears at most once in the integration list. Indeed, a thread may access data from a cell and from its neighbors but the writing occurs only within the considered cell.

The difference with [START_REF] Guivarch | Parallelisation of MACOPA, a multi-physics asynchronous solver[END_REF] lies within the fact that for first order time integration, the action list is treated then rescheduled within the DTS. For ARK2, the action list is treated once to compute the new mid-point slope but is not rescheduled and is treated a second time instead. An action can be treated the second time provided that for each cell, the mid-point slopes are already available from the neighbor cells (from previous simulation times or from another thread or from the thread itself).

A multi-threaded scaling speed-up test has been performed on a single node of the Eos supercomputer at Calmip [6].

The test case, which is strongly asynchronous, is the same as the one described on Fig. 4 (left), with a 101 × 101 mesh, generated using the polynomial law with α = 0.01. The speed-up obtained on one node, with thread numbers ranging from 1 to 8, is represented on Fig. 8 and shows that the openMP parallel version of the asynchronous algorithm has a good performance with respect to the ideal scaling. However, the used parallel implementation is only a first step and there is still room for improvements. A much more detailed study of the parallel paradigm for the ARK2 algorithm as well as in-depth performance analysis for both distributed and shared memory parallelization still remains to be made, this will be the subject of future works.

Conclusion

In the present paper, we derived an asynchronous RK2 scheme from the standard Runge-Kutta scheme which permits the selection of independent time steps in each mesh element. We presented the ARK2 algorithm to discretize a given linear partial differential equation with an arbitrary space discretization. At each time step, only one cell is involved and only the neighbors need to be updated, which makes the method easily applicable to any arbitrary mesh. We proved that the ARK2 scheme is a second order convergent scheme and thus it preserves the accuracy of the original RK2 scheme. Numerical experiments assess the convergence rate of the ARK2 scheme in one dimensional case, have already been published in [START_REF] Toumi | Asynchronous numerical scheme for modeling hyperbolic systems[END_REF]. In this paper we focused in two dimensional tests. The numerical study shows that the ARK2 scheme is second order convergent independently of the chosen mesh, which validates the theoretical results. Moreover, the numerical tests indicate that the use of the optimal CFL condition imposed on each cell independently of other cells, does not impact the accuracy and efficiency of the classical RK2 scheme. Furthermore, the results demonstrate that, comparing with the classical integration, the asynchronous scheme is effective in terms of computation time. In this paper, a simple academic two-dimensional test has been treated. More representative physical tests that validate the gain in computation time may be found in [START_REF] Toumi | Methode numerique asynchrone pour la modelisation de phenomenes multi-echelles[END_REF]. To increase the speedup efficiency of the asynchronous integration, a parallelization strategy using shared memory paradigm has been explained.

For the sake of simplicity, we imposed the linearity of the treated partial differential equation. The asynchronous method is likely to prove useful for non-linear systems. In the present work, we only focused on order two, mainly for reasons of readability and clarity. Extension of the algorithm in the case of ARKp (p > 2) will be the subject of future works.

Appendix A. Proof of Theorem 1

Proof. The proof is done by using the mathematical induction. We suppose that the formula (6) is true for t n q < t q+1 and that t n+1 q = t r i +1 i,q . Then according to Proposition 1 we have

y t r i +1 i,q i = y t r i i,q i + m+1 k=1 δ k i V i ⎛ ⎝ y t r i i,q i + t i 2 y t t r i i,q i ⎞ ⎠ + m+1 k=1 δ k i j∈N + (i) F + i, j ⎛ ⎝ y t r i i,q i + t k mid,i, j -r i t i y t t r i i,q i ⎞ ⎠ + m+1 k=1 δ k i j∈N -(i) F - i, j ⎛ ⎝ y t r j j,q j + t k mid,i, j -r k j t j y t t r j j,q j ⎞ ⎠ y t r i +1 i,q i
is then given by:

y t r i +1 i,q i = y t r i i,q i + t i ⎛ ⎝ V i y t r i i,q i + j∈N + (i) F + i, j y t r i i,q i ⎞ ⎠ + ( t i ) 2 2 V i y t t r i i,q i + j∈N + (i) F + i, j m+1 k=1 δ k i × t k mid,i, j -r i ( t i ) 2 y t t r i i,q i + m+1 k=1 δ k i j∈N -(i) F - i, j ⎛ ⎝ y t r j j,q j + t k mid,i, j -r k j t j y t t r j j,q j ⎞ ⎠ .
According to Proposition 2, we have

m+1 k=1 δ k i × t k mid,i, j = (2r i +1)( t i ) 2 2 , ∀ j ∈ N (i). As a result, m+1 k=1 δ k i × t k mid,i, j -r i ( t i ) 2 = ( t i ) 2 2 , ∀ j ∈ N (i). Then, y t r i +1 i,q i = y t r i i,q i + t i B P i y t r i i,q i + ( t i ) 2 2 B P i y t t r i i,q i (15) + m+1 k=1 δ k i j∈N -(i) F - i, j ⎛ ⎝ y t r j j,q j + t k mid,i, j -r k j t j y t t r j j,q j ⎞ ⎠ . ( 16
)
Now we use the Proposition 3 and the induction hypothesis to get the following expression of the slope i y t ( t j ) n α n,r k j (B P j ) n-2 (B 2 y t q ) j + t k mid,i, jr k j t j (r k j t j )(B 2 y t q ) j + (r k j t j ) 2 2

t r i i,q i = (B y t q ) i + (r i t i )(B 2 y t q ) i + (r i t i ) 2 2 B P i (B 2 y t q ) i + 2r i n=3 ( t i ) n α n,r i (B P i ) n-1 (B 2 y t q ) i + Rt r i i,q i,asyn , where Rt r i i,q i,asyn = (S R) t r i i,q i + (R J) t r i i,q i + (R Z) t r i i,q i (S R
)B P j (B 2 y t q ) j + ( t j ) n α n,r k j (B P j ) n-1 (B 2 y t q ) j .

To conclude, it remains to be noted that R r i i,asyn and Rr i i,asyn are bounded using the induction hypothesis, in particular, the formula [START_REF] Hundsdorfer | Monotonicity conditions for multirate and partitioned explicit Runge-Kutta schemes[END_REF]. Then the formula (6) is verified at t

r i +1 i,q with R t r i +1 i,q i,asyn ∞ ≤ M r i +1 i n=3 C n,r i +1 ( t i ) n B n ∞ y 0 ∞ , M r i +1 i = max(2(r i + 1), h r i +1 i,q ). 2

Fig. 1 .

 1 Fig. 1. The four steps of the second phase of the asynchronous algorithm with five elements.

Fig. 2 .

 2 Fig. 2. Notation: the temporal evolution of the system between two instants t ri i and t ri +1 i where the elements j and n are the neighbors of i (left). Dependency of the element K i on the other elements of the mesh (right). (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

Fig. 4 .

 4 Fig. 4. The computational domain = [0, 1] × [0, 1] with the refined region in the middle (left) and with polynomial mesh in X and Y directions; example N x = N y = 10 and α = 0.01 (right). Each rectangle is then splitted into two triangles.

Fig. 5 .

 5 Fig. 5. Second order ARK2 scheme for mesh 1 and mesh 2.

Fig. 6 .

 6 Fig. 6. Second order ARK2 scheme for the polynomial mesh where the parameter α = 0.01.

Fig. 7 .

 7 Fig. 7. The computational domain = [0, λ] × [0, λ] with circle in the middle of radius r (left) and the multi-scale mesh of (right).

Fig. 8 .

 8 Fig. 8. The speed-up versus the number of threads. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.) list to execute. In this context, an action corresponds to the refreshment of one cell of the mesh (Steps 1 to 4 in the ARK2 algorithm, see section 2.2).Thread safety is ensured since each action belongs to a single thread and is therefore rescheduled within its own DTS. Data safety is also ensured because the algorithm has been split in two parts. The first part consists in computing intermediate variables such as interface fluxes or source terms and generating an "integration list", which is the lists of all cells for which the time derivatives have to be recomputed at the current time. The second part consists in integrating the equations in the dedicated cells and recomputing the time derivatives using the conservation equation (flux balance and source term) according to the integration list. This list is then split within the threads. Hence, date safety is guaranteed provided that a single cell appears at most once in the integration list. Indeed, a thread may access data from a cell and from its neighbors but the writing occurs only within the considered cell.The difference with[START_REF] Guivarch | Parallelisation of MACOPA, a multi-physics asynchronous solver[END_REF] lies within the fact that for first order time integration, the action list is treated then rescheduled within the DTS. For ARK2, the action list is treated once to compute the new mid-point slope but is not rescheduled and is treated a second time instead. An action can be treated the second time provided that for each cell, the mid-point slopes are already available from the neighbor cells (from previous simulation times or from another thread or from the thread itself).A multi-threaded scaling speed-up test has been performed on a single node of the Eos supercomputer at Calmip [6].

  Pre-computed update times correspond to multiples of t i in each K i . The time advance is performed using the Discrete Time Scheduler (DTS) algorithm as defined in[START_REF] Unfer | An asynchronous scheme with local time-stepping for multi-scale transport problems: application to gas discharges[END_REF], t cur is the current simulation time.At t cur such that t cur / t i is an integer p, we have to perform the following operations 1. Compute the tentative dofs at t cur ỹi (t cur ) = ỹi (t pre,i ) + (t curt pre,i ) ỹ

	We denote by y m i the approximation of y i (m t i ) and by		y t	i m	the approximation of	dy i
										t i	(t mid,i ),
	where t mid,i = t cur -	t cur -t pre,i 2		and t pre,i is the last time where the tentative dofs has been updated. The tentative
	mid-point slope is defined as follows
	ỹ t i	(t mid,i ) = V i y	p-1 i	+	t i 2		y t	p-1 i	+	j∈N + (i)	F + i, j y	p-1 i	+ t mid,i, j -(p -1) t i	y t	p-1 i
			+	j∈N -(i)	F -i, j y	n j j + t mid,i, j -n j t j	y t	n j j	,
	where n j =	t pre,i t j	, t mid,i, j = t next,i, j -	t next,i, j -t pre,i, j 2		with t pre,i, j = max(n j t j , (p -1) t i ) and t next,i, j =
	min((n j + 1) t j , p t i ).

dt (m t i ).

  r {i,z} t {i,z} < t

				r i
				r i i is given by the following formula
	y t	r i i	= V i y	r i i +

i and δ {i,z-1},{i,z} = r {i,z-1} t {i,z-1}r {i,z} t {i,z} > 0 where t s is the time of the last synchronization of the solution and h r i {i,t s } is the number of dependencies presented by the element K i between t r i i and t s . We also supposed that K {i,0} ≡ K j , j ∈ N (i) are the neighbors of K i .

Proof. The expression of the slope i at t j∈N + (i)

  is the last time, multiple of t {i,1} , where the element K {i,1} has been updated. Then δ j,{i,1} = r j t jr {i,1} t {i,1} > 0 and

											we have the
	following result
	y t	r j j	= B P j y r j	j +	{i,1}∈N -( j)	F -j,{i,1} y	r {i,1} {i,1} + δ j,{i,1}	y t	r {i,1} {i,1}
											r {i,1} {i,1} < t	r j
	y t	r {i,1} {i,1}	= B P {i,1} y r {i,1}	{i,1} +	{i,2}∈N -({i,1})	F -{i,1},{i,2} y	r {i,2} {i,2} + δ {i,1},{i,2}	y t	r {i,2} {i,2}	.
	Likewise,
	y t	r {i,2} {i,2}	= B P {i,2} y r {i,2}	{i,2} +	{i,3}∈N -({i,2})	F -{i,2},{i,3} y	r {i,3} {i,3} + δ {i,2},{i,3}	y t	r {i,3} {i,3}

, where for all j ∈ N (i) and for all {i, 1} ∈ N ( j), r {i,1} t {i,1} = t j

  ) n α n,r j (B P j ) n-2 (B 2 y t q ) j + δ i, j B P j (r j t j )(B y t q ) j + (r j t j ) 2 ) n α n,r j (B P j ) n-2 (B 2 y t q r {i,1} ) n α n,r {i,1} (B P {i,1} ) n-2 (B 2 y t q P {i,z} y t q {i,z} + (r {i,z} t r {i,z} )(By t q ) {i,z} + (r {i,z} t i ) 2 r {i,z} ) n α n,r {i,z} (B P {i,z} ) n-2 (B 2 y t q F - {i,z},{i,z+1} y t q {i,z+1} + (r {i,z+1} t r {i,z+1} )(B y t q ) {i,z+1} + (r {i,z+1} t i ) 2 r {i,z+1} ) n α n,r {i,z+1} (B P {i,z+1}) n-2 (B 2 y t q Which gives the following result ) n α n,r i (B P i ) n-2 (B 2 y t q ) n α n,r i (B P i ) n-2 (B 2 y t q ,asyn + t k mid,i, jr k j t j (B y t q ) j + (r k j t j )(B 2 y t q ) j (B 2 y t q ) i + (r i ( t i ) 2 )B P i (B y t q i ) n α n,r i+1 (B P i ) n-2 (B 2 y t q

	(R J) t i + R + (r k r i i,q r k j t j,q = j t j ) 2 j∈N -(i) 2 B P j (B 2 y t q F -i, j	(r j t j ) 2 2 ) j + 2r k j n=3 ( t j ) n α n,r k (B 2 y t q ) j 2r j n=3 ( t j 2 (B 2 y t q ) j j (B P j ) n-1 (B 2 y t q	) j +	r k j j,q j,asyn Rt	.
	2r j We use Proposition 2 to compute the following term
	+ (r i t i ) 2 n=3 ( t j ) j + {i,1}∈N -( j) 2 ) i + ( t i ) 2 F -{i,1}, j r {i,1} t r {i,1} (B y t q 2 B P i (B y t q ) i	) {i,1}
	+	(r {i,1} t i ) 2 2 + m+1 k=1 δ k i	(B 2 y t q j∈N -(i)	) {i,1} + F -i, j (r k 2r {i,1} j t j )(B y t q	) j + t k mid,i, j -r k j t j (B y t q	) {i,1} ) j =	((r i + 1) t i ) 2 2	(B 2 y t q	) i .
	and Therefore,
	(R Z) t i y t r i +1 r i i,q i,q i = y t q = h t r i i,q -1 i + ((r i + 1) t i )(B y t q ⎛ ⎝ z δ {i,n-1},{i,n}	) i +	F -{i,n-1},{i,n} ((r i + 1) t i ) 2 2 (B 2 y t q ) i + ⎞ ⎠ × 2(r i +1)	) i + R	r i +1 i,q i,asyn , t
	where, R	z=1 i,asyn = R t r i +1 i,q t r i i,q i,asyn + t i B P i R n=1	{i,n}∈N -({i,n-1}) 2 i,asyn + t r i i,q ( t i ) 2 2 B P i Rt r i i,q i,asyn + (B 2 y t q k m+1	) {i,z} + δ k i ⎛ ⎝ R r k 2r {i,z} j t j,q j,asyn + t k mid,i, j -r k j t j	r k j j,q j,asyn Rt	⎞ ⎠	) {i,z}
	+ + {i,z+1}∈N -({i,z}) m+1 k δ k i j∈N -(i)	F -i, j	(r k j t j ) 2 2	(B 2 y t q	) j +	2r k j n=3	2	(B 2 y t q	) {i,z+1}
				2r {i,z+1}
	+							) {i,z+1}
	(S R) t i r i t i,q	r i i,q	t is a function of the asynchronous rests, (R J)	r i i,q
	Now, we replace replace the values of y	y t i and y r i in the formula of y i r i r j	r i +1 i	given by the equation (15) and we use the induction hypothesis to
	y	t i	r i +1 i,q	= y t q i + (r i t i )(B y t q	) i +	(r i t i ) 2 2	(B 2 y t q	) i
	+	2r i	( t i ) i + R	t i,asyn + t i B P i y t q r i i,q i + (r i t i )(B y t q	) i
				n=3
	+	(r i t i ) 2 2	(B 2 y t q	) i +	n=3 2r i	( t i ) i + R	r i i,q i,asyn t
	+	) t i ( t i ) 2 r i i,q = B P i R 2 B P i (B y t q t r i i,q i,asyn + ) i + (r i t i )(B 2 y t q j∈N -(i) F -i, j R	t j,asyn + δ i, j B P j R r j j,q ) i + (r i t i ) 2 2 B P i (B 2 y t q t r j j,q j,asyn + ) i	{i,1}∈N -( j)	F -{i,1}, j R	t {i,1},asyn r {i,1} {i,1},q
	+ + 2r i h r i i,t q -1 z=1	⎛ ⎝ n=1 z	δ {i,n-1},{i,n}	{i,n}∈N -({i,n-1}) ) i + Rt r i i,q i,asyn + F -{i,n-1},{i,n} k=1 m+1 δ k i	j∈N -(i) ⎞ ⎠	F -i, j	y t q j
	× B P {i,z} R + (r k j t j )(B y t q t r {i,z} {i,z},q {i,z},asyn + ) j + (r k j t j ) 2 2	(B 2 y t q	F -{i,z},{i,z+1} R ) j + 2r k j	t {i,z+1},asyn . r {i,z+1} {i,z+1},q

{i,z+1}∈N -({i,z})

n=3 ( t B n=3 ( t n=3 ( t i

is the rest coming from the neighbors j where j ∈ N -(i) and (R Z) i is the rest coming from the dependency of the element i on the other elements of the mesh. j for all j ∈ N (i).

n=3 ( t i ) n α n,r i (B P i ) n-1 (B 2 y t q n=3 ( t j ) n α n,r k j (B P j ) n-2 (B 2 y t q ) j jn=3 ( t

✩ This work was funded by the French National Research Agency under contract No. ANR-11-MONU-0019.


Julie Quilliot

Marion Couderc

Charlotte Giraud

Martin Soubrier

Sylvain Mathieu

Efficacy of intra-articular hyaluronic acid injection in knee osteoarthritis in everyday life