Nicolas Crouseilles

Matthieu Kuhn

Guillaume Latu

Comparison of numerical solvers for anisotropic diffusion equations arising in plasma physics

Keywords: plasma physics, anisotropic diffusion, semi-implicit time schemes, high order methods physique des plasmas, diffusion anisotrope, schéma semi-implicite, méthodes d'ordre élevé

This work is devoted to the comparison of numerical schemes to approximate anisotropic diffusion problems arising in tokamak plasma physics. We focus on the spatial approximation by using finite volume method and on the time discretization. This latter point is delicate since the use of explicit integrators leads to a severe restriction on the time step. Then, implicit and semi-implicit schemes are coupled to finite volumes space discretization and are compared for some classical problems relevant for magnetically confined plasmas. It appears that the semi-implicit approaches (using ARK methods or directional splitting) turn out to be the most efficient on the numerical results, especially when nonlinear problems are studied on refined meshes, using high order methods in space.

Introduction

Diffusion problems occur in a large number of applications such that image processing [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF], resonance imaging [START_REF] Basser | Diffusion-tensor MRI: theory, experimental design and data analysis -a technical review[END_REF] or transport in geological formations [START_REF] Berkowitz | Characterizing flow and transport in fractured geological media: A review[END_REF]. In this work, we focus on magnetically confined plasmas applications for which anisotropy is generated by the magnetic field [START_REF] Braginskii | Transport Processes in a Plasma[END_REF]. In this context, anisotropic phenomena arise between parallel (to the magnetic field) direction and perpendicular direction. This leads to diffusion models in which the ratio between perpendicular χ ⊥ and parallel χ diffusion coefficients can be small (typically χ ⊥ /χ ∼ 10 -3 ). As a consequence, at least two kinds of problems occur in the numeral simulations. First, the time step has to follow the smallest scales, which makes the computational cost prohibitive. Second, numerical errors pollute the perpendicular diffusion or transport [START_REF] Gunter | Finite element and higher order difference formulations for modelling heat transport in magnetised plasmas[END_REF][START_REF] Gunter | Modelling of heat transport in magnetized plasmas using non-aligned coordinates[END_REF][START_REF] Sharma | A fast semi-implicit method for anisotropic diffusion[END_REF][START_REF] Sharma | Preserving monotonicity in anisotropic diffusion[END_REF], due to the fast parallel dynamics.

Several works have been proposed to improve the existent numerical schemes, essentially by means of two aspects. The first one consists in the time integration; indeed, due to the severe restriction on the time step, it seems useful to deal with implicit or semi-implicit methods allowing unconditional stability. For example in [START_REF] Sharma | A fast semi-implicit method for anisotropic diffusion[END_REF], a semi-implicit scheme based on a directional splitting is proposed. In [START_REF] Filbet | Numerical study of a nonlinear heat equation for plasma physics[END_REF], a penalisation technique coupled with a directional splitting enables to invert a simple one-dimensional Laplacian instead of the full nonlinear diffusion problem. In [START_REF] Mentrelli | Asymptotic-Preserving scheme for highly anisotropic non-linear diffusion equations[END_REF], the authors proposed an Asymptotic Preserving method to deal with the stiff nonlinear diffusion problem. In [START_REF] Gunter | Finite element and higher order difference formulations for modelling heat transport in magnetised plasmas[END_REF], a fully implicit scheme is coupled with a high order finite element solvers.

The second aspect concerns the numerical method for the spatial approximation. Most of the time, second order finite volume schemes are used. To capture strong gradients and to decrease the perpendicular diffusion, it seems natural to consider higher order spatial discretizations and observe the influence on the numerical solution, as shown in [START_REF] Gunter | Finite element and higher order difference formulations for modelling heat transport in magnetised plasmas[END_REF] where high order finite elements or finite differences are considered. But, several models involving anisotropic diffusion operator also involve transport terms for which finite volume method are well suited. Another important point concerns the maximum principle and the monotonicity preservation along the direction of diffusion. We refer to [START_REF] Sharma | A fast semi-implicit method for anisotropic diffusion[END_REF][START_REF] Sharma | Preserving monotonicity in anisotropic diffusion[END_REF][START_REF] Nordbotten | Monotonicity of control volume methods[END_REF] for more details.

Our main goal in this work is to combine, test, validate and compare different recent numerical methods for anisotropic diffusion type models. On the one hand, the spatial discretization is performed using a finite volume method of second and fourth orders (see [START_REF] Zhang | A fourth-order accurate finite-volume method with structured adaptive mesh refinement for solving the advection-diffusion equation[END_REF]). On the second hand, we focus on the time discretization: since a key limitation of all explicit methods is that the time step is limited by the usual severe CFL condition, it is desirable to use implicit methods. In this work, we try and compare some approaches which are CFL free: a full implicit method and some semi-implicit methods recently introduced in the literature (in [START_REF] Sharma | A fast semi-implicit method for anisotropic diffusion[END_REF] and in [START_REF] Filbet | Numerical study of a nonlinear heat equation for plasma physics[END_REF]).

More precisely, the different approaches cited above are restricted to first or second order in time. The use of high order time discretization has not been yet considered so often in this domain. To achieve this task, Additive Runge-Kutta methods is an attractive strategy since it enables to reach arbitrary high order. In this work, we combine the ideas of [START_REF] Zhang | A fourth-order accurate finite-volume method with structured adaptive mesh refinement for solving the advection-diffusion equation[END_REF] and [START_REF] Filbet | Numerical study of a nonlinear heat equation for plasma physics[END_REF] to derive a resulting method which is fourth order accurate (in space and time) with no stability constraint related to the time step. The main idea due to [START_REF] Filbet | Numerical study of a nonlinear heat equation for plasma physics[END_REF][START_REF] Filbet | A class of asymptotic preserving schemes for kinetic equations and related problems with stiff sources[END_REF] is to penalize the (non)linear diffusion term by a Laplacian. The Laplacian is then considered implicit whereas the rest (the (non)linear diffusion term minus the Laplacian) is considered explicit. At the end, the cost of the obtained scheme is reduced to the inversion at each time step of a Laplacian (for which computational cost and complexity of the solver are reduced) instead of a nonlinear diffusion term. This can decrease the computational cost to a factor 10 (see [START_REF] Filbet | Numerical study of a nonlinear heat equation for plasma physics[END_REF]). A first order (in time) version would be the IMEX scheme, but the generalization to higher order is furnished by Additive Runge-Kutta methods. Let us remark that Asymptotic Preserving schemes have been developed in this framework (see [START_REF] Lozinski | Highly anisotropic temperature balance equation and its asymptotic-preserving resolution[END_REF][START_REF] Mentrelli | Asymptotic-Preserving scheme for highly anisotropic non-linear diffusion equations[END_REF][START_REF] Narski | Asymptotic Preserving scheme for strongly anisotropic parabolic equations for arbitrary anisotropy direction[END_REF][START_REF] Chacón | An asymptotic-preserving semi-Lagrangian algorithm for the time-dependent anisotropic heat transport equation[END_REF][START_REF] Van Es | Finite-difference schemes for anisotropic diffusion[END_REF]), which leads to a uniform accuracy with respect to the degree of anisotropy. Even if these techniques ensure the very interesting uniform accuracy property, the price to pay is to solve large linear systems which can induce problems of memory consumption in high dimensions.

Considering the spatial discretization, a good framework is the finite volume methods as proposed in [START_REF] Zhang | A fourth-order accurate finite-volume method with structured adaptive mesh refinement for solving the advection-diffusion equation[END_REF]. In particular, when one wants for example to couple with transport, it seems natural to use finite volume methods for both diffusion and transport operators [START_REF] Zhang | A fourth-order accurate finite-volume method with structured adaptive mesh refinement for solving the advection-diffusion equation[END_REF]. Moreover, when non periodic boundary conditions are assumed, spectral methods fail whereas several techniques enable finite volumes to reach high order accuracy even at the boundary (see [START_REF] Zhang | A fourth-order accurate finite-volume method with structured adaptive mesh refinement for solving the advection-diffusion equation[END_REF][START_REF] Filbet | Numerical study of a nonlinear heat equation for plasma physics[END_REF]). Second order finite volume methods are obtained by using midpoint quadrature to evaluate the integrals over faces. Recently, higher order quadratures have been used to achieve fourth order accuracy. The extension to nonlinear transport problems has been proposed in [START_REF] Mccorquodale | A High-Order Finite-Volume Method for Conservation Laws on Locally Refined Grids[END_REF] in the hyperbolic framework and we adapt it here for nonlinear diffusion problems.

In this work, we review the methods introduced above and propose high order extensions of some of them. Moreover, in order to get some insight on the behavior of these methods for modeling two or three dimensional physical problems, we compare them for some classical problems relevant for tokamak edge plasmas (see [START_REF] Beyer | Nonlinear dynamics of transport barrier relaxations in tokamak edge plasmas[END_REF][START_REF] Fuhr | Evidence from numerical simulations of transport-barrier relaxation in tokamak edge plasmas in the presence of electromagnetic fluctuations[END_REF][START_REF] Tamain | TOKAM-3D: a 3D fluid code for transport and turbulence in the edge plasma of tokamaks[END_REF]). More precisely, numerical simulations through electromagnetic fluid turbulence code require high spatial and time accuracy to recover physical relevant phenomena such as transport barrier relaxation (see [START_REF] Beyer | Nonlinear dynamics of transport barrier relaxations in tokamak edge plasmas[END_REF][START_REF] Fuhr | Evidence from numerical simulations of transport-barrier relaxation in tokamak edge plasmas in the presence of electromagnetic fluctuations[END_REF]). Plus, explicit time Inria integrators have difficulties (due to their prohibitive computational cost) to describe such long time phenomena so that implicit or semi-implicit numerical schemes are needed, and a comparative study on the efficiency and the performances applied to simpler problem are interesting in themself.

The use of high order semi-implicit time integrator (ARK2 or ARK4) coupled with fourth order finite volume method for the spatial discretization lead to an accurate solver in terms of perpendicular diffusion in both linear and nonlinear case. We will discuss the "domain of validity" of the time step to recover the desired order. Indeed, even if arbitrary large time step can be used, when the time step is too large, the order of accuracy is lost. However, the use of high order method enables to improve this. Finally, in view of coupling with transport term to describe complex phenomena (see [START_REF] Beyer | Nonlinear dynamics of transport barrier relaxations in tokamak edge plasmas[END_REF][START_REF] Fuhr | Evidence from numerical simulations of transport-barrier relaxation in tokamak edge plasmas in the presence of electromagnetic fluctuations[END_REF]), efficiency is required; hence, fully implicit solvers turn out to be too costly and the solvers we propose in this work, by restricting to the solving of one or two-dimensional discrete Laplacians, offer an interesting alternative between accuracy and efficiency.

We will focus on the numerical simulation of the following anisotropic diffusion equation satisfied by the function T (x, y), with (x, y)

∈ [0, L x ] × [0, L y ], (L x , L y > 0) ∂ t T = -∇ • Q, Q = -B ∇T (1.1)
where ∇ = (∂ x , ∂ y ) denotes the spatial gradient, Q the flux and B denotes a symmetric semidefinite matrix.

The matrix can take several forms such as B = ( b ⊗ b); in this case, Q is the heat flux along magnetic field lines and b = b(x, y) = (b x , b y ) is the magnetic field unit vector. The diffusion matrix B can also depend on the unknown T leading to a nonlinear diffusion problem. This situation adds technical problems since a fully implicit scheme needs to deal with a costly fix point algorithm.

The rest of the work is organized as follows. First, we present different time integration techniques. Then, we focus on the spatial approximations by recalling the finite volume framework adapted to the targeted model. A linear stability analysis is presented in section 4. Finally, section 5 is dedicated to a set of numerical test cases.

Time discretization

In this section, we present and recall some semi-implicit time integrators which allows to get unconditionally stable results for the Equations (1.1).

2.1 Semi-implicit method SH for Sharma-Hammett (see [START_REF] Sharma | A fast semi-implicit method for anisotropic diffusion[END_REF])

We detail the diffusion equation satisfied by

T ∂T ∂t = ∇ • (B∇T ) = [∂ x (b xx ∂ x T ) + ∂ x (b xy ∂ y T ) + ∂ y (b xy ∂ x T ) + ∂ y (b yy ∂ y T )] ,
where B is given by

B = b xx b xy b xy b yy , with b xx b yy ≥ b 2 xy .
The time discretization is based on a directional splitting and is composed of the following two steps; first, the x direction is considered

T ⋆ -T n ∆t = ∂ x (b xx ∂ x T ⋆ ) + ∂ x (b xy ∂ y T n ),
which can be rewritten as, using the notation

D ij = -∂ i (b ij ∂ j ) T ⋆ = (1 + ∆tD xx ) -1 (1 -∆tD xy )T n .
The second step deals with the y direction:

T n+1 -T ⋆ ∆t = ∂ y (b xy ∂ x T ⋆ ) + ∂ y (b yy ∂ y T n+1 ),
which can be rewritten as

T n+1 = (1 + ∆tD yy ) -1 (1 -∆tD yx )T ⋆ .
We finally get T n+1 as a function of T n in the following compact form

T n+1 = (1 + ∆tD yy ) -1 (1 -∆tD yx )(1 + ∆tD xx ) -1 (1 -∆tD xy )T n .
One of the main advantage of this method is that it involves only one dimensional spatial operator to invert, which is very efficient from a computational point of view. However, its generalization to the three dimensional case and to the nonlinear case is not straightforward. Indeed, as mentioned in [START_REF] Sharma | A fast semi-implicit method for anisotropic diffusion[END_REF], the splitting needs a modification in the three dimensional case to ensure stability. Moreover, in the nonlinear case B = B(T ),

B(T ) = b xx (T ) b xy (T ) b xy (T ) b yy (T ) ,
a direct generalization would lead to consider fully implicit the diagonal terms which requires a fixed point algorithm. The first step would be

T ⋆ -T n ∆t = ∂ x (b xx (T ⋆ )∂ x T ⋆ ) + ∂ x (b xy (T n )∂ y T n ), which rewrites as F(T ⋆ ) = T n + ∆t∂ x (b xy (T n )∂ y T n ) with F(T ⋆ ) = T ⋆ -∆t∂ x (b xx (T ⋆ )∂ x T ⋆ ).
Solving this nonlinear equation requires a fixed point or Newton type algorithm which turns out to be very costly when high order methods in space are used on refined meshes. An issue would be to consider a semi-implicit version:

∂ x (b xx (T n )∂ x )T n+1 and ∂ y (b yy (T n )∂ y )T n+1
. This will be used in the nonlinear tests presented in the last section. Finally, its generalization to higher order is not straightforward and requires a large number of stages.

Additive Runge-Kutta (ARK) scheme

This time integration can be viewed as an adaptation of the penalisation method used in [START_REF] Filbet | Numerical study of a nonlinear heat equation for plasma physics[END_REF] and introduced for stiff kinetic equations in [START_REF] Filbet | A class of asymptotic preserving schemes for kinetic equations and related problems with stiff sources[END_REF]. Let us recall the strategy in our context. Starting from the continuous equation, ∂T ∂t = ∇ • (B∇T ) ,

Inria

we introduce the parameter λ ∈ R + such that the largest eigenvalue |µ| of the matrix B satisfies |µ| < λ. Hence we introduce the operator λ∆ and apply the following time discretization (see [START_REF] Filbet | Numerical study of a nonlinear heat equation for plasma physics[END_REF])

T n+1 -T n ∆t = ∇ • (B∇T n ) -λ∆T n + λ∆T n+1 .
This method is attractive since the operator to invert is cheap and independent on space and time, which can be a great advantage from a memory consumption point of view. Moreover, its generalization to higher order (in time) can be performed using Additive Runge Kutta method (ARK). Using the notation T (1) = T n , one can compute the subsequent stage values T (s) , s = 2, 3, 4, 5, 6 by solving

(I -∆tγλ∆) T (s) = T n + ∆t L, where L = s-1 j=1 a s,j [∇ • (B∇) -λ∆] T (s) ,
with γ = a s,s can be imposed to be constant for all stages. In particular, we select two methods introduced in [START_REF] Liu | Some new additive Runge Kutta methods and their applications[END_REF] (see also [START_REF] Kennedy | Additive Runge-Kutta schemes for convectiondiffusion-reaction equations[END_REF][START_REF] Zhang | A fourth-order accurate finite-volume method with structured adaptive mesh refinement for solving the advection-diffusion equation[END_REF]): second order method (ARK2) corresponds to ARK.2.A.1 (s = 2) whereas fourth order method (ARK4) corresponds to ARK4.A.1 (s = 6). For both methods, the diffusion term L (and consequently ∇ • (B∇T )) is considered explicitly whereas an implicit treatment of the term ∆ is performed. The main advantages is that the matrix to invert is simple and sparse. Other choices of ARK methods can be done, to get A-stability or L-stability.

Remark 2.1 As discussed in [START_REF] Filbet | A class of asymptotic preserving schemes for kinetic equations and related problems with stiff sources[END_REF][START_REF] Lemou | Implicit schemes for the Fokker-Planck-Landau equation[END_REF][START_REF] Filbet | Numerical study of a nonlinear heat equation for plasma physics[END_REF], the choice of the "penalisation" operator is large. One way to observe this is to go back to a fully implicit method

T n+1 -T n ∆t = ∇ • B(T n+1 )∇T n+1 .
A standard fix point algorithm would be, with the initial guess T n+1,0 = T n I -∆t∇ • B(T n+1,s )∇ T n+1,s+1 = T n , s ≥ 0, until a criteria |T n+1,s+1 -T n+1,s | < ε is satisfied for a given threshold ε. This algorithm can be very prohibitive and one wants to decrease the computational cost. Another strategy is to linearize the nonlinear term and approximate the "linearized" term

∇ • B(T n+1 )∇T n+1 ≃ ∇ • B(T n )∇T n+1 = ∇ • (B(T n )∇T n ) + ∇ • B(T n )∇(T n+1 -T n ) ≃ ∇ • (B(T n )∇T n ) + ∇ • λI ∇(T n+1 -T n ) = ∇ • (B(T n )∇T n ) + λ∆(T n+1 -T n ) = ∇ • ([B(T n ) -λI]∇T n ) + λ∆T n+1
Then, the strategy proposed here can be interpreted as only one iteration of a fixed point algorithm, in which the nonlinear diffusion matrix B(T n ) has been approximated in such a way that efficient computations can be achieved. Hence, we choose here a Laplacian but other choices such as a non constant diagonal diffusion matrix can be considered. A deeper study would be interesting to observe the influence of this approximation (see [START_REF] Lemou | Implicit schemes for the Fokker-Planck-Landau equation[END_REF]). One can also refer to [START_REF] Crouzeix | An implicit-Explicit multistep method of the approximation of parabolic equations[END_REF][START_REF] Ch | Linearly implicit time discretization of nonlinear parabolic equations[END_REF] for related studies.

Remark 2.2 Let us remark that the combination of the two last approaches (SH and ARK) does not lead to an unconditionally stable method. Indeed, when the matrix B is not diagonal, it is not possible to stabilize the extra diagonal terms with λ∂ xx .

Spatial discretization

Let us denote by x i+1/2 = x min + (i + 1/2)∆x with ∆x = (x max -x min )/N x and y j+1/2 = y min + (j + 1/2)∆y with ∆y = (y max -y min )/N y , N x , N y ∈ N.

We recall and present two kinds of space discretization, both based on finite volume. We denote by T i,j the average of the unknown on a control volume V i,j

T n i,j = 1 ∆x∆y Vi,j
T (t n , x, y)dxdy.

As usual, we integrate the equation on a control volume

V i,j = C x i × C y j := [x i-1/2 , x i+1/2 ] × [y j-1/2 , y j+1/2 ] to obtain ∂ t T i,j = 1 ∆x∆y Vi,j [∇ • (B∇T )]dxdy that is ∂ t T i,j = 1 ∆x∆y C y j [(b xx ∂ x T )(x i+1/2 , y) -(b xx ∂ x T )(x i-1/2 , y)]dy + 1 ∆x∆y C y j [(b xy ∂ y T )(x i+1/2 , y) -(b xy ∂ y T )(x i-1/2 , y)]dy + 1 ∆x∆y C x i [(b yy ∂ y T )(x, y j+1/2 ) -(b yy ∂ y T )(x, y j-1/2 )]dx + 1 ∆x∆y C x i [(b xy ∂ x T )(x, y j+1/2 ) -(b xy ∂ x T )(x, y j-1/2 )]dx (3.2)
Up to now, the relationship is exact. We have to approximate the different terms, i.e. provide expressions linking face averaged quantities to the cell averaged unknown T i,j . We recall some useful formulas in Appendix 7.1 (see also [START_REF] Zhang | A fourth-order accurate finite-volume method with structured adaptive mesh refinement for solving the advection-diffusion equation[END_REF]). Let us remark that the finite volume framework enables to include naturally advective terms or to consider non periodic boundary conditions (see [START_REF] Zhang | A fourth-order accurate finite-volume method with structured adaptive mesh refinement for solving the advection-diffusion equation[END_REF]). The first step consists in expressing the integral of the product as a function of a product of integrals, using Proposition 7.1 (Appendix 7.1 page 35), at the desired order.

Then, we have to manage to different types of terms. The first line of (3.2) leads to integrals of the form C y )dy (which are approximated following Proposition 7.2). At the end, (3.2) becomes a system of ODE satisfied by T i,j since all the right hand side has been expressed as a function of the unknown T i,j , i = 0, . . . , N x , j = 0, . . . , N y .

When nonlinear diffusion problems are considered, one needs to evaluate integrals of the form is a given function (linear diffusion problems), one can easily approximate this integral at the Inria desired order (by standard quadrature). However, for nonlinear problems, it is more delicate. Following the strategy used in [START_REF] Mccorquodale | A High-Order Finite-Volume Method for Conservation Laws on Locally Refined Grids[END_REF] for nonlinear hyperbolic problems, we propose in the following a strategy to approximate at the fourth order C y j b xx (x i+1/2 , y)dy from T i,j . Let us recall the strategy in our context.

First, from the equality (point to cell)

T i,j = T (x i , y j ) + ∆x 2 24 (∂ 2 x T )(x i , y j ) + ∆y 2 24 (∂ 2 y T )(x i , y j ) + O(∆x 4 ) + O(∆y 4 ),
we can reconstruct a fourth order accurate pointwise value T (x i , y j ) using second order finite difference of cell averaged values up to fourth order accuracy (cell to point)

T (x i , y j ) = T i,j - ∆x 2 24 T i+1,j -2T i,j + T i-1,j ∆x 2 - ∆y 2 24 T i,j+1 -2T i,j + T i,j-1 ∆y 2 +O(∆x 4 ) + O(∆y 4 ).
Then the diffusion coefficients can be evaluated at (x i , y j ) to get for example b xx (T )(x i , y j ).

Once the point values diffusion coefficients are known at each i, j, one can reconstruct with a high order accuracy the cell averaged value b xx (T ) i,j with (point to cell)

b xx (T ) i,j := 1 ∆x∆y C x i C y j b xx (T )(x, y)dxdy = b xx (T )(x i , y j ) + ∆x 2 24 b xx (T ) i+1,j -2b xx (T ) i,j + b xx (T ) i+1,j ∆x 2 + ∆y 2 24 b xx (T ) i,j+1 -2b xx (T ) i,j + b xx (T ) i,j-1 ∆y 2 +O(∆x 4 ) + O(∆y 4 ).
(T ) i,j = b xx (T )(x i , y j ) + O(∆x 2 ) + O(∆y 2 ).
Obviously, the initial and final conditions have to be considered as a cell averaged quantity to ensure the correct orders of accuracy. Then, for a given pointwise initial (or exact) solution T (x i , y j ), one has to transform it in a cell averaged quantity following

T i,j = T (x i , y j ) + ∆x 2 24 T (x i+1 , y j ) -2T (x i , y j ) + T (x i-1 , y j ) ∆x 2 + ∆y 2 24 T (x i , y j+1 ) -2T (x i , y j ) + T (x i , y j-1 ) ∆y 2 ,
which is exact up to terms of order 4.

Remark 3.1 The extension of the SH scheme to the fourth order in space breaks the locality of the y direction for the first step (and in the x direction for the second step), hence inducing a full two-dimensional system to invert.

Remark 3.2

The choice of the λ parameter for the nonlinear case needs some discussions since the spectral radius may increase during the simulation. We follow the strategy proposed in [START_REF] Filbet | Numerical study of a nonlinear heat equation for plasma physics[END_REF].

It consists in setting λ as two time the spectral radius of the diffusion matrix, and updating this value when needed ( i.e. if the spectral radius increases and becomes greater than λ).

Linear stability analysis

In this section, we perform the Von Neumann stability analysis of the numerical schemes derived above. To do that, we consider a single mode T (t, x, y) = r(t) exp(-i(k x x + k y y) where r(t) is the amplification factor in time, and k x , k y denote the wave numbers in the x and y directions. We consider the simple case for which the matrix B is constant and of the form

B = b xx b xy b xy b yy , with b xx b yy ≥ b 2 xy .
Then, the numerical scheme becomes, for the second order case

dT i,j dt = 1 ∆x 2 b xx [T i+1,j -2T i,j + T i-1,j ] + 2 ∆x∆y b xy [T i+1,j+1 -T i+1,j-1 -T i-1,j+1 + T i-1,j-1 ] + 1 ∆y 2 b yy [T i,j+1 -2T i,j + T i,j-1 ] ,
whereas in the four order case, we get

dT i,j dt = 1 24∆x 2 b xx [-T i+2,j + 16T i+1,j -30T i,j + 16T i-1,j -T i-2,j ] + 4 ∆x∆y b xy (D x D y T ) i,j + 1 24∆y 2 b yy [-T i,j+2 + 16T i,j+1 -30T i,j + 16T i,j-1 -T i,j-2 ] ,
where the notation (D x D y T ) i,j is defined in Equations (7.4) and (7.5). The computations are detailed in the second order space discretization whereas the fourth order case is left to Appendix 7.2 page 39.

Semi-implicit scheme (SH)

First we focus on the stability analysis of the SH scheme. The first step of the scheme writes

T ⋆ i,j = T n i,j + ∆t ∆x 2 b xx T ⋆ i+1,j -2T ⋆ i,j + T ⋆ i-1,j + ∆t ∆x∆y b xy T n i+1,j+1 -T n i+1,j-1 -T n i-1,j+1 + T n i-1,j-1 ,
whereas the second step is

T n+1 i,j = T ⋆ i,j + ∆t ∆y 2 b yy T n+1 i,j+1 -2T n+1 i,j + T n+1 i,j-1 + ∆t ∆x∆y b xy T ⋆ i+1,j+1 -T ⋆ i+1,j-1 -T ⋆ i-1,j+1 + T ⋆ i-1,j-1 .
Following the Von Neuman analysis, we inject a plane wave solution: T (t, x, y) = r(t) exp -i(kxx+kyy) to get the amplification factor for the first step

r 1 = 1 -∆t ∆x 2 b xy sin(k y ∆y) sin(k x ∆x) 1 + 4∆t ∆x 2 b xx sin 2 (k x ∆x/2)
.

Inria

In the same way, the amplification factor r 2 for the second step writes

r 2 = 1 -∆t ∆x 2 b xy sin(k y ∆y) sin(k x ∆x) 1 + 4∆t ∆y 2 b yy sin 2 (k y ∆y/2)
, and then the total amplification factor is r = r 1 r 2 . We introduce the CFL number ncfl

∆t = ncfl ∆x 2 4 ,
so that the amplification factor writes, with A x = sin(k x ∆x/2) and A y = sin(k y ∆y/2)

r := r(ncfl, k x ∆x, k y ∆y) = (1 -(ncfl/2)b xy A x A y cos(k x ∆x/2) cos(k y ∆y/2)) 2 (1 + (ncfl/2)b xx A 2 x )(1 + (ncfl/2)b yy A 2 y )
. As expected, we observe that it is always smaller that one, which illustrates the unconditional stability. Moreover, as ncfl decreases (which means that ∆t decreases), one gets closer to the exact amplification factor, plotted in Figure 2:

One can prove that |r| ≤ 1 |r| ≤ (1 + (ncfl/2)|b xy ||A x ||A y |) 2 (1 + (ncfl/2)b xx A 2 x )(1 + (ncfl/2)b yy A 2 y ) ≤ (1 + (ncfl/2) √ b xx b yy |A x ||A y |) 2 (1 + (ncfl/2)b xx A 2 x )(1 + (ncfl/2)b yy A 2 y ) = (1 + (ncfl 2 /4)b xx b yy A 2 x A 2 y + ncfl √ b xx b yy |A x ||A y |) (1 + (ncfl 2 /4)b xx b yy A 2 x A 2 y + (ncfl/2)(b xx A 2 x + b yy A 2 y )) ≤ (1 + (ncfl 2 /4)b xx b yy A 2 x A 2 y + (ncfl/2)(b xx A 2 x + b yy A 2 y )) (1 + (ncfl 2 /4)b xx b yy A 2 x A 2 y + (ncfl/2)(b xx A 2 x + b yy A 2 y )) = 1,
exp[-(b xx k 2 x ∆x 2 + b yy k 2 y ∆y 2 + 2b xy k x ∆xk y ∆y)ncfl/4].
However, when ncfl becomes larger, the damping rate is not sufficiently strong in the perpendicular direction (first and third quadrants). Finally, we plot on Figure 3 the coefficient factor using a fourth order discretization in space. One can observe the stability property also with the fourth order. In the first and third quadrants, the damping rate appears to be lower that in the second order case. However, on the diagonal (second and fourth quadrants), the amplification factor is closer to the analytical solution which is a good property of high order scheme.

Additive Runge Kutta

This subsection is dedicated to the Von Neumann stability analysis of the ARK scheme. We shall restrict to the first order. The analysis can be extended in a straightforward way to arbitrary ARK time integrators. Let us recall the numerical scheme

T n+1 = T n + ∆t∇ • (B∇T n ) -∆tλ∆T n + ∆tλ∆T n+1 ,

RR n°8560

or in a more detailed way

1 -∆tλ(∂ 2 x + ∂ 2 y ) T n+1 = 1 + ∆t(b xx ∂ 2 x + b yy ∂ 2 y + 2b xy ∂ 2 xy ) -∆tλ(∂ 2 x + ∂ 2 y ) T n .
The penalisation coefficient λ is set to be greater than the spectral radius of matrix B, which means that λ ≥ |X 1,2 | where X 1,2 are the two eigenvalues of B

X 1,2 = 1 2 (Tr(B) ± (b xx -b yy ) 2 + 4b 2 xy ).
We consider the second order spatial discretization coupled with a first order time discretization. The numerical scheme then reads

T n+1 i,j - ∆tλ ∆x 2 [T n+1 i+1,j + T n+1 i,j+1 -4T n+1 i,j + T n+1 i-1,j + T n+1 i,j-1 ] = T n i,j + ∆t(b xx -λ) ∆x 2 [T n i+1,j -2T n i,j + T n i-1,j ] + ∆t(b yy -λ) ∆x 2 [T n i,j+1 -2T n i,j + T n i,j-1 ] +2 ∆tb xy 4∆x 2 (T n i+1,j+1 -T n i+1,j-1 -T n i-1,j+1 + T n i-1,j-1 ).
Inserting the plane wave solution T (t, x, y) = r(t) exp -i(kxx+kyy) , we obtain the following amplification factor r = r 1 /r 2 with (denoting ncfl= (4∆t)/∆x 2 )

r 1 = 1 - 4∆t(b xx -λ) ∆x 2 sin 2 (k x ∆x/2) - 4∆t(b yy -λ) ∆x 2 sin 2 (k y ∆x/2)
-2 ∆tb xy ∆x 2 sin(k x ∆x) sin(k y ∆x),

r 2 = 1 + 4∆tλ ∆x 2 (sin 2 (k x ∆x/2) + sin 2 (k y ∆x/2)).
Introducing A x = sin(k x ∆x/2) and A y = sin(k y ∆y/2), we get

r 1 r 2 = 1 - ncfl (b xx A 2 x + b yy A 2 y -2b xy A x A y cos(k x ∆x/2) cos(k y ∆x/2)) 1 + ncfl λ(A 2 x + A 2 y )
.

We now need to check the inequalities -1 < r 1 /r 2 < 1. The inequality r 1 /r 2 < 1 is always ensured since the numerator is non negative:

b xx A 2 x + b yy A 2 y -2b xy A x A y cos(k x ∆x/2) cos(k y ∆x/2) ≥ b xx A 2 x + b yy A 2 y -2b xy A x A y , and
b xx A 2 x + b yy A 2 y -2b xy A x A y = 1 b yy b yy b xx A 2 x + b 2 yy A 2 y -2b yy b xy A x A y ≥ 1 b yy (b xy A x -b yy A y ) 2 .

Inria

The inequality -1 < r 1 /r 2 rewrites as

ncfl (b xx A 2 x + b yy A 2 y -2b xy A x A y cos(k x ∆x/2) cos(k y ∆x/2)) 1 + ncfl λ(A 2 x + A 2 y ) ≤ 2.
But, with the following estimate on the numerator

b xx A 2 x + b yy A 2 y -2b xy A x A y cos(k x ∆x/2) cos(k y ∆x/2)) ≤ b xx A 2 x + b yy A 2 y + 2 |b xx | |b yy |A x A y ≤ b xx A 2 x + b yy A 2 y + b xx b yy A 2 x A 2 y ≤ 2 b xx A 2 x + b yy A 2 y , we deduce ncfl (b xx A 2 x + b yy A 2 y ) 1 + ncfl λ(A 2 x + A 2 y ) ≤ 1, which is always true since λ = max(|X 1 |, |X 2 |) ≥ max(b xx , b yy ).
In Figures 4 and5, we plot the amplification factor for different ncfl for the second and the fourth order case. The computations for the fourth order case are detailed in Appendix 7.2 page 39. As previously, we choose b xx = b yy = b xy = 1/2. As expected, the amplification factor is always smaller than 1 which traduces the unconditional stability. Finally, in Figures 6, we plot the amplification factor corresponding to the approach combining the two methods (SH and ARK). One can observe, for the given example, that it is not always stable (only when ncfl≤ 1).

Test cases

We propose here several test cases to validate our numerical methods. These tests can be classified into two categories: analytical test cases and physical test cases. With the analytical ones, we intend to observe numerical errors magnitude and order, whereas the physical test cases are evaluated qualitatively. Those tests are also of increasing difficulty, from constant to non linear diffusion problem, including non constant in space (but constant in time) diffusion. We intend to observe here the benefits of order 4 in space methods and high order in time on the numerical perpendicular diffusion error. Also, we aim to highlight the advantages and the reduced costs of semi-implicit time schemes (e.g. SH, ARK2 or ARK4), observing the impact of the ncfl parameter on the error magnitude and order, to be balanced with the execution time gain. All the tests are 2D, but generalization to 3D test case is most of the time straightforward.

In this whole section, ∆t is the time step, ∆x is the spatial steps on a N 3 x spatial grid, and the time step is defined by: ∆t = ncfl∆x 2 /(4ρ(B)),

where B is the diffusion matrix, ρ(B) is the spectral radius of B and ncfl is a number which will be determined in each numerical test.

Analytic test case

We consider the following solution T (t, x, y) = e -10t sin(πx) cos(πy), x, y ∈ [-1, 1], t ≥ 0. Then, computing ∇ • (B∇T ) with leads to

B = 1 1 1 1 ,
∇ • (B∇T ) = ∂ 2 x T + 2∂ 2 x,y T + ∂ 2 y T = -2π 2 T -2π 2 e -10t cos(πx) sin(πy).
Hence, T (t, x, y) is a solution of

∂ t T = ∇ • (B∇T ) + Q, with Q(t, x, y) = -(10 + 2π 2 )T (t,
x, y) -2π 2 e -10t cos(πx) sin(πy).

For this equation, we compare the different methods presented before. For all these methods, the error between the analytical and numerical solution behaves as O(∆t p ) + O(∆x q ) where p is the order in time and q the order in space. According to the size of ∆t, one interesting point is to observe the total error. In particular, when ∆t is taken of order ∆x 2 , one can see that considering p = 2 and q = 4 leads to a total error in O(∆x 4 ). Otherwise, if ∆t = O(∆x), this is not true and fourth order accuracy is needed in time also to target error order 4. This is first illustrated by Figure 7, in which we show the error of the scheme for ncfl= 1 at t max = 0.05, with second and fourth order in space and different semi-implicit and implicit time schemes of order 1 (ARK1, SH, implicit), 2 (ARK2) and 4 (ARK4). The error is calculated by taking the L 2 norm of the difference between the analytical and the numerical solution. The graph on the left hand side, corresponding to a second order accurate in space, indicates results are second order accurate whatever the time scheme's order. This result is fully coherent taking into account that, for our diffusion problem, ∆t p = O(∆x 2p ) (see Equation (5.3) with ncfl= 1).

Here, this leads to ∆t p = O(∆x 2 ), independently of p ∈ N * . The graph on the right hand side, corresponding to a fourth order accuracy in space, shows results of second order accuracy when using time schemes of order 1, whereas time schemes of order 2 and 4 show order 4 accurate results. Again, the relation ∆t p = O(∆x 2p ) with p = 1 confirms that, with order 1 accuracy in time, we are only able to obtain order 2 accuracy results in space.

Then, another interesting point to observe here is the impact of ncfl on the error magnitude and its order. The stability condition for our diffusion problem in case of an explicit time scheme writes:

∆t < 1 4ρ(B)
× min(∆x 2 , ∆y 2 ).

Using ∆t values larger than the explicit time scheme's stability condition (i.e. ncfl > 1) is possible in case of semi-implicit and implicit time schemes, but leads to some drawbacks. The price to pay is increasing the error magnitude and potentially degrading the desired accuracy order. In Table 1: KCFL maximum to preserve order 1 in space. Here, L = 2, N x is the number of mesh points in both directions and n denotes the time scheme order. Assuming the time scheme is order p accurate and the spatial scheme is of order q, we have, by (5.3)

N x = 32 N x = 64 N x = 128 N x = 256 N x = 512 N x =
Err(N ) = O(∆t p + ∆x q ) = O ncfl × ∆x 2 4ρ(B) p + ∆x q .
Now, if we want the total error to be order q accurate, we have to satisfy the relation:

ncfl × ∆x 2 4ρ(B) p ≤ C∆x q ⇔ ncfl 4C 1/p ρ(B) ≤ ∆x (q-2p) 1/p ,
where C is a positive constant. Replacing ∆x by its value and setting KCFL = ncfl 4C 1/q ρ(B) leads to:

KCFL ≤ L N x (p-2q) 1/q .
Tables 1, 2 and 3 shows maximum KCFL values to preserve the order of approximation 1, 2, 3 and 4 in space respectively. These tables show the higher are the spatial resolution and the time scheme's order, the higher ncfl number we are able to choose. Let us see two examples:

• To preserve an order 4 accuracy with N x = 128, we look at the column corresponding to N x = 128 in Table 3. An order 4 time scheme allows to use KCFL = 64, whereas an order 2 one allows KCFL = 1 only. Additionally, with an order 1 time scheme, we see KCFL is restricted to approximately 2 × 10 -4 , leading to a very small time step. This observation motivates the use of high order time schemes.

• If we now look at the second line of Table 2 (time sheme of order 2), corresponding to an order 2 accuracy preservation, we see that multiplying the number of point by 2 allows to double the value of KCFL. Even more interesting, looking at the third line (time scheme of order 4), we report it is possible to (almost) triple KCFL value when doubling the number of points. This motivates again the use of high order in time, and also highlights semi-implicit scheme's benefits.

These last information about theoritical KCFL are cross checked by the test case we will consider in this section. Figure 8 shows the same type of plots than Figure 7, with fourth order accuracy in space only, and for different ncfl numbers. Again, the error the L 2 norm of the difference between the analytical and the numerical solution of the considered test case. With ncfl = 1, we see both ARK2 and ARK4 are fourth order accurate, whereas the other time schemes, of order 1, are of global order 2 only (explained before). When increasing ncfl parameter, we Inria Table 2: KCFL maximum to preserve order 2 in space. Here, L = 2, N x is the number of mesh points in both directions and n denotes the time scheme order. Table 3: KCFL maximum to preserve order 4 in space. Here, L = 2, N x is the number of mesh points in both directions and n denotes the time scheme order. Moreover, we observe that, in most of cases, increasing the number of points in space allows to improve the order accuracy of the measured error. We also see we are able to get good results in terms of error magnitude for ncfl up to 100 with order 4 in space accuracy and order 2 in time.

N x = 32 N x = 64 N x = 128 N x = 256 N x = 512 N x =
N x = 32 N x = 64 N x = 128 N x = 256 N x = 512 N x =
As a conclusion to these numerical tests, we clearly verified our methods on an analytical constant in space and time diffusion case. In addition, we motivated the use of high order methods in space and time when attempting to reach large time steps, while keeping control on the additional error that it produces.

Analytic test case: evaluation of numerical perpendicular diffusion error

Here, we consider a second analytical test case, designed to measure the numerical diffusion error for the different methods introduced in this work. This error is produced by the parallel diffusion (parallel to magnetic field lines), and may be amplified by the anisotropy. Hence, it might pollute the effective perpendicular diffusion in case of models such as Emedge3D one (see [START_REF] Beyer | Nonlinear dynamics of transport barrier relaxations in tokamak edge plasmas[END_REF][START_REF] Fuhr | Evidence from numerical simulations of transport-barrier relaxation in tokamak edge plasmas in the presence of electromagnetic fluctuations[END_REF]). The key idea is to initialize the temperature T and the magnetic field (b x , b y ) in such a way that the spatial operator cancels each other, leading to a theoretical constant temperature over time. However, the approximation methods used produce numerical diffusion error in most of the cases. Thus, a clear observable error magnitude between the analytical and numerical solution is the anomalous diffusion error. We consider the periodic 2D domain [-1, 1] × [-1, 1] and the equation to solve is The initial state is given by the relation:

∂ t T = ∇ • (B∇T ),
T (t = 0, x, y) = cos(mπx -nπy) = cos(µx -νy)
where µ = mπ and ν = nπ.

Let us know determine B such that ∇ • (B∇T ) = 0:

∇ • (B∇T ) = b 2 x ∂ 2 x cos(µx -νy) + b 2 y ∂ 2 y cos(µx -νy) + 2b x b y ∂ x ∂ y cos(µx -νy) = -b 2 x µ 2 cos(µx -νy) -b 2 y ν 2 cos(µx -νy) + 2b x b y µν cos(µx -νy) = -cos(µx -νy) (b x µ) 2 -2b x µb y ν + (b y ν) 2 = -T (b x µ -b y ν) 2
Hence, setting b x = ν and b y = µ gives ∇ • (B∇T ) = 0. In nuclear fusion codes, it is common to define the safety factor q(r) as the ratio m/n, typically varying from 1 at the plasma center to 3 or 4 at the edge. For this test, we consider m = 2 and n = 1 fixed, leading to a safety factor q = 2. As we want to measure only perpendicular diffusion error, we consider a unique value of ncfl= 1.

Figure 9 shows plots for numerical perpendicular thermal diffusion coefficient for RK2, ARK2, ARK4 and SH time schemes, different domain sizes and for order 2 and 4 spatial accuracy. As there should not be any diffusion, the better coefficient number is the lower. Except in the case of order 4 in space with SH time scheme, we see the thermal diffusion coefficients of the numerical diffusion error is independent of the chosen time scheme. We also clearly see the benefits of order 4 in space over order 2 on these two plots. Plus, SH time scheme is outperformed by other time schemes in the case of order 4 accuracy in space (right plot). This is because S.H. scheme is order 1 accuracy in time, which is not enough to benefit of order 4 in space.

Diffusion on a ring

In this part, we consider another case, with a non constant (in space) diffusion matrix. Let us define r = x 2 + y 2 . The magnetic field is defined by: with the initial temperature given by: T (t = 0, x, y) = 0.1 + 10e -((x-0.6) 2 +y 2 )/0.02 .

(b x , b y ) = -y r , x r
In this case, the temperature will reach a steady state over time, as magnetic field lines are closed describing concentric circles. This test case is also performed in [START_REF] Sharma | A fast semi-implicit method for anisotropic diffusion[END_REF]. The main change is the initial condition. Indeed, for our order 4 spatial scheme, second derivative need to be fully defined everywhere on the grid to avoid oscillation problems (arising with Dirac type initial conditions). As a second modification, our equivalent parallel diffusion coefficient χ , defining the speed of the diffusion in the parallel direction, is set to 1. For this test case, reference runs have been performed using a classical RK2 method, for two final time values: t max = 0.2 and 5 (see Figure 10).

Figures 11 shows results for this test case at t max = 0.2, for different semi-implicit and implicit time schemes and with varying ncfl numbers. The grid size is 256 × 256. Notice that, with χ = 1, reaching t max = 0.2 is equivalent to perform the test detailed in [START_REF] Sharma | A fast semi-implicit method for anisotropic diffusion[END_REF] with t max = 20 and χ = 0.01 (as there is no explicit perpendicular diffusion). Our results for the different time schemes are very similar to the reference run (on Figure 10). We observe typically the same half close ring shapes and temperature values than those observed in [START_REF] Sharma | A fast semi-implicit method for anisotropic diffusion[END_REF], enabling us to validate the different methods on this non constant in space diffusion problem for ncfl = 10. However, results for the two ARK penalization methods seems to be more sensitive to ncfl augmentation (see Figure 11 for ncfl = 1000): the diffusion speed appears to be slower. A possible reason to this last observation is the additional parameter λ needed for the penalization method, generating an extra error amplified by the ncfl factor. Notice that the ARK4 scheme is closer to the reference Inria solution than the ARK2 version, motivating the choice of a higher time order for the penalization method. We also remark on these figures that order 2 and order 4 in space methods are very close in terms of shapes and values. This is due to the low value of t max , not allowing the perpendicular error to develop.

Figure 12 shows results for this ring diffusion test case at a higher final time t max = 5 (equivalent to t max = 500 with χ = 0.01 in [START_REF] Sharma | A fast semi-implicit method for anisotropic diffusion[END_REF]), with ncfl = 20. It is now possible to see the impact of using an order 4 method in space, coupled together with a time scheme of order greater than 1. Indeed, we observe perpendicular diffusion error, produced by the parallel one (parallel here means parallel to magnetic field lines) is less important considering the order 4 method in space for ARK2 and ARK4 time schemes. The closed ring produced is thinner, indicating a better control of the error in the perpendicular direction. Also notice the maximum temperature value is higher for the order 4 in space methods (close to 0.8). The same result is shown on the reference run on Figure 10. These improvements are not noticeable on the implicit and the SH versions because these two schemes are only first order in time, not allowing to preserve the order 4 in space as showed in Section 5.1 (see also Figure 8).

Endly, Figures 13 and14 show results for a fixed ncfl = 1000 number, with a growing grid size. These graphs confirm observations made in Section 5.1 on the influence of ncfl on the error magnitude. Indeed we see that high CFL numbers are reachable only if the considered grid size is big enough to compensate the additional error in time, produced by a higher ∆t.

Constant diffusion on a periodic band

This test case aims to validate and to compare the different methods presented on a linear and constant case, i.e. the diffusion matrix is constant in space and time. It is to be considered as a preamble test to the nonlinear tests, as some characteristics are quite similar.

The domain we consider is [-1, 1] × [-1, 1] and the equation to solve is 2 ). The initial state is given by the relation:

∂ t T = ∇ • (B∇T ),
T (t = 0, x, y) = 1 + 3e -2r 2 if r < 2π 5 , 1
otherwise, where r = x 2 + y 2 . Periodic boundary conditions are used in both x and y directions, so that the temperature propagates on a periodic diagonal band in time.

This test is well designed to observe the numerical perpendicular diffusion error. Indeed, without explicit perpendicular diffusion, the temperature is supposed to reach a steady state, becoming constant along magnetic field lines. However, perpendicular diffusion occurs because of numerical errors committed by spatial scheme's approximations. For this test case, reference runs have been performed using again classical RK2 method, for two final time values: t max = 0.2 and 5 (see Figure 15).

For this test, conclusions are close to the ring test ones (in Section 5.3). Again, plots for short runs (see Figure 16) are very similar to reference runs, with penalization (ARK2 and ARK4) methods more sensitive to ncfl increase. Also, order 4 benefits is highlighted on Figure 17 depicting long runs are less diffusive in the perpendicular direction considering appropriate schemes, i.e. order 2 at least in time and order 4 in space.

Nonlinear analytic test case

The domain we consider is [-0.5, 0.5] × [-0.5, 0.5] and the equation to solve is

∂ t T = ∇ • (B(T )∇T ) + Q, with B(T ) = T 5/2 ǫ b b + (I -b b) and b b = b 2 x b x b y b x b y b 2 y , where b = (b x , b y ) = (1, 1
2 ), and Q ins a source term that will be described afterwards. As the initial condition, we consider the following function:

T (t, x, y) = c 1 + c 2 (sin(2πx) + ǫ cos(2πx) sin(2πy)) e -c3t ,
that we inject in the equation. Computing ∇ • (B(T )∇T ) leads to:

∇ • (B(T )∇T ) = 5 2ǫ T 3 2 (∂ x T ) 2 + 1 ǫ T 5 2 ∂ 2 x T + 5 2ǫ T 3 2 ∂ x T ∂ y T + 1 ǫ T 5 2 -1 ∂ x ∂ y T + 1 4 5 2ǫ T 3 2 (∂ y T ) 2 + 1 ǫ T 5 2 -3 ∂ 2 y T with ∂ x T = -2πc 2 ǫ sin(2πx) sin(2πy)e -c3t ∂ y T = 2πc 2 cos(2πy)(1 + ǫ cos(2πx))e -c3t ∂ 2 x T = -4π 2 c 2 ǫ cos(2πx) sin(2πy)e -c3t ∂ 2 y T = -4π 2 c 2 sin(2πy)(1 + ǫ cos(2πx))e -c3t ∂ x ∂ y T = -4π 2 c 2 ǫ sin(2πx) cos(2πy)e -c3t ,
and hence, defining Q(t, x, y) = -∇ • (B(T (t, x, y))∇T (t, x, y)) enables us to get an analytical solution of the equation. More than validating the different methods presented in this work on a full nonlinear diffusion equation, this test case also allows one to measure the numerical perpendicular diffusion produced by the parallel diffusion. Indeed, it is possible to emphasize this last one by setting ǫ as a small value (typically ǫ = 10 -3 , 10 -6 or 10 -9 , see [START_REF] Gunter | Finite element and higher order difference formulations for modelling heat transport in magnetised plasmas[END_REF][START_REF] Lozinski | Highly anisotropic temperature balance equation and its asymptotic-preserving resolution[END_REF]). Notice that papers dealing with the same numerical perpendicular diffusion problem often study the influence of the ratio χ ⊥ /χ . We consider here a normalized formulation with χ ⊥ set to 1 and ǫ = χ , allowing us to tackle the same problem. For the following tests, we only consider ARK2 and ARK4 time schemes. Indeed, the penalization methods allow to be more powerful in terms of execution time when considering order 4 in space, mostly because we are able to invert the implicit part of the scheme only one time and apply it then at each time step.

Figure 18 shows error orders in L 2 norm for order 2 and order 4 test cases. Both graphs highlight proper errors at orders 2 and 4, hence validating the methods for our nonlinear test case.

Figure 19 shows error as a function of ǫ value for ARK2 and ARK4 time schemes, order 2 and 4 in space and with different ncfl values. Here, the time step ∆t is evaluated independently of the value of ǫ. Hence, we are only able to consider unconditionally stable methods (i.e. semi-implicit or implicit methods). However, we observe again versions with spatial order 4 accuracy performs in terms of error magnitude. Also, it is clear that our schemes are not asymptotic preserving (AP) as in [START_REF] Lozinski | Highly anisotropic temperature balance equation and its asymptotic-preserving resolution[END_REF][START_REF] Mentrelli | Asymptotic-Preserving scheme for highly anisotropic non-linear diffusion equations[END_REF], but our aim is to design and compare efficient methods for a fixed and not high anisotropy (χ ⊥ /χ < 10 3 ) . 2 ). The initial condition is:

Nonlinear diffusion on a periodic band

T (t = 0, x, y) = 1 + 3e -2r 2 if r < 2π 5 , 1
otherwise.

where r = x 2 + y 2 . Periodic boundary conditions are used in both x and y directions, so that the temperature propagates on a periodic diagonal band in time. Figure 20 shows results for different ǫ values with fixed ARK2 time scheme, N x = 128, t max = 0.005 and ncfl = 1. On these plots, we observe ǫ values drive the speed of the parallel diffusion as desired. Indeed, for ǫ = 0.01, the high temperature band is more diffused as the one produced with ǫ = 0.5 for instance. We also observe the diffusion band becomes larger when using very low ǫ values. This is because, here again, the time step is chosen independently of ǫ, hence generating a higher error for low ǫ values and emphasizing the perpendicular numerical diffusion. Despite the error produced, we are able to reach a steady state regarding ǫ = 2 • 10 -3 , 10 -4 , 2 • 10 -4 , 10 -4 .

Conclusion

In this work, we propose, recall and compare some numerical schemes for diffusion problems. The semi-implicit scheme proposed in [START_REF] Sharma | A fast semi-implicit method for anisotropic diffusion[END_REF], an extension of the scheme proposed in [START_REF] Filbet | Numerical study of a nonlinear heat equation for plasma physics[END_REF] and the full implicit scheme are compared using a second and fourth order finite volume method for the spatial discretization. It appears that the SH scheme is efficient in its second order (in space) version but appears less efficient compared to the fourth order (in space) ARK with penalisation methods. Indeed, the numerical perpendicular diffusion error is decreased when higher order methods are used for the numerical tests we studied. Moreover, high order (in time) enables to recover precision and theoretical order even or large ncfl numbers. This is of great interest when long time simulations are needed, in order to deal with larger values of ∆t. Indeed, we were able to reach ncfl values up to some dozens, while keeping satisfactory error magnitudes values.

Several extensions to this work can be envisaged. First, a spatial depending λ can be tested in order to follow in the best way the anisotropy. Second, non periodic boundary conditions enabling to consider more physical test cases have to be included. This will be done using the strategy in [START_REF] Zhang | A fourth-order accurate finite-volume method with structured adaptive mesh refinement for solving the advection-diffusion equation[END_REF]. Then, the coupling with Asymptotic Preserving strategies by reformulating the initial problem would be interesting to ensure uniform precision, independently of the ratio χ ⊥ /χ , hence allowing to reach larger anisotropic diffusion problems. Finally, a fine study on the computational costs of the penalisation / ARK method could be performed, to measure the potential benefits on the execution times of this diffusion problem.

Appendix

Relations for spatial discretization

In this part, we recall some relations useful for the derivation of the numerical scheme. We consider one dimensional functions f := f (x) and h := h(x) and two-dimensional function g := g(x, y). A uniform one dimensional generic mesh is used: z k+1/2 = z k-1/2 + ∆z, for k ∈ Z. We denote by f k the cell averaged quantity

f k = 1 ∆z

Inria

The following proposition enables to express the face averaged derivative from cell averaged values, up to order 4. and the fourth order approximation

f ′ (z k+1/2 ) = 5 4∆z [f k+1 -f k ] - 1 12∆z [f k+2 -f k-1 ] + O(∆z 4 ).
Then, we deduce, for a function g = g(x, y), C y j ∂ x g(x i+1/2 , y)dy = 1 ∆x (g i+1,j -g i,j ) + O(∆x 2 ), for the second order = 5 4∆x (g i+1,j -g i,j ) -1 12∆x

(g i+2,j -g i-1,j ) + O(∆x 4 ),

for the fourth order and C x i ∂ y g(x, y j+1/2 )dx = 1 ∆y (g i,j+1 -g i,j ) + O(∆y 2 ), for the second order = 5 4∆y (g i,j+1 -g i,j ) -1 12∆y

(g i,j+2 -g i,j-1 ) + O(∆y 4 ),

for the fourth order.

Finally, we present some relations useful for the numerical approximations of the cross derivatives.

Proposition 7.4 For a function g = g(x, y), we get C x i ∂ x g(x, y j+1/2 )dx = g(x i+1/2 , y j+1/2 ) -g(x i-1/2 , y j+1/2 ) = (D y D x g) i,j , and C y j ∂ y g(x i+1/2 , y)dy = g(x i+1/2 , y j+1/2 ) -g(x i+1/2 , y j-1/2 ) = (D x D y g) i,j . We used the notations (D x g) i,j = 1 ∆x (g i+1,j -g i,j ) + O(∆x 2 ), for the second order = 7 12∆x (g i+1,j + g i,j ) -1 12∆x (g i+2,j + g i-1,j ) + O(∆x 4 ),

for the fourth order.

Inria

j

  ∂ x T (x i+1/2 , y)dy and C y j b xx (x i+1/2 , y)dy. They are respectively approximated using Proposition 7.3 and Proposition 7.2. The same is true for the third line which leads to integral of the form C x i ∂ y T (x, y j+1/2 )dx (which needs Proposition 7.3) and C x i b yy (x, y j+1/2 )dx (which needs Proposition 7.2). For the terms of the second line, integrals of type C y j ∂ y T (x i+1/2 , y)dy require Proposition 7.4 whereas terms of type C y j b xy (x i+1/2 , y)dy require Proposition 7.2. Finally, the fourth line of (3.2) makes appear terms like C x i ∂ x T (x, y j+1/2 )dx (which are approximated following Proposition 7.4 38) and C x i b xy (x, y j+1/2

b

  xx (x i+1/2 , y)dy where b xx depends on the cell averaged unknown T i,j . Obviously, when b xx

  Then, using Proposition 7.2, one can compute face averaged values C y j b xx (T )(x i+1/2 , y)dy, as a function of cell averaged values b xx (T ) i,j . Let us remark that second order approximations do not involve these computations since we have b xx

  which proves the unconditional stability of the scheme. With the choice b xx = b yy = b xy = 1/2, we plot in Figure 1 the coefficient |r| = |r 1 r 2 | for different values of ncfl.

Figure 1 :InriaFigure 2 :

 12 Figure 1: Amplification factor |r(t)| for ncfl= 0.01, 0.1, 10, 100, 1000 for SH scheme of order 2.

Figure 3 :

 3 Figure 3: Amplification factor |r(t)| for ncfl= 1, 10, 100, 1000 for SH scheme of order 4.

Figure 4 :

 4 Figure 4: Amplification factor |r(t)| for ncfl= 1, 10, 100, 1000 for ARK scheme of order 2 in space.

Figure 5 :

 5 Figure 5: Amplification factor |r(t)| for ncfl= 1, 10, 100, 1000 for ARK scheme of order 4 in space.

Figure 6 :InriaFigure 7 :

 67 Figure 6: Amplification factor |r(t)| for ncfl= 1, 10, 100, 1000 "SH-ARK" scheme of order 2 in space.

with B = b b constant in space and time and b b = b 2 x

 2 b x b y b x b y b 2 y .

Figure 8 :

 8 Figure 8: Error orders for semi-implicit and implicit time schemes and for ncfl = 1, 10, 100, 500 and 1000 with order 4 accuracy in space.

Figure 9 :

 9 Figure 9: Numerical perpendicular thermal diffusion coefficient : spatial order 2 (left) and 4 (right).

Figure 10 :

 10 Figure 10: Reference diffusion on a ring at t max = 0.2 and 5, with RK2 time scheme.

with B = b b and b b = b 2 x b x b y b x b y b 2 y.

 22 where (b x , b y ) = (1, 1

Figure 11 :

 11 Figure 11: Diffusion on a ring at t max = 0.2, for semi implicit and implicit time schemes, with N x = 256 and ncfl= 10 and 1000.

Inria

Figure 12 :

 12 Figure 12: Diffusion on a ring at t max = 5, for different time scheme, with N x = 128 and ncfl = 20. These plots allow one to observe the benefit of order 4 spatial scheme. It reduces numerical perpendicular diffusion error.

Figure 13 :

 13 Figure 13: Diffusion on a ring at t max = 0.2, for implicit and SH time schemes, with N x = 64, 128, 256 and ncfl= 1000 fixed.

Inria

Figure 14 :

 14 Figure 14: Diffusion on a ring at t max = 0.2, for ARK2 and ARK4 time schemes, with N x = 64, 128, 256 and ncfl= 1000 fixed.

Figure 15 :

 15 Figure 15: Reference diffusion on a periodic band at t max = 0.2 and 5, with RK2 time scheme.

Figure 16 :

 16 Figure 16: Diffusion on a ring at t max = 0.2, for semi implicit and implicit time schemes, with N x = 256 and ncfl= 100 and 1000.

Figure 17 :

 17 Figure 17: Diffusion on a periodic band at t max = 5, for different time scheme, with N x = 128 and ncfl = 20.

Inria

Figure 18 :

 18 Figure 18: Error orders for ARK2 and ARK4 time schemes, ncfl = 1 with order 2 and 4 accuracy in space and ǫ = 1.

Figure 19 :

 19 Figure 19: Error for ARK2 and ARK4 time schemes, ncfl = 1 with order 2 and 4 accuracy in space and different ǫ values.

= b 2 x b x b y b x b y b 2 y,

 22 This test case is very close to the one presented in Section 5.4, except we consider now the same nonlinear matrix diffusion as defined in Section 5.5. The domain we consider is [-1, 1] × [-1, 1] and the equation to solve is∂ t T = ∇ • (B(T )∇T ) + Q, with B(T ) = T5/2 ǫ b b + (I -b b) and b b Inria where (b x , b y ) = (1, 1

z k+1/ 2 z k-1/ 2 f

 22 (z)dz.

Figure 20 :

 20 Figure 20: Nonlinear diffusion on a periodic band at t max = 0.005, for different ǫ values, with N x = 128 and ncfl = 1.

Proposition 7 . 3

 73 For any function f , we have the following second order approximationf ′ (z k+1/2 ) = 1 ∆z (f k+1 -f k ) + O(∆z 2 ).

  we assume ∆x = ∆y = L/N x with L the domain length and N x the number of points, equal in each direction. Let us now consider the error Err(N ) := ||T N -T (t max )|| at time t max = N ∆t (where N is the number of iterations).

							1024
	n= 1	16.0	32.0	64.0	128.0	256.0	512.0
	n= 2	64.0	181.0	512.0	1448.2	4096.0	11585.2
	n= 4	128.0	430.5	1448.2	4871.0	16384.0	55109.0
	the following,

RR n°8560

Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

We recall the notations introduced above C x i = [x i-1/2 , x i+1/2 ] and C y j = [y j-1/2 , y j+1/2 ], with ∆x and ∆y the mesh size in the directions x and y so that the cell averaged writes g i,j = 1 ∆x∆y C x i C y j g(x, y)dxdy.

The first one enables to express the product of the integral of two functions as a product of the integral of the function.

whereas a fourth order approximation of I gives

The following proposition enables to express the face averaged values from the cell averaged values, up to order 4. Proposition 7.2 A second order approximation of the face average of f from cell averaged

whereas a fourth order approximation gives

Then, we deduce, for a function g = g(x, y),

), for the second order

for the fourth order and

), for the second order

for the fourth order.

RR n°8560

and (D y g) i,j = 1 ∆y (g i,j+1 -g i,j ) + O(∆y 2 ), for the second order

for the fourth order, so that

for the second order (7.4)

+ O(∆x 4 ) + O(∆y 4 ), for the fourth order. (7.5)

Linear stability analysis: fourth order

We detail in this part the computations of the linear stability analysis of the two time integrators we focus on: SH and ARK schemes.

Semi-implicit scheme (SH)

We write the numerical scheme of order four in space and couple with the semi-implicit scheme.

We then define the first step of the scheme as (using the notation above)

whereas the second step writes

where

and (D X D y T ) i,j given by (7.5).

RR n°8560

Following the Von Neumann analysis, we inject a plane wave solution T (t, x, y) = r(t) exp -i(kxx+kyy) to get the amplification factor r 1 for the first step (with ncfl= 4∆t/∆x 2 )

,

whereas the amplification factor r 2 for the second step is .

The amplification factor r = r 1 r 2 is plotted in Figure 3 for different values of ncfl.

Additive Runge Kutta

We write the numerical scheme of order four in space and couple with the penalisation technique. We then define the numerical scheme as

where (D x xT ) i,j , (D x xT ) i,j are given by (7.6) and (7.7), whereas (D x D y T n ) i,j is given by (7.5).

Following the Von Neumann analysis, we inject a plane wave solution T (t, x, y) = r(t) exp -i(kxx+kyy) to get the amplification factor r = r 1 /r 2 (with ncfl= 4∆t/∆x 2 ) with The amplification factor r = r 1 /r 2 is plotted in Figure 5 for different values of ncfl.