Fabrice Deluzet

Jacek Narski

A TWO FIELD ITERATED ASYMPTOTIC-PRESERVING METHOD FOR HIGHLY ANISOTROPIC ELLIPTIC EQUATIONS

Keywords: Anisotropic diffusion, asymptotic preserving scheme, iterative method AMS subject classifications. 65N30

A new two field iterated Asymptotic-Preserving method is introduced for the numerical resolution of strongly anisotropic elliptic equations. This method does not rely on any integration of the field defining the anisotropy. It rather harnesses an auxiliary variable removing any stiffness from the equation. Compared to precedent realizations using the same approach, the iterated method allows for the resolution of each field independently within an iterative process to converge the two unknowns. This brings advantages in the computational efficiency of the method for large meshes, a better scaling of the matrices condition number with respect to the mesh refinement as well as the ability to address complex anisotropy topology including closed field lines.

1. Introduction. The present paper is aimed at introducing a new Asymptotic-Preserving scheme for the resolution of singular perturbation problems stemming from strongly anisotropic elliptic equations. This type of equations are representative of plasma physics evolution under large magnetic fields such as Tokamak plasmas [START_REF] Degond | Asymptotic-preserving schemes for fluid models of plasmas[END_REF][START_REF] Degond | Asymptotic-Preserving methods and multiscale models for plasma physics[END_REF]. Here the focus is made on a simplified model problem containing the main difficulty characterizing these equations but without all the complexity of the physical background. This simplified context allows the construction of analytic solutions which are used to assess the effectiveness of the numerical method introduced herein. Let b denote the vector field providing the direction of the magnetic field, b verifying |b| = 1, the model problem writes

         -∇ • (A ε ∇u ε ) = f ε in Ω, n • A ε ∇u ε = 0 on Γ N , u ε = 0 on Γ D , (1)
where n is the outward normal to the domain, Γ N ∪ Γ D the domain boundary, with b • n = 0 on Γ D and b • n = 0 on Γ N . The anisotropy of the problem is defined by the diffusion matrix A ε related to the vector field b by two positive functions A and A ⊥ with

A ε = 1 ε A b ⊗ b + (Id -b ⊗ b)A ⊥ (Id -b ⊗ b) . (2)
In this equation Id is the identity matrix, the tensor product being denoted ⊗. The parameter ε -1 defines the strength of the anisotropy.

The difficulty addressed in this paper is related to the singular nature of the problem. Indeed in the limit of infinite anisotropy strength (ε → 0) the system (1) is degenerate. Indeed, the differential operator in the elliptic equation reduces to the dominant operator (the derivatives carried by ε -1 in (1)) which is supplied with Neumann like boundary conditions. This degenerate system admits an infinite amount of solutions, any function with no gradient along b being in the kernel of the dominant operator. The derivation of efficient numerical methods for the approximation of this class of problems is a difficult task. The straight discretization of [START_REF] Amestoy | A Fully Asynchronous Multifrontal Solver Using Distributed Dynamic Scheduling[END_REF] gives rise to system matrices with condition number blowing up with the increase of the anisotropy strength. This is outlined in precedent works (see [START_REF] Degond | Duality-based asymptotic-preserving method for highly anisotropic diffusion equations[END_REF] for numerical investigations or [START_REF] Lozinski | Numerical analysis of an asymptotic-preserving scheme for anisotropic elliptic equations[END_REF] and [START_REF] Yang | Iterative solvers for elliptic problems with arbitrary anisotropy strengths[END_REF] for an analysis). Therefore these approaches are limited to reduced anisotropy strength.

A way to circumvent this difficulty is to develop Asymptotic-Preserving methods as introduced in [START_REF] Jin | Asymptotic preserving (AP) schemes for multiscale kinetic and hyperbolic equations: a review, Lecture Notes for Summer School on "Methods and Models of Kinetic Theory[END_REF] for a different context. Actually a well posed system can be derived to compute uniquely the solution in the limit of infinite parallel diffusion. The aim of such method is to guarantee that the discrete system is consistent with this well posed problem for vanishing ε rather than with the degenerate one. This ensures that the condition number of the system matrix remains bounded independently of the anisotropy strength.

In precedent works, different AP schemes have been derived for this class of problems. The first iterations were devoted to anisotropy directions aligned with one coordinate [START_REF] Degond | Duality-based asymptotic-preserving method for highly anisotropic diffusion equations[END_REF][START_REF] Besse | Efficient numerical methods for strongly anisotropic elliptic equations[END_REF]. This requirement has been released in [START_REF] Degond | An asymptotic-preserving method for highly anisotropic elliptic equations based on a micro-macro decomposition[END_REF] and extended to closed field lines in [START_REF] Narski | Asymptotic Preserving scheme for strongly anisotropic parabolic equations for arbitrary anisotropy direction[END_REF]. In all these works the problem is reformulated into a two field problem based on a decomposition of the solution into a microscopic and a macroscopic component. This reformulated two field system offers the advantage of embedding the limit problem. Hence, the limit ε → 0 is a regular limit in this set of equations. However the decomposition of the solution is not unique and different numerical methods can be derived according to the choices implemented in this reformulated system. The present work aims at exploring further the possibilities offered by a different decomposition. Note that the method developed herein does not rely on any geometrical procedure, requiring an integration along the b-field lines as proposed by other authors [START_REF] Chacón | An asymptotic-preserving semilagrangian algorithm for the time-dependent anisotropic heat transport equation[END_REF][START_REF] Tang | An asymptotic preserving method for strongly anisotropic diffusion equations based on field line integration[END_REF][START_REF] Wang | Uniformly Convergent Scheme for Strongly Anisotropic Diffusion Equations with Closed Field Lines[END_REF].

The main goal of the present work is to correct some of the weaknesses of the precedent realizations, namely the Micro-Macro and the stabilized Micro-Macro methods. The first weakness is related to the structure of the discrete system issued from the discretization of the reformulated problem. So far, this system strongly couples the equations providing both components and is therefore solved at once. We propose a different method referred to as "two field iterated Asymptotic-Preserving" method which offers the ability to solve each component independently, in an iterative process. The system solved for each component is the same mildly anisotropic problem parameterized by a numerical parameter ε 0 ε with a different source term for every component. This gains an improved efficiency in terms of computational resources compared to the direct resolution of the two field system implemented for the Micro-Macro methods. The second advantage of this new method is related to the conditioning of the system matrix. The linear systems issued from precedent AP (Micro-Macro) methods [START_REF] Degond | An asymptotic-preserving method for highly anisotropic elliptic equations based on a micro-macro decomposition[END_REF][START_REF] Lozinski | Numerical analysis of an asymptotic-preserving scheme for anisotropic elliptic equations[END_REF] have a condition number scaling as 1/h 4 , h denoting the typical mesh size. The two field iterated method introduced herein requires only the resolution of linear systems with a condition number scaling as 1/(ε 0 h 2 ). An additional advantage is the ability to carry out numerical approximations with closed field lines. This is a difficulty that can not be addressed by the Micro-Macro AP scheme [START_REF] Degond | An asymptotic-preserving method for highly anisotropic elliptic equations based on a micro-macro decomposition[END_REF]. Indeed this numerical method requires that all the field lines cross the domain boundary. It should be pointed out that a solution has been proposed in [START_REF] Narski | Asymptotic Preserving scheme for strongly anisotropic parabolic equations for arbitrary anisotropy direction[END_REF] in the frame of the "stabilized" Micro-Macro scheme. It consists in introducing a stabilization operator small enough not to deteriorate the precision of the numerical method. The main difficulty with this approach lies in the choice of the stabilization operator scaling. Indeed, it should be kept large enough to preserve a good condition-ing of the system matrix but small enough to be comparable to the truncation error of the discretizations. This prevents from using the stabilized method with high order methods. Another approach is proposed in [START_REF] Wang | Uniformly Convergent Scheme for Strongly Anisotropic Diffusion Equations with Closed Field Lines[END_REF] relying on an integration along the magnetic field lines. This method is derived without coordinates or meshes related to the magnetic field but requires precise integration along the field lines only implemented with analytic expressions of the magnetic field. The two field iterated method is free from these weaknesses. We also show that, the parameter ε 0 can be chosen in a wide range of values preventing the so-called locking effect [START_REF] Babuška | On Locking and Robustness in the Finite Element Method[END_REF] and securing a fast convergence of the iterations as well as a good conditioning of the linear systems. The numerical method is also free from the perpendicular dynamic pollution by the parallel one, reported by other authors in very similar frameworks [START_REF] Günter | Finite element and higher order difference formulations for modelling heat transport in magnetised plasmas[END_REF][START_REF] Günter | Modelling of heat transport in magnetised plasmas using non-aligned coordinates[END_REF][START_REF] Van Es | Finite-difference schemes for anisotropic diffusion[END_REF][START_REF] Van Es | Finite-volume scheme for anisotropic diffusion[END_REF].

The outlines of the paper are the following. The problem at hand in the present work is stated in Section 2 with highlights on the singular nature of the limit ε → 0. The two field iterated AP method is then introduced and the convergence of the iterative procedure is demonstrated. Finally, emphasizes are made on how this new method compares to precedent works. Numerical investigations are carried out in Section 3. Different setups are proposed to asses the effectiveness of the method. The locking effect is first investigated and the robustness of the method with respect to this classical issue is outlined. The efficiency of the two field iterated method is benchmarked against the Micro-Macro scheme. This demonstrates tremendous gains for large meshes. Two other test cases are finally proposed. The second one is a diffusion in a ring similarly to computations performed in [START_REF] Crouseilles | Comparison of Numerical Solvers for Anisotropic Diffusion Equations Arising in Plasma Physics[END_REF][START_REF] Sharma | A fast semi-implicit method for anisotropic diffusion[END_REF] but proposed here with anisotropy strength much more severe. The last test case is aimed at demonstrating the ability to carry out accurate numerical approximations in frameworks including closed field lines.

2.

The anisotropic problem and its asymptotic-preserving formulation.

Γ D = {x ∈ Γ | b(x) • n = 0}, Γ N = Γ \ Γ D , (3)
where n is the outward normal to Ω. Let us also decompose any vector v ∈ R d , gradients ∇φ, with φ(x) a scalar function, and divergence ∇ • v into a part parallel to the anisotropy direction and a part perpendicular to it with:

v := (v • b)b , v ⊥ := (Id -b ⊗ b)v , such that v = v + v ⊥ , ∇ φ := (v•∇φ)∇u , ∇ ⊥ φ := (Id -b⊗b)∇φ , such that ∇φ = ∇ φ + ∇ ⊥ φ , ∇ • v := ∇ • v , ∇ ⊥ • v := ∇ • v ⊥ , such that ∇•v = ∇ •v + ∇ ⊥ •v ,
where we denoted Id the identity matrix and ⊗ the vector tensor product. The following notations and definitions will be helpful in the sequel.

Definition 1. Let V and G be the functional spaces defined by

V = {v ∈ H 1 (Ω) : v| Γ D = 0} , (4) G = {v ∈ V : ∇ v = 0} . (5) For any function φ ∈ V, ε 0 ∈ R, ε 0 > 0, A ∈ C ∞ ( Ω) a positive function and A ⊥ ∈ M d×d (C ∞ ( Ω)) a matrix satisfying A 0 ||v|| 2 ≤ v T A ⊥ v ≤ A 1 ||v|| 2 , ∀v ∈ R d (6)
for some positive constants A 0 and A 1 , we introduce the operators ∆ , ∆ ⊥ and ∆ ε0 defined as

∆ φ = ∇ • A ∇ φ , (7a) ∆ ⊥ φ = ∇ ⊥ • (A ⊥ ∇ ⊥ φ) , (7b) ∆ ε0 φ = ∆ φ + ε 0 ∆ ⊥ φ ; (7c) and for (u, v) ∈ V × V the associated bilinear forms a (u, v) = Ω A ∇ u • ∇ vdx , (8a) a ⊥ (u, v) = Ω (A ⊥ ∇ ⊥ u) • ∇ ⊥ vdx . ( 8b
)
Finally the matrix A ε0 is introduced with

(9a) A ε0 = A (b ⊗ b) + ε 0 (Id -b ⊗ b) A ⊥ (Id -b ⊗ b)
and the induced norm

||u|| 2 ε0 = a (u, u) + ε 0 a ⊥ (u, u) . (9b)
2.2. The singular perturbation problem. The problem studied in this paper is the following: find u ε such that

         - 1 ε ∆ u ε -∆ ⊥ u ε = f in Ω, 1 ε n • A ∇ u ε + n ⊥ • (A ⊥ ∇ ⊥ u ε ) = 0 on Γ N , u ε = 0 on Γ D , (10)
This problem is refereed to as a singular perturbation problem, because of its degeneracy for vanishing ε. Indeed, setting ε to 0 in [START_REF] Degond | Duality-based asymptotic-preserving method for highly anisotropic diffusion equations[END_REF], the problem reduces to

         -∆ u 0 = 0 in Ω, n • A ∇ u 0 = 0 on Γ N , u 0 = 0 on Γ D , (11)
which admits an infinite number of solutions as any function v that is constant in the direction of anisotropy (v ∈ G) solves this problem. The limit of the solution can be computed by multiplying [START_REF] Degond | Duality-based asymptotic-preserving method for highly anisotropic diffusion equations[END_REF] by a test function v ∈ G, integrating by parts over Ω and then let ε → 0. This leads to the following, well posed problem: find u 0 ∈ G such that

Ω (A ⊥ ∇ ⊥ )u 0 • ∇ ⊥ v = Ω f v , ∀v ∈ G , (12)
which defines a weak formulation of the limit problem. The difficulty when dealing with the numerical approximation of (1) consists in imposing the consistency of the scheme with the limit problem [START_REF] Günter | Finite element and higher order difference formulations for modelling heat transport in magnetised plasmas[END_REF] rather than the degenerate one [START_REF] Degond | An asymptotic-preserving method for highly anisotropic elliptic equations based on a micro-macro decomposition[END_REF] when ε → 0. Standard discretizations of the problem [START_REF] Degond | Duality-based asymptotic-preserving method for highly anisotropic diffusion equations[END_REF] are not compliant with this property. The condition number of the associated system matrices are increasing with the anisotropy strength ε -1 . This translates that the numerical methods provide a discretization of the degenerate problem for vanishing ε-values. To address this issue, the philosophy of Asymptotic Preserving schemes relies on a discretization of a suitable reformulated problem. This system is equivalent to the problem [START_REF] Degond | Duality-based asymptotic-preserving method for highly anisotropic diffusion equations[END_REF] for ε > 0, however the limit problem [START_REF] Günter | Finite element and higher order difference formulations for modelling heat transport in magnetised plasmas[END_REF] is recovered from the reformulated system when ε is set to 0.

Another difficulty encountered when dealing with the numerical resolution of anisotropic problems is the so-called locking phenomenon [START_REF] Babuška | On Locking and Robustness in the Finite Element Method[END_REF]. To highlight this issue in the present framework, let us again consider the reduced problem [START_REF] Degond | An asymptotic-preserving method for highly anisotropic elliptic equations based on a micro-macro decomposition[END_REF]. This problem states that the solution has no gradient along b for vanishing ε. If the discrete space does not contain functions that are constant in the direction of the anisotropy, then the numerical approximation of this problem does not converge to the solution of the problem. It is important to point out that the locking is not related to the fact that the reduced problem is ill posed on the continuous level but to the coarse approximation properties of the discrete functional space. That is the case, for example, when either unstructured (triangular) meshes or rectangular Cartesian grids with variable anisotropy directions are used with low order numerical methods. For small non zero values of ε, large enough to preserve a good conditioning of the matrix related to the discretized version of the [START_REF] Degond | Duality-based asymptotic-preserving method for highly anisotropic diffusion equations[END_REF], the locking phenomenon is manifested in the discrete solution converging to zero as ε gets smaller. This feature will be illustrated in the section devoted to the numerical investigations.

A two field iterated Asymptotic-Preserving method.

Let us now propose a two step iterative method to solve the singular perturbation problem [START_REF] Degond | Duality-based asymptotic-preserving method for highly anisotropic diffusion equations[END_REF]. Let us consider ε0 smaller than one but big enough so that the singular perturbation problem for ε = ε0 is not yet singular nor the discretized system suffers from locking. Let us define ε 0 = max{ε 0 , ε} so that ε 0 is never smaller than ε.

Let us first observe that the source of the numerical issues in the resolution of the original problem [START_REF] Degond | Duality-based asymptotic-preserving method for highly anisotropic diffusion equations[END_REF] is the dominant derivative, multiplied by ε -1 , in the direction of the anisotropy. The idea behind the herein proposed scheme relies on the introduction of an additional variable that fulfils the following relation: ε∆ q = ∆ u. This operation allows to eliminate the stiff term from the equation, preventing by this means the degeneracy of the equation. The two field system becomes:

-∆ q -∆ ⊥ u = f , -∆ u = -ε∆ q , (13)
supplied with the boundary conditions precised in [START_REF] Degond | Duality-based asymptotic-preserving method for highly anisotropic diffusion equations[END_REF] for both u and q. This system does not have a unique solution as q is defined up to a function constant in the direction of the anisotropy. Let us now multiply the first equation by ε 0 and add it to the second one to get: [START_REF] Jin | Asymptotic preserving (AP) schemes for multiscale kinetic and hyperbolic equations: a review, Lecture Notes for Summer School on "Methods and Models of Kinetic Theory[END_REF] -

∆ ε0 u = ε 0 f + (ε 0 -ε)∆ q,
allowing to compute u uniquely if q is known. The next step consists in decoupling the problems and solve the two resulting equations in an iterative manner, finding first an approximation to u using q computed in the previous step, then recompute q and repeat until convergence. This yields the following iterations

-∆ ε0 u n+1 = ε 0 f + (ε 0 -ε)∆ q n , -∆ q n+1 = f + ∆ ⊥ u n+1 . ( 15
)
The second equation of this iterative scheme is not yet invertible. Let us now add the term -ε 0 ∆ ⊥ q n+1 to the left hand side and subtract its equivalent for q n from the left hand side. The resulting problem for q n+1 has a unique solution for given q n and u n+1 . Finally, the two field iterated method is defined in the following way:

         -∆ ε0 u n+1 = ε 0 f + (ε 0 -ε)∆ q n in Ω, n • A ε0 ∇u n+1 = -(ε 0 -ε)n • A ∇ q n on Γ N , u n+1 = 0 on Γ D , (16)
         -∆ ε0 q n+1 = f + ∆ ⊥ (u n+1 -ε 0 q n ) in Ω, n • A ε0 ∇q n+1 = -n • A ⊥ ∇ ⊥ u n+1 -ε 0 q n on Γ N , q n+1 = 0 on Γ D , (17)
where q n+1 is an auxiliary variable and u n+1 the approximation to u ε . In this method, the original strongly anisotropic elliptic problem ( 10) is replaced by a set of two only mildly anisotropic equations parameterized by ε 0 ε. Moreover, the matrix to be inverted in the first step [START_REF] Narski | Asymptotic Preserving scheme for strongly anisotropic parabolic equations for arbitrary anisotropy direction[END_REF] of the iterative method is the same as in the final step [START_REF] Ratnani | Anisotropic Diffusion in Toroidal geometries[END_REF], the only difference is in the right hand side of the equation. That is to say, the matrix has to be factorized only once, the rest of the iterative scheme is a fast triangular system solve. This method does note require any discretization of the space G (functions constant in the direction of the anisotropy), which can be complicated for generic field b. To be complete, the variational formulation of the iterative scheme [START_REF] Narski | Asymptotic Preserving scheme for strongly anisotropic parabolic equations for arbitrary anisotropy direction[END_REF][START_REF] Ratnani | Anisotropic Diffusion in Toroidal geometries[END_REF] is stated:

Find (q n+1 , u n+1 ) ∈ V × V such that a (u n+1 , v) + ε 0 a ⊥ (u n+1 , v) = ε 0 (f, v) -(ε 0 -ε)a (q n , v), ∀v ∈ V, (18)
a (q n+1 , w) + ε 0 a ⊥ (q n+1 , w) = (f, w) -a ⊥ (u n+1 -ε 0 q n , w), ∀w ∈ V. ( 19
)
Let us now prove that the iterative scheme [START_REF] Narski | Asymptotic Preserving scheme for strongly anisotropic parabolic equations for arbitrary anisotropy direction[END_REF][START_REF] Ratnani | Anisotropic Diffusion in Toroidal geometries[END_REF] converges and that the limit solution solves the original singular perturbation problem.

Theorem 2. For any (q 0 , u 0 ) ∈ V × V, the sequence (q n , u n ) n>0 defined by the iterative method [START_REF] Narski | Asymptotic Preserving scheme for strongly anisotropic parabolic equations for arbitrary anisotropy direction[END_REF][START_REF] Ratnani | Anisotropic Diffusion in Toroidal geometries[END_REF] converges to a solution (q, ū). The component ū of the stationary point solves uniquely the initial singular perturbation problem [START_REF] Degond | Duality-based asymptotic-preserving method for highly anisotropic diffusion equations[END_REF] for ε > 0 and the limit problem [START_REF] Günter | Finite element and higher order difference formulations for modelling heat transport in magnetised plasmas[END_REF] when ε = 0.

To prove Theorem 2, the following lemmas and proposition are necessary.

Lemma 3. The operator ∆ ε0 is invertible on V. The eigenvalues of the operator ∆ -1 ε0 ∆ are real non negative and bounded by 1. The eigen functions ν 0 associated to the null eigenvalue belong to the kernel of the operator ∆ : ν 0 ∈ G. Lemma 4. The iterative method defined by Eqs. [START_REF] Narski | Asymptotic Preserving scheme for strongly anisotropic parabolic equations for arbitrary anisotropy direction[END_REF][START_REF] Ratnani | Anisotropic Diffusion in Toroidal geometries[END_REF] yields the following recurrence

q n+1 = A I q n -∆ -1 ε0 ∆ ∆ -1 ε0 f (20)
for n ≥ 1, the iteration operator A I being defined as

A I = 1 - ε ε 0 ∆ -1 ε0 ∆ - ε 0 -ε ε 0 ∆ -1 ε0 ∆ 2 . ( 21
)
The eigenvalues of A I , denoted i , are real with i ∈ [0, 1]. The eigenfunctions associated to the largest eigenvalue i = 1 is in the kernel of the operator ∆ .

Lemma 5 (Orthogonality of q n+1 -q n with respect to w ∈ G). For any q 0 ∈ V all functions in the sequence (q n ) n≥0 issued from the iterative method [START_REF] Narski | Asymptotic Preserving scheme for strongly anisotropic parabolic equations for arbitrary anisotropy direction[END_REF][START_REF] Ratnani | Anisotropic Diffusion in Toroidal geometries[END_REF] differ from each other only by a function orthogonal to G, the space of functions constant in the direction of anisotropy with respect to the H 1 seminorm. That is to say, for any i, j ≥ 0 the difference q j -q i is orthogonal to G with respect to the H 1 seminorm. Moreover, if ∇ ⊥ q 0 = 0 than q n is orthogonal to G with respect to the H 1 seminorm for all n ≥ 0. Proposition 6. For any fixed point (ū, q) of the iterative method defined by Eqs. [START_REF] Narski | Asymptotic Preserving scheme for strongly anisotropic parabolic equations for arbitrary anisotropy direction[END_REF][START_REF] Ratnani | Anisotropic Diffusion in Toroidal geometries[END_REF], the component ū is the solution of the singular perturbation problem [START_REF] Degond | Duality-based asymptotic-preserving method for highly anisotropic diffusion equations[END_REF] for ε > 0 and of the limit problem [START_REF] Günter | Finite element and higher order difference formulations for modelling heat transport in magnetised plasmas[END_REF] for ε = 0.

Proof of lemma 3. The operator ∆ ε0 is invertible due to standard elliptic arguments. The eigensystem of the operator ∆ -1 ε0 ∆ is defined by the problem:

Find λ i ∈ R and q i ∈ V such that ∆ -1 ε0 ∆ q i = λ i q i , (22)
or equivalently

∆ q i = λ i ∆ ε0 q i . ( 23
)
Multiplication by q i (or by q * i , if q i is assumed to have complex values and λ i ∈ C) and integration by parts yield

λ i = a (q i , q i ) a (q i , q i ) + ε 0 a ⊥ (q i , q i ) . (24)
Clearly all eigenvalues are real and between 0 and 1.

Proof of lemma 4. Thanks to Eq. ( 16) it follows that, on the one hand

u n+1 = -ε 0 ∆ -1 ε0 f -(ε 0 -ε)∆ -1 ε0 ∆ q n (25)
and, on the other hand

f + ∆ ⊥ u n+1 = 1 ε 0 -∆ u n+1 -(ε 0 -ε)∆ q n .
Plugging this identity into Eq. ( 17) yields, (26)

∆ ε0 q n+1 = ∆ ε0 q n + 1 ε 0 ∆ u n+1 - ε ε 0 ∆ q n .
After some algebra and using Eq. ( 25), the recurrence relation ( 20) between q n+1 and q n is recovered. The eigenvalues of the iteration operator A I defined by Eq. ( 21)

verify i = 1 - ε ε 0 λ i - ε 0 -ε ε 0 λ 2 i ( 27
)
where λ i are the eigenvalues of the operator ∆ -1 ε0 ∆ characterized in lemma 3. Note that i is a decreasing function of λ i for ε 0 > ε, with i = 1 for λ i = 0 and i = 0 for λ i = 1.

Proof of lemma 5. Let us first prove that q n+1 -q n is orthogonal to the space G with respect to the H 1 seminorm. Let us plug w ∈ G into [START_REF] Tang | An asymptotic preserving method for strongly anisotropic diffusion equations based on field line integration[END_REF] to get

ε 0 a ⊥ (q n+1 , w) = (f, w) -a ⊥ (u n+1 , w) + ε 0 a ⊥ (q n , w) , ∀w ∈ G. (28)
Owing to the identity a ⊥ (u n+1 , w) = (f, w) and thanks to [START_REF] Sharma | A fast semi-implicit method for anisotropic diffusion[END_REF] evaluated with v = w ∈ G, the following relation is derived

a ⊥ (q n+1 -q n , w) = 0 , ∀w ∈ G, (29)
which finally leads to

a (q n+1 -q n , w) + a ⊥ (q n+1 -q n , w) = 0 , ∀w ∈ G. ( 30
)
This proves that q n+1 -q n is orthogonal to G with respect to the H 1 seminorm. It follows immediately that q n+l -q n is also orthogonal to G for any l ≥ 1 and n ≥ 0:

(31) a (q n+l -q n , w) + a ⊥ (q n+l -q n , w)

= n+l-1 i=n a (q i+1 -q i , w) + a ⊥ (q i+1 -q i , w) = 0 , ∀w ∈ G. Moreover, if ∇ ⊥ q 0 = 0 than q n is orthogonal to G in the H 1 seminorm.
Proof of Proposition 6. Let (q, ū) be the stationary point of the iterative scheme. Eqs. ( 16) and ( 17) yield

-∆ q -∆ ⊥ ū = f (32) -∆ ū -ε 0 ∆ ⊥ ū = -ε∆ q + ε 0 f + ∆ q (33)
which gives ∆ ū = ε∆ q, (34) a relation that couples the parallel gradient of ū with that of q. Combining this again with (32) one obtains the initial singular perturbation problem:

- 1 ε ∆ ū -∆ ⊥ ū = f. (35)
The boundary conditions become:

n • ∇ q = -n • ∇ ⊥ ū (36) and n • A ε0 ∇ū = -(ε 0 -ε)n • ∇ q = (ε 0 -ε)n • ∇ ⊥ ū, ( 37
)
which proves that the boundary conditions for the original singular perturbation problem are recovered for the converged solution of the iterative scheme.

This problem admits a unique solution ū for ε = 0, independent on u 0 . If ε = 0 then Eqs. (32) and (35) provide the following system:

-∆ q -∆ ⊥ ū = f , -∆ ū = 0 . (38)
The second equation forces ū to belong to the space G of functions constant in the direction of anisotropy and the first equation defines the strong formulation of the limit problem [START_REF] Günter | Finite element and higher order difference formulations for modelling heat transport in magnetised plasmas[END_REF] with q acting as a Lagrange multiplier.

Proof of Theorem 2. Let us write [START_REF] Tang | An asymptotic preserving method for strongly anisotropic diffusion equations based on field line integration[END_REF] for u n with n ≥ 1, subtract it from the equation for u n+1 and choose v = u n+1 -u n . One obtains

(39) a u n+1 -u n , u n+1 -u n + ε 0 a ⊥ u n+1 -u n , u n+1 -u n = -(ε 0 -ε)a q n -q n-1 , u n+1 -u n
The Cauchy-Schwartz inequality yields

||u n+1 -u n || ε0 ≤ (ε 0 -ε)||∇ (q n -q n-1 )||. ( 40
)
It is now sufficient to prove that the sequence (q n ) n>0 converges. Thanks to Lemma 4, it follows

q n+1 -q n = A I q n -q n-1 (41)
for n ≥ 1, A I being the iteration operator defined by Eq. ( 21). The eigenvalues of A I are real and non negative (see lemma 4). Moreover, the largest eigenvalue is equal to 1 with the associated eigenfunctions belonging to G, the kernel of the operator ∆ . Thanks to lemma 5, in particular Eq. ( 29), we conclude that q n -q n-1 does not contain any non trivial function from G which concludes the proof.

Remark 7 (Non uniqueness of q). The stationary point (q, ū) of the iterative method [START_REF] Narski | Asymptotic Preserving scheme for strongly anisotropic parabolic equations for arbitrary anisotropy direction[END_REF][START_REF] Ratnani | Anisotropic Diffusion in Toroidal geometries[END_REF] is one of the solutions of the following problem: find (q, u)

∈ V × V such that a ⊥ (u, v) + a (q, v) = (f, v) , ∀v ∈ V , a (u, w) -εa (q, w) = 0 , ∀w ∈ V . ( 42
)
The above system does not admit a unique solution. Indeed, if (q, u) solves this problem, then (q + g, u) also does for any function g ∈ G. Therefore q is not uniquely determined by Eqs. (42). The TFI method consists in the choice q = q. Note that q depends on the initial value q 0 and is obtained as the stationary point of the well posed set of equations [START_REF] Narski | Asymptotic Preserving scheme for strongly anisotropic parabolic equations for arbitrary anisotropy direction[END_REF][START_REF] Ratnani | Anisotropic Diffusion in Toroidal geometries[END_REF].

2.4.

A comparison with the Micro-Macro AP-schemes. The two field iterated scheme bares some similarities with the Asymptotic Preserving scheme based on Micro-Macro decomposition (MMAPP) proposed in [START_REF] Degond | An asymptotic-preserving method for highly anisotropic elliptic equations based on a micro-macro decomposition[END_REF]. Indeed, the MMAP scheme couples Eqs. (32) and (34) but with a different choice for q and no iterative process. The weak formulation of the MMAP scheme writes:

Find (q, u) ∈ V in × V such that a ⊥ (u, v) + a (q, v) = (f, v) , ∀v ∈ V , a (u, w) -εa (q, w) = 0 , ∀w ∈ V in . V in = {v ∈ V : v = 0 on Γ in } , (43)
Γ in being the part of the boundary where b • n > 0. Note that, from this formulation, the two equations can be hardly decoupled to define an iterative process. Hence, the MMAP scheme is solved at once for the two fields (u, q). The uniqueness of q is strongly related to the assumption that all the field lines intersect the domain boundary, hence the definition of V in . The MMAP method is therefore restricted to anisotropy fields that do not contain closed lines. However, one can tackle this problem by introducing a stabilization operator [START_REF] Narski | Asymptotic Preserving scheme for strongly anisotropic parabolic equations for arbitrary anisotropy direction[END_REF][START_REF] Lozinski | Numerical analysis of an asymptotic-preserving scheme for anisotropic elliptic equations[END_REF]) yielding the weak formulation of the stabilized MMAP

Find (q, u) ∈ V × V such that a ⊥ (u, v) + a (q, v) = (f, v) , ∀v ∈ V , a (u, w) -εa (q, w) = σh k (q, w) , ∀w ∈ V . ( 44
)
where σ and k are stabilization parameters chosen in order to match the scale of the scheme approximation error. Precisely, k = 2 for P 1 -FEM and k = 3 for P 2 -FEM. The difficulty here lies in the calibration of the stabilization parameters in order not to alter the precision of the scheme and to preserve a moderate condition number of the system matrix. The conditioning of the matrix for the discrete MMAP formulation depends on 1/h 4 and 1/(σh 2+k ) for the stabilized version [START_REF] Lozinski | Numerical analysis of an asymptotic-preserving scheme for anisotropic elliptic equations[END_REF]. It is therefore ε independent.

The method here introduced is well defined for all anisotropy topologies including closed field lines. The matrix stemming from the discretization of the operator ∆ ε0 is indeed invertible regardless of the anisotropy direction b. Moreover, the condition number of the two scalar systems are not only ε independent, but it also scales more favorably, as 1/(ε 0 h 2 ), independently of the precision of the numerical method. The two field iterated method may appear similar to the stabilized MMAP scheme. However, at convergence this new formulation is completely equivalent to the original set of equations with no condition on ε 0 . This is a crucial feature that allows to choose ε 0 in a large range of values. Contrariwise, this choice is tightly constrained for the stabilization parameters and of course test case dependent. This new method hence permits to overcome the limitations of the stabilization methods for high order methods.

Numerical investigations.

3.1. Introduction. The goal of this section is to present some validation tests for the proposed method. We study the finite element formulation of the problem in different two dimensional settings, finite elements being usually a method of choice when dealing with elliptic problems. We consider two frameworks. The first one is a first order P 1 -FEM on unstructured triangular grids. The second one relies on a Cartesian rectangular grid with a second order Q 2 -FEM discretization.

Three configurations are considered: 1. A rectangular domain with open field lines and oscillating anisotropy directions; 2. A ring shaped domain with closed circular field lines; 3. A rectangular domain with both open and closed lines. The first test is performed with both P 1 (unstructured grids) and Q 2 (Cartesian Mesh) finite elements. The second test case is carried out with P 1 -FEM and the last one with Q 2 -FEM.

The iterative scheme presented herein relies on the iterative resolution of a mildly anisotropic problem for both u and q. The discretization of such problems by FEM may suffer from locking [START_REF] Babuška | Locking effects in the finite element approximation of elasticity problems[END_REF][START_REF] Babuška | On Locking and Robustness in the Finite Element Method[END_REF]. We therefore propose to analyse how the locking effect can be circumvented thanks to the choice of the numerical method as well as the value of the parameter ε 0 . This later parameter defines the strength of the anisotropy in the problem solved at each step of the iterative process. It is thus a key point in parameterizing the efficiency of the numerical method.

Let us first focus on the finite element discretization of the iterative scheme. The finite element space V h denotes either the P 1 or the Q 2 elements defined on a discretization of the domain Ω with a mesh cell of typical size h. Let the homogeneous Dirichlet boundary conditions on Γ D be enforced in the definition of V h , i.e. V h ⊂ V. A discrete formulation of the scheme reads: find (q h , ūh ) ∈ V h × V h , the stationary point of the sequence (q n+1 h , u n+1 h

α = 0 α = 2, m = 1 α = 2, m = 10
) ∈ V h × V h , n ≥ 1, solution to a (u n+1 h , v h ) + ε 0 a ⊥ (u n+1 h , v h ) = ε 0 (f, v h ) -(ε 0 -ε)a (q n h , v h ), ∀v h ∈ V h , a (q n+1 h , w h ) + ε 0 a ⊥ (q n+1 h , w h ) = (f, w h ) -a ⊥ (u n+1 h -εq n h , w h ), ∀w h ∈ V h , (45)
In all the numerical investigations conducted in the sequel, the iterative method (45) is initiated with q 0 h = u 0 h = 0. The manufactured solution method is implemented in order to define the different setups. An analytic anisotropy direction is defined by means of a vector field b. The analytic expression of the problem solution u ε is used together with that of b to compute the source term f accordingly to

f = -∆ ⊥ u ε - 1 ε ∆ u ε .
This expression is introduced in the system (45) to carry out the numerical approximation (q h , ūh ). The component ūh is thus compared against the exact analytic expression of the problem solution to evaluate the effectiveness of this new numerical method.

Test problem 1:

Open field lines with oscillating anisotropy directions. Let Ω = [0, 1] × [0, 1] be the square computational domain. Let us consider the anisotropy direction defined by

b = B |B| , B = α(2y -1) cos(mπx) + π παm(y 2 -y) sin(mπx) , ( 46
)
where m/2 is the number of oscillation periods in the computational domain and α its amplitude. For α = 0 this vector field is constant and aligned along the direction of x. When α > 0 the field oscillates in the domain. The analytic solution of the problem is given by u ε = sin πy + α(y 2 -y) cos(mπx) + ε cos (2πx) sin (πy) , (47) Three configurations will be investigated. A constant anisotropy direction aligned the x-direction. This setup is defined by α = 0. An anisotropy direction slowly varying in the computational domain, parametrized by α = 2, m = 1. Finally an anisotropy direction with fast oscillations, defined by α = 2, m = 10. For these computations, the anisotropy ratio is set to ε = 10 -15 . Therefore the only variations of the problem solution occur along the direction defined by b. The plots displayed on Fig. 1 relate the solution as well as the anisotropy direction for the configurations precised above.

For these three anisotropies, the numerical method (45) is performed on 30 iterations to define the numerical approximation ūh carried out with different values for the parameter ε 0 , on eight different meshes with h ranging from 1/10 to 1/1280. The corresponding number of mesh vertices varies from 153 for the coarsest mesh to approximately 2 • 10 6 for the most refined mesh.

P -FEM, Unstructured triangular meshes. The relative L 2 and H 1 errors are displayed on Figs. 2, 3 and 4 for the aligned, slowly and rapidly varying anisotropy directions defined above.

For the aligned anisotropy direction (α = 0, see Fig. 2) and ε 0 = 10 -1 the convergence of the iterative method in the L 2 norm is obtained after at most five iterations (for the finest mesh) and after at most three iterations in the H 1 -norm. For ε 0 = 10 -2 the convergence is even faster with two iterations being sufficient for the H 1 -norm and three for the L 2 -norm. The results are however less precise than for ε 0 = 10 -1 . Moreover, for the coarsest meshes the divergence of the iterations is observed. This is due to the locking phenomena, as explained in the next lines. This effect is even more visible with ε 0 = 10 -3 .

For the slowly variable direction of anisotropy (α = 2, m = 1, Fig. 3) the convergence is slow for ε 0 = 10 -1 . The stationary point can not be reached in 30 iterations in the L 2 -norm for the most refined meshes. For intermediate and coarse meshes the convergence is however obtained in less than 10 iterations. The locking is causing slow divergence of the numerical solution for the coarsest mesh (h = 10 -1 ). For ε 0 = 10 -2 the stationary point is reached in at most 3 iterations for both norms. Some locking effects are manifested in small augmentation of the error in course of the iterations. This is observed on the L 2 norm evolution for the coarsest mesh and for all meshes using the H 1 norm. For ε 0 = 10 -3 the stationary point is obtained in just two iterations for both norms. The precision is however worse compared to ε 0 = 10 -2 and the locking causes the relative error to blow up for two coarsest meshes. For these computations, the norm of the numerical solution is converging towards zero. This feature characterizes the locking.

In the most demanding test case with rapidly oscillating anisotropy direction (α = 2, m = 10, Figure 4) the optimal value of ε 0 is again 10 -2 : the stationary point is thus obtained after two iterations in both norms. The same convergence rate is obtained for ε 0 = 10 -3 but the numerical error is approximately ten times bigger with this setting. The locking allows accurate computations only on the finest meshes for this test case. For ε 0 = 10 -1 the convergence is very slow and the stationary point is not obtained for fine meshes in 30 iterations.

Q 2 -FEM, Cartesian meshes. The results related to these computations are gathered on Figs. 5-10. The use of Cartesian grids eliminates the locking phenomenon for the anisotropy aligned with one coordinate (see Figs 5 and6 related to α = 0). The stationary point is reached in 8 iterations for ε 0 = 10 -1 , in 4 for ε 0 = 10 -2 , 3 for ε 0 = 10 -3 and 2 for ε 0 = 10 -4 . The precision remains the same whatever the values of ε 0 for the H 1 norm (see Fig. 6). An increase of the L 2 error norm is observed for the most refined meshes and the smallest ε 0 -value (10 -4 ). One can indeed observe on Fig. 5 that the L 2 error increases when the number of cells ranges from 320 × 320 (h = 0.003125) to 640 × 640 (h = 0.0015625) and then to 1280 × 1280 (h = 0.00078125). Similar conclusions hold true for the computations carried out on the two most refined meshes with ε 0 = 10 -3 . The rate of convergence is almost independent of the mesh size. However more iterations are mandatory on the finest meshes to converge the solution to the precision accessible on these grids.

This loss of precision is explained by the conditioning of the matrix (stemming from the discretization of the operator ∆ ε0 ), which is proportional to 1/ε 0 h 2 . For the most refined meshes and the smallest values of ε 0 , the condition number of this matrix (computed by MUMPS [START_REF] Amestoy | A Fully Asynchronous Multifrontal Solver Using Distributed Dynamic Scheduling[END_REF][START_REF] Amestoy | Hybrid scheduling for the parallel solution of linear systems[END_REF]) is estimated as large as 10 9 . Therefore computing a numerical approximation with a precision larger than 10 -6 is out of reach. The condition number estimated for ε 0 = 10 -1 is of the order of 10 6 which accounts for the improved precision (10 -9 ) obtained with this value of the parameter.

For the varying anisotropy directions the second order finite elements help to prevent the locking. This is a feature documented in the literature [START_REF] Babuška | Locking effects in the finite element approximation of elasticity problems[END_REF][START_REF] Babuška | On Locking and Robustness in the Finite Element Method[END_REF]. For the slowly varying case (Figs. 7 and8) the numerical solution converges even for the coarse meshes except for the smallest value ε 0 = 10 -4 . However, even in this case, no blow up of the error is observed. For ε 0 = 10 -1 the stationary point is reached in up to 12 iterations for mesh sizes smaller than or equal to 80 × 80 (h ≤ 0.025) for the L 2 norm. For finer meshes the algorithm does not converge within 30 iterations in the L 2 norm (Fig. 7). In the H 1 norm (see Fig. 8) the convergence is obtained for mesh sizes smaller than or equal to 320×320 (h ≤ 0.003125). The best performance for intermediate and refined meshes is obtained for ε 0 -values in the range [10 -4 , 10 -3 ]. The convergence is thus obtained after three iterations only. For coarse meshes however some locking effects are still observed with the deterioration of the precision, more apparent for the smallest ε 0 -values.

When the anisotropy direction is varying rapidly in the computational domain (Figs. 9 and 10), the locking is causing the blow up of the numerical error for coarse meshes. Here also, the norm of the numerical approximation is observed to converge towards zero. For intermediate and refined meshes, the convergence is observed for all values of ε 0 . The convergence rate increases with the vanishing of ε 0 . The best numerical precision is obtained for ε 0 = 10 -3 .

Partial conclusions and comments. This first test case is intended to asses the importance of the parameter ε 0 and the robustness of the method with respect to the choice of its value.

The convergence rate of the iterative method increases with the diminishing of ε 0 . However the values of this parameter are submitted to certain constraints. For low order polynomial representations, small values of ε 0 produce the so-called locking which destroys the accuracy of the approximation on the coarsest meshes, this effect being more pronounced on unstructured grids. The locking is significantly mitigated thanks to higher order polynomial representation, as reported by the literature [START_REF] Babuška | Locking effects in the finite element approximation of elasticity problems[END_REF][START_REF] Babuška | On Locking and Robustness in the Finite Element Method[END_REF] and demonstrated by the computations carried out on Cartesian meshes (since only a P 1 -FEM method is implemented on triangular meshes, this can not be illustrated on triangular grids). The comparisons of the first and second order methods demonstrate that the locking can be avoided thanks to the use of high order discretizations. This effect remains only for the coarsest meshes with Q 2 -FEM. However this feature should be put into perspective. Indeed coarse meshes do not contain enough points to resolve accurately the variations of the anisotropy. Hence any numerical method can hardly yield acceptable results. With high order methods (Q 2 -FEM), the iterative method introduced in this paper is robust with respect to the choice of the parameter ε 0 . The convergence is secured for the values of this parameter between 10 -3 and 10 -2 for all the computations carried out in the frame of this first test case. With ε 0 = 10 -4 the convergence is obtained in less than 10 iterations for all the investigations conducted. However for computations carried on the most refined meshes with this small value of ε 0 , the precision of the method is degraded. This is due to the deterioration of the system matrix conditioning, proportional to 1/(ε 0 h 2 ). This weakness shall be overcome thanks to the use of preconditioners to offset the deterioration of the matrix condition number for large meshes. These investigations are deferred to future works. Note that the results reported in the precedent figures are related to computations carrying out anisotropy strengths as large as 10 15 . No significant differences have been observed over the range of ε-values in [10 -20 , 1] regarding the method precision, convergence properties and optimal choice of the parameter ε 0 . This is illustrated on Figure 11 presenting the L 2 and H 1 errors as functions of the anisotropy strength ε as well as the number of iterations required to convergence. These computations are related to the test problem 1, with mildly varying anisotropy (α = 2, m = 1) carried out on a mesh with 320 2 points with the Q 2 -FEM and for ε 0 = 10 -3 . The method is observed to be precise irrespectively to the anisotropy strength on the range ε ∈ [10 -15 , 1].

Comparisons with MMAP scheme. The MMAP scheme, introduced in [START_REF] Degond | An asymptotic-preserving method for highly anisotropic elliptic equations based on a micro-macro decomposition[END_REF], consists in solving the two fields (u, q) problem (42). In this system, the uniqueness of the auxiliary variable q is provided by demanding additionally that q = 0 on the part of the boundary where the field lines enter the computational domain (b • n > 0).

A comparison of the precision of the two field iterative method and the MMAP method is reported in Tab. 1. Note that the conditioning of the matrix associated with 1 Test Problem 1: Comparisons of the precision of the MMAP and the iterative method for a Q 2 -FEM discretization on a mesh with 1280 × 1280 points (h = 0.00078125). The number of iterations required to obtain the smallest relative error in L 2 and H 1 norms is reported for the different ε 0values parameterizing the iterative method ("-" meaning that the method has not converged in 30 iterations). the MMAP scheme is proportional to 1/h 4 . It is estimated as large as 10 12 for the aligned (α = 0) and rapidly varying anisotropy directions (α = 2, m = 10) and 10 10 for slowly varying directions (α = 2, m = 1) for the most refined mesh considered so far (1280 × 1280, h = 0.00078125). The large matrix conditioning deteriorates the precision of the method for the aligned case for the finest mesh and the optimal convergence rate is lost. That explains the fact that the iterative scheme is 100 more precise than the MMAP method in this configuration. For less refined meshes, the MMAP scheme and the iterative scheme with ε 0 ∼ 10 -3 yield similar precision.

α = 0 α = 2, m = 1 α = 2, m = 10 method L 2 # H 1 # L 2 # H 1 # L 2 # H 1 # MMAP 9.
The numerical efficiency of the two methods are now compared. It may seem at first glance that the iterative scheme is more time consuming that the MMAP method as it requires several resolutions of a linear system. However, the system related to the iterative scheme is twice as small and hence its resolution is faster and requires less memory in comparison to the MMAP scheme. Moreover the iteration number to convergence is small when ε 0 is close to the optimal range of values. In Tab. 2 a comparison of the computational time for both methods is proposed. The same sparse direct solver, namely the MUMPS package [START_REF] Amestoy | A Fully Asynchronous Multifrontal Solver Using Distributed Dynamic Scheduling[END_REF][START_REF] Amestoy | Hybrid scheduling for the parallel solution of linear systems[END_REF] is used to implement the LU matrix factorization and solve the linear systems involved in any method. These results show that the MMAP method is approximately twice as fast on coarse and intermediate meshes. On the 640×640 mesh it is the iterative scheme that performs better. Finally, for the finest mesh (1280 × 1280), the difference is clearly in favour of the iterative scheme which turns to be seven times faster. The computational cost required for the resolution of a linear system twice as large explains the poor efficiency of the MMAP compared to the iterative method for the largest mesh. Solving few times a linear system with a small size is more efficient than solving once a large linear system. Comparison of the runtime of the iterative scheme (total runtime for 30 iterations, time per iteration, number of iterations for convergence and time to convergence) with the runtime obtained for the MMAP scheme for Q 2 -FEM on different mesh. Runtimes obtained on the MacBook Pro laptop equipped with a 3.1 GHz Intel Core i7 dual core processor, 16GB of RAM and a Solid State Drive. The code is written in fortran compiled with gfortran-5.4.0 with -Ofast -march=corei7 optimization flags.

3.3.

Test problem 2 -diffusion in a ring. This test case reproduces the framework proposed in [START_REF] Crouseilles | Comparison of Numerical Solvers for Anisotropic Diffusion Equations Arising in Plasma Physics[END_REF] and [START_REF] Ratnani | Anisotropic Diffusion in Toroidal geometries[END_REF] investigating anisotropic diffusion problems in a Torus. It consists in simulating the diffusion in a circular domain, a context representative of magetized plasma simulation for Tokamaks. The computational domain is defined by Ω = {(x, y) ∈ R × R|0.25 ≤ x 2 + y 2 ≤ 1} and the anisotropy direction is given by the field b provided in polar coordinates (r, θ):

b = cos θ -sin θ . (48)
The analytic solution of the original problem, as represented on Fig. 12, is given by

u ε = -sin(2πr) + ε sin(2πr) cos θ . ( 49
)
These simulations are only performed on unstructured meshes (triangles and P 1 -FEM) with ε = 10 -15 defining a severe anisotropy. It is important to notice, that standard discretization of this problem, although much more elaborated than the one implemented herein (see for instance [START_REF] Crouseilles | Comparison of Numerical Solvers for Anisotropic Diffusion Equations Arising in Plasma Physics[END_REF][START_REF] Ratnani | Anisotropic Diffusion in Toroidal geometries[END_REF][START_REF] Sharma | A fast semi-implicit method for anisotropic diffusion[END_REF]) cannot handle anisotropy strengthes ε -1 larger than ∼ 10 4 , this ratio being limited to 10 -2 in [START_REF] Crouseilles | Comparison of Numerical Solvers for Anisotropic Diffusion Equations Arising in Plasma Physics[END_REF][START_REF] Sharma | A fast semi-implicit method for anisotropic diffusion[END_REF]. It is important also to point out that the elliptic problem addressed in the present paper is much more demanding, from the numerical point of view, than the diffusion problem considered by other authors. Indeed the discretization of the time derivative of the solution introduces a mass matrix offsetting partially the anisotropy. This effect is more significant The source term of the problem is analytically computed according to the preceding definitions of the anisotropy direction and solutions, in order to implement the manufactured solution technique. However, the magnetic field data are discrete and carried out on the same mesh. In contrast the method proposed in [START_REF] Wang | Uniformly Convergent Scheme for Strongly Anisotropic Diffusion Equations with Closed Field Lines[END_REF] relies on an analytic magnetic field in order to perform integrations along the field lines with a resolution much greater than that of the mesh. The numerical convergence of the iterative scheme for intermediate and refined meshes and values of ε 0 equal to 10 -3 and 10 -4 is presented in Figs. 15 and16. With the largest value of ε 0 the convergence is very slow and for coarse meshes the locking prevents the convergence. Even on fine meshes the scheme has not converged in 30 iterations in both slowly and rapidly oscilating variants. We did not observe any significant difference in the convergence speed for both setups. Except for the two coarsest meshes, the convergence is observed for the computations carried out with the smallest value of ε 0 , with a rate hardly dependent of the mesh size.

Conclusions.

In this paper a new Asymptotic-Preserving scheme is introduced for the efficient resolution of anisotropic elliptic equations. This method consists in iterating the resolution of two one field problems which require the solution of the same linear system. This system is issued from the discretization of a mildly anisotropic problem, parameterized by a numerical parameter ε 0 ε, where ε -1 is the strength of the anisotropy. The advantages of this new scheme are three folds. First the method can address any topology of anisotropies including closed field lines with no specific requirement on the grid carrying out either the unknown or the magnetic field. Second, the condition number of the linear systems solved for the iterated method scales better than that of other asymptotic-preserving (Micro-Macro) methods. Third, the computational efficiency of the method may be substantially improved with respect to these later methods. This is already demonstrated for large meshes in two dimensional frameworks. More substantial gains can be anticipated for three dimensional computations since the linear systems at hand are issued from classical elliptic problems for which very efficient solvers can be used. This issue will be investigated in subsequent works. The method already appears to be robust with respect to the choice of ε 0 and do not suffer from the locking effect provided that high order methods and meshes resolving the magnetic field variations are used. The convergence of the iterations is improved for small ε 0 -values, however at the price of a deterioration of the matrix conditioning. Future works will also be devoted to the construction of preconditioners, in order to offset the deterioration of the matrix conditioning when increasing the anisotropy of the inner problems. The extension of the method to anisotropies with varying strengths is also an issue that will be investigated in subsequent works. the project AAPFR17 2TTB "APREMADAD: Asymptotic-PREserving methods for MAgnetized Plasma simulations : Drift fluid limit and Anisotropic Diffusion problems" (2017-2018) is also acknowledged.

2. 1 .

 1 Introduction and notations. Let b ∈ (C ∞ ( Ω)) d be a smooth vector field in a domain Ω ⊂ R d , with d = 2, 3 and |b(x)| = 1 for all x ∈ Ω. Let us also decompose the boundary Γ = ∂Ω into two parts: Γ D parallel to b and its complement Γ N :

Fig. 1 .

 1 Fig. 1. Test problem 1: Exact solution for three sets of parameters defining the anositopry directions.

Fig. 2 .

 2 Fig. 2. Test Problem 1 (P 1 -FEM, unstructured mesh): Relative L 2 (left) and H 1 (right) errors as functions of the iteration number, for an anisotropy direction aligned with one coordinate (α = 0) and different ε 0 -values.

Fig. 3 .

 3 Fig. 3. Test Problem 1 (P 1 -FEM, unstructured mesh): Relative L 2 (left) and H 1 (right) errors as functions of the iteration number, for a slowly varying anisotropy direction (α = 2, m = 1) and for different ε 0 -values.

Fig. 4 .

 4 Fig. 4. Test Problem 1 (P 1 -FEM, unstructured mesh): Relative L 2 (left) and H 1 (right) errors as functions of the iteration number, for a rapidly varying anisotropy direction (α = 2, m = 10) for different ε 0 -values.

  ε 0 = 10 -4

Fig. 5 .

 5 Fig. 5. Test Problem 1 (Q 2 -FEM, Cartesian grid): Relative L 2 error as a function of the iteration number for an anisotropy direction aligned with one coordinate (α = 0) and for different ε 0 -values.

  ε 0 = 10 -4

Fig. 6 .

 6 Fig. 6. Test Problem 1 (Q 2 -FEM, Cartesian grid): Relative H 1 error as a function of the iteration number for an anisotropy direction aligned with one coordinate (α = 0) and for different ε 0 -values.

  ε 0 = 10 -4

Fig. 7 .

 7 Fig. 7. Test Problem 1 (Q 2 -FEM, Cartesian grid): Relative L 2 error for a slowly varying anisotropy direction (α = 2, m = 1) and different ε 0 -values.

- 4 Fig. 8 .

 48 Fig. 8. Test Problem 1 (Q 2 -FEM, Cartesian grid): Relative H 1 error for a slowly varying anisotropy direction (α = 2, m = 1) and different ε 0 -values.

  ε 0 = 10 -4

Fig. 9 .

 9 Fig. 9. Test Problem 1 (Q 2 -FEM, Cartesian grid): Relative L 2 error for a rapidly varying anisotropy direction (α = 2, m = 10) and different ε 0 -values.

  ε 0 = 10 -4

Fig. 10 .

 10 Fig. 10. Test Problem 1 (Q 2 -FEM, Cartesian grid): Relative H 1 error for a rapidly varying anisotropy direction (α = 2, m = 10) and different ε 0 -values.

Fig. 11 .

 11 Fig. 11. Test Problem 1 (Q 2 -FEM, cartesian mesh): Relative L 2 (left) and H 1 (center) errors as well as the number of iterations to convergence in the L 2 norm (left) as functions of the anisotropy strength ε, for a slowly varying anisotropy direction (α = 2, m = 1) for ε 0 = 10 -3 and h = 1/320.

Fig. 13 .

 13 Fig. 13. Test problem 2: Relative L 2 (left) and H 1 (right) error norms for different values of ε 0 and a P 1 -FEM on different mesh resolutions.

Fig. 14 .

 14 Fig. 14. Test problem 3: Exact solutions as defined by Eqs. (51) (left) and (52) (right) and anisotropy direction.

Fig. 15 .

 15 Fig. 15. Test problem 3: Relative L 2 (left) and H 1 (right) errors as functions of the number of iterations for the slowly oscillating solution carried out with a Q 2 -FEM on several meshes and different values of ε 0 .

Fig. 16 .

 16 Fig. 16. Test problem 3: Relative L 2 (left) and H 1 (right) errors as functions of the number of iterations for the rapidly oscillating solution carried out with a Q 2 -FEM on several meshes and different values of ε 0 .

  Table

		68 10 -8	8.52 10 -5	1.47 10 -9		1.46 10 -6		4.31 -7	1.36 10 -4
	ε 0 = 10 -1	6.85 10 -10 8 8.98 10 -7 6	4.11 10 -6	-	2.86 10 -5	-	1.72 -3	-4.84 10 -3	-
	ε 0 = 10 -2	9.36 10 -10 5 8.98 10 -7 3	6.28 10 -8	-	1.74 10 -6	-	1.19 -5	-1.52 10 -4	-
	ε 0 = 10 -3	7.11 10 -9	3 8.98 10 -7 3	1.23 10 -9 12 1.42 10 -6 4	1.81 -6 5 1.36 10 -4 3
	ε 0 = 10 -4	1.23 10 -7	2 9.07 10 -7 2	1.74 10 -9 3 1.43 10 -6 2	1.78 -5 2 1.38 10 -4 2

Table 2

 2

Acknowledgements. This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. This work has been supported by the french "Agence Nationale pour la Recherche (ANR)" in the frame of the contract ANR-11-MONU-009-01 "MOONRISE: MOdels, Oscillations and NumeRIcal SchEmes" (2015-2019) as well as the "labex CIMI" (International Centre for Mathematics and Computer Science in Toulouse) in the frame of the project "SCANISO: SCalable solvers for ANISOtropic equations arising in magnetized plasma simulations" (2017-2019). Support from the "Fédération de Fusion pour la Recherche par Confinement Magnétique" (FrFCM) in the frame of

than the time step values are small. This artefact is not present in the system at hand in the present work, addressing the stationary problem. This guarantees that the numerical parameters can be set accordingly to the physics of interest rather than to prevent the deterioration of the matrix conditioning. The convergence results of the two field iterated method are presented in Fig. 13.

The scheme behaves well for this test case too. The conclusions drawn from the preceding investigations hold true for this setup. The solution converges rapidly to the stationary point for ε 0 = 10 -2 and ε 0 = 10 -3 . For ε 0 = 10 -1 the convergence is very slow and the stationary point is reached for the coarsest meshes only. The best precision is obtained for ε 0 = 10 -2 .

3.4. Test problem 3 -magnetic islands. The last test case is also related to the physics describing hot plasmas in Tokamaks. The main difficulty of this test case is the presence of two so-called magnetic islands. They consist of closed magnetic field lines in some specific regions of the domain. Some of the magnetic field lines are open and reconnect the boundaries of the domain, the other being closed. In the sequel, the typical size of these structures will be parametrized by a (in our simulations a = 0.05). The computational domain is square Ω The first solution is mildly oscillating in the domain and while the second defines a highly oscillatory solution, which is challenging for a numerical method to capture. The analytical solutions as well as the anisotropy direction are presented on Fig. 14.