Asymptotic Preserving scheme for strongly anisotropic parabolic equations for arbitrary anisotropy direction

Computer Physics Communications 185 (2014) 3189–3203

Contents lists available at ScienceDirect

Computer Physics Communications journal homepage: www.elsevier.com/locate/cpc

Asymptotic Preserving scheme for strongly anisotropic parabolic equations for arbitrary anisotropy direction Jacek Narski a , Maurizio Ottaviani b,∗ a

Université de Toulouse, Institut de Mathématiques de Toulouse, 118 route de Narbonne, F-31062 Toulouse, France

b

CEA, IRFM, F-13108 Saint-Paul-lez-Durance, France

article

info

Article history: Received 30 May 2013 Received in revised form 21 August 2014 Accepted 25 August 2014 Available online 4 September 2014 Keywords: Anisotropic parabolic equation Ill-conditioned problem Singular perturbation model Limit model Asymptotic Preserving scheme Magnetic island

abstract This paper deals with the numerical study of a strongly anisotropic heat equation. The use of standard schemes in this situation leads to poor results, due to high anisotropy. Furthermore, the recently proposed Asymptotic-Preserving method (Lozinski et al., 2012) allows one to perform simulations regardless of the anisotropy strength but its application is limited to the case where the anisotropy direction is given by a field whose lines are all open. In this paper we introduce a new Asymptotic-Preserving method, which overcomes those limitations without any loss of precision or increase in computational costs. The convergence of the method is shown to be independent of the anisotropy parameter 0 < ε < 1 for fixed coarse Cartesian grids, and for variable anisotropy directions. The context of this work is magnetically confined fusion plasmas. © 2014 Elsevier B.V. All rights reserved.

1. Introduction This work deals with the efficient numerical treatment of heat transport in a strongly anisotropic medium. In particular, we address models of magnetized plasma with magnetic field perturbations, such as those produced by tearing modes and magnetic islands. In classical transport theory of strongly magnetized plasmas, the ratio of the parallel (χ∥ ) to the perpendicular (χ⊥ ) heat conductivity of a given species (electrons or ions) scales like (Ωc τc )2 , where Ωc is the cyclotron frequency (the rotation frequency around the field lines) and τc the collision frequency. This product is very large; in the case of electrons, it can be of the order of 1010 in the weakly collisional plasmas typical of thermonuclear reactors such as ITER [1]. In plasma theory and simulations, one usually considers a reference (primary) magnetic field. It is usually, but not necessarily, the solution of a plasma equilibrium equation expressing the balance between electromagnetic and pressure forces. Often, the primary magnetic field has a symmetry, such as the axisymmetry of a tokamak device [2]. In three dimensions, a symmetry implies the existence of magnetic surfaces on which field lines lie. Plasma transport across these surfaces depends on the perpendicular transport coefficients, and since they are small, good plasma confinement would be expected when these surfaces are closed. Conversely, transport on a magnetic surface, being controlled by parallel transport coefficients, is comparatively fast. As a result, plasma states are characterized by almost constant values of the various physical quantities (density, temperatures, etc.) on a given magnetic surface, with gradients directed mainly across such surfaces. However, it turns out that in many applications the symmetry of the primary magnetic field configuration is broken either by dynamical effects, often the consequence of instabilities, or by external perturbations. As a result, closed magnetic surfaces can be broken, leading to increased transport. This is the case with magnetic islands [3]. They can be seen as secondary magnetic structures occurring when a perturbation with a small component of the magnetic field pointing in the direction perpendicular to the original magnetic surfaces produces another set of magnetic surfaces with distinct topology.

∗

Corresponding author. Tel.: +33 442252234. E-mail address: [email protected] (M. Ottaviani).

http://dx.doi.org/10.1016/j.cpc.2014.08.018 0010-4655/© 2014 Elsevier B.V. All rights reserved.

3190

J. Narski, M. Ottaviani / Computer Physics Communications 185 (2014) 3189–3203

To be more precise, let us consider perhaps the simplest interesting case in two dimensions (2D), in some rectangular domain with Cartesian coordinates (x, y). Assume a primary field By (x) = x. The original magnetic surfaces (reduced to lines in this 2D case) are given by x = const. Transport in the x direction is limited by the small perpendicular conductivity, whereas transport along y is fast. A quantity like the temperature T would settle very quickly into a state such that T = T (x), some function of x. Add now a small perturbation with a component pointing along x, such that Bx (y) = Acos(y), with A a (usually small) constant. One can see that, regardless the value of A, the domain is split into three regions of which the one comprised within the curve x2 = 2A[1 + sin(y)] (called the separatrix) has closed field lines. This region is called a ‘‘magnetic island’’. Outside the island the field lines are only bent, whereas inside they form closed loops which can be thought of as the result of reconnecting two former field lines initially on either side of the x = 0 line and with magnetic field vectors pointing in the opposite directions. Hence the name of ‘‘magnetic reconnection’’ given to this process. Inside the island, and along each closed field line, a quantity like the temperature would be roughly constant, as a consequence of the strong parallel conductivity. In steady state one expects that the temperature would be roughly constant across the whole island, unless energy is fed inside the island. Thus islands act as an effective short circuit for radial transport and they are therefore often unwanted effects in confinement devices. Theories of the formation of magnetic islands rely on various ingredients. In the regime where tearing modes (TM) [4] are linearly unstable, magnetic islands are the result of TM evolution [5] and saturation [6,7]. When, however, TMs are stable, magnetic islands can still occur through a mechanism of self-sustainment [8]. In this regime, a key element of the island dynamics is the competition between the parallel and the perpendicular heat fluxes, depending in particular on the ratio χ∥ /χ⊥ (≡1/ε ), which may ultimately determine whether the island grows or is suppressed [9]. One can then see the interest of studying the heat (or diffusion) equation with very strong anisotropy, whether in the context of magnetic island dynamics or of other applications. This equation can occur in numerical simulations as a standalone problem, or as part of a more complex problem. In the latter case it would occur as a consequence of splitting the full problem in sub-problems for more effective coding. For a given quantity u, representing for example density or temperature, the heat equation studied in this paper can be written as 1

∂t u − ∇∥ · (A∥ ∇∥ u) − ∇⊥ · (A⊥ ∇⊥ u) = 0, ε

(1)

where ∇∥ and ∇⊥ denote the gradient in the direction parallel (respectively perpendicular) to the magnetic field. Depending on boundary conditions, this problem can become ill-posed in the limit of ε → 0. Conventional numerical methods usually fail to give accurate results with reasonable computational resources for realistic physical parameters. Indeed, when 1/ε is of the order of 1010 the problem becomes extremely anisotropic and a standard discretization leads to very badly conditioned linear systems with condition number proportional to 1/(ε h2 ) (h being the spatial discretization step). It is therefore important to develop a numerical scheme that can address the problem and give accurate results independently of the value of ε . The original motivation of this work comes from fusion plasma physics, but similar anisotropic problems are encountered in many other fields. For example one can mention image processing [10,11], transport modelling in fractured geological structures [12] or semiconductor modelling [13]. Numerical resolution of strongly anisotropic problems has been addressed by many authors. For example, suitable coordinates are often used in the context of plasma simulations [14–16]. However, this approach can be difficult to implement, especially when the magnetic field is variable in time. This is why it is preferable to choose a method which does not require mesh or coordinate adaptation, like in [17]. Another approach relies on numerical schemes specially developed for anisotropic problems. Finite difference schemes were investigated in [18–20]. A high order finite element method was proposed in [21]. Multigrid methods [22] can sometimes be beneficial. These methods are usually efficient for a selected range of ε but do not behave well in the limit ε → 0. A different way to overcome this difficulty (adopted in this paper) is to apply the so called Asymptotic Preserving scheme introduced first in [23] to deal with singularly perturbed kinetic models. The idea is to reformulate the initial problem into an equivalent form, which remains well-posed in the limit of infinite anisotropy strength. A similar reformulation was applied to the anisotropic stationary diffusion equation in [24] and then to the nonlinear anisotropic heat equation in [25,26]. The method presented therein is based on the introduction of an auxiliary variable, which serves to eliminate from the equation the dominant part, i.e. the one multiplied by 1/ε . The choice of the auxiliary variable presented in those papers allows one to solve the problem regardless of the anisotropy strength but imposes serious limitations on the magnetic field. In particular, the case of magnetic islands cannot be treated by those schemes. In this paper we propose a new method which overcomes this limitation. The difference between this new scheme and the one presented in [26] lies in the choice of the auxiliary variable q. In the previous work q is defined by a parallel gradient of u: ∇∥ q = 1ε ∇∥ u. In order to ensure uniqueness, the value of q is set to 0 on the part of the boundary where field lines enter the domain. This clearly limits the application of the method only to the domains with open field lines. Indeed, if a closed line is present, then q is defined up to a function constant along this closed line. In order to overcome this drawback we propose an alternative approach. The auxiliary function still satisfies the parallel gradient relation. However this time it is required to be in the space of functions with zero average along all field lines. This requirement provides uniqueness of the auxiliary variable independently of the magnetic field topology and thus clearly resolves the limitations of the previous method. However, in the general case of arbitrary magnetic fields, this kind of space is far from being easy to discretize. This problem is overcome with a penalty stabilization technique well known in the context of Stokes equations. The idea is to relax the parallel gradient relation between u and q with a suitable stabilization term. This term ensures that the auxiliary variable has zero average along the field lines at the price of some (small) additional error, which does not alter the convergence rate. The plan of the article is as follows. In Section 2 the mathematical problem is presented. Section 3 is devoted to the description of the numerical method. The numerical tests with known analytic solutions and the application to the problem of transport in a magnetic island are presented in Section 4.

J. Narski, M. Ottaviani / Computer Physics Communications 185 (2014) 3189–3203

3191

2. Description of the mathematical problem We are interested in the solution of an anisotropic, two or three dimensional heat problem defined on a domain Ω . Let the anisotropy direction be given by a smooth and normalized vector field b, |b| = 1 and let the computational domain Ω be a bounded and sufficiently smooth two or three dimensional subset of Rd with d = 2, 3. The domain Ω is equipped with a boundary Γ , which is decomposed according to the boundary conditions into two parts: ΓD and ΓN = ∂ Ω \ ΓD with the Dirichlet and Neumann boundary condition imposed respectively. The boundaries ΓN and ΓD are not linked with the anisotropy direction. In this way one is able to deal with a wide range of physical problems. It is convenient to decompose vectors v ∈ Rd , gradients ∇φ , with φ(x) a scalar function, and divergences ∇ · v , with v(x) a vector field, into a part parallel to the anisotropy direction b and a part perpendicular to it. These parts are defined as follows:

v∥ := (v · b)b, ∇∥ φ := (b · ∇φ)b, ∇∥ · v := ∇ · v∥ ,

v⊥ := (Id − b ⊗ b)v, ∇⊥ φ := (Id − b ⊗ b)∇φ, ∇⊥ · v := ∇ · v⊥ ,

such that v = v∥ + v⊥ , such that ∇φ = ∇∥ φ + ∇⊥ φ, such that ∇ · v = ∇∥ · v + ∇⊥ · v,

where ⊗ denotes the vector tensor product. The mathematical problem, we are interested in, reads: find the particle temperature u(t , x), solution of the evolution equation

 1  ∂t u − ∇∥ · (A∥ ∇∥ u) − ∇⊥ · (A⊥ ∇⊥ u) = 0, in [0, T ] × Ω ,    ε 1 n · ( A∥ ∇∥ u(t , ·)) + n⊥ · (A⊥ ∇⊥ u(t , ·)) = gN (t , ·), on [0, T ] × ΓN , (PH ) ∥  ε  u(t , ·) = gD (t , ·), on [0, T ] × ΓD ,   u(0, ·) = u0 (·), in Ω . The diffusion coefficients A∥ and A⊥ are bounded and of the same order of magnitude, satisfying 0 < A0 ≤ A∥ (x) ≤ A1 ,

for almost all x ∈ Ω ,

A0 ∥v∥ ≤ v A⊥ (x)v ≤ A1 ∥v∥2 , 2

t

∀v ∈ Rd and for almost all x ∈ Ω ,

with 0 < A0 < A1 some positive constants. The anisotropy of the problem is characterized by a parameter ε , which can be very small and causes substantial difficulties in the limit ε → 0. Indeed, putting formally ε = 0 in (PH) leads to the following reduced problem

 −∇ · (A∥ ∇∥ u) = 0, in [0, T ] × Ω ,   ∥ n∥ · (A∥ ∇∥ u(t , ·)) = 0, on [0, T ] × ΓN ,  u(t , ·) = gD0(t , ·), on [0, T ] × ΓD , u(0, ·) = u (·), in Ω . This problem may be ill-posed, depending on the boundary conditions and on the anisotropy field b. For example, when some field lines of b are closed in Ω the system would admit infinitely many solutions, as any function constant along the closed lines of b (meaning ∇∥ u ≡ 0) and satisfying the boundary conditions would solve the reduced problem. The same problem occurs when the field lines are open but do not pass through a boundary supplied with the Dirichlet conditions. This argument applies also to the case, where periodic (instead of Neumann) boundary conditions are imposed on ΓN . This is the case of numerical simulations related to tokamak fusion plasmas, where the computational domain is topologically equivalent to a torus. The numerical discretization of the original (PH) problem in the limit ε → 0 can therefore lead to a very badly conditioned linear system. Indeed, the condition number is proportional to 1/ε. The goal of this work is to introduce a new numerical scheme which permits to solve the original singular perturbation problem (PH) independently of the value of ε by the use of easy-to-implement numerical tools without the need to adapt the mesh to the anisotropy strength/direction. The method presented here generalizes the numerical scheme introduced in [26]. It removes the limitations of the cited numerical scheme while keeping all the advantages. It allows solving the initial singular perturbation problem (PH) accurately on a simple Cartesian grid, which does not need to be specially tailored to suit the topology of the anisotropy field lines and whose size does not depend on the anisotropy strength. The method is developed in the framework of Asymptotic-Preserving schemes. That is to say, the method is consistent with a socalled Limit problem as ε → 0 and it is stable independently of the small parameter ε . The initial singular perturbation problem (PH) is reformulated in a suitable way. Introduction of an auxiliary variable allows one to remove the terms proportional to 1/ε from the equation. The resulting system is equivalent to the initial one and it is well-posed in the limit ε → 0. The modification introduced here allows extending the applicability of the numerical scheme proposed in [26] to the case where the field b possesses closed lines without any loss of accuracy. 3. Numerical method 3.1. Semi-discretization in space Let us consider for simplicity a homogeneous Dirichlet case gD = 0, the nonhomogeneous case being a simple extension. Let us write the variational formulation of the singular perturbation problem (PH): find u(t , ·) ∈ V := H 1 (Ω ) such that

(P ) ⟨∂t u(t , ·), v⟩

V ∗ ,V

+

1

ε

 Ω

A∥ ∇∥ u(t , ·) · ∇∥ v dx +

 Ω

A⊥ ∇⊥ u(t , ·) · ∇⊥ v dx −

 ΓN

gN (t , ·)v = 0,

∀v ∈ V

(2)

3192

J. Narski, M. Ottaviani / Computer Physics Communications 185 (2014) 3189–3203

for almost every t ∈ (0, T ). When dealing with nonhomogeneous Dirichlet boundary conditions the variational formulation is written for a function w = u − g˜ , where g˜ is a suitably chosen function such that g˜ ∈ V and g˜ |ΓD = gD and hence w|ΓD = 0. The weak formulation (2) is then obtained with additional terms including g˜ and test functions v . As already discussed in the previous section, taking the formal limit of ε → 0 leads to an ill-posed problem:

 Ω

A∥ ∇∥ u(t , ·) · ∇∥ v dx = 0

with any function belonging to the vector space of functions constant in the anisotropy direction:

G := {p ∈ V /∇∥ p = 0 in Ω } being a solution. As proposed in [26], a correct Limit problem can be established by seeking a solution in the subspace G instead of V . In this case the leading order term (containing the parallel gradient) is eliminated from the equation and we are left with the following Limit problem: find u(t , ·) ∈ G such that

(L) ⟨∂t u(t , ·), v⟩V ∗ ,V +

 Ω

A⊥ ∇⊥ u(t , ·) · ∇⊥ v dx −

 ΓN

gN (t , ·)v = 0,

∀v ∈ G

for almost every t ∈ (0, T ). An Asymptotic Preserving scheme designed for the problem (P) should give accurate results and be stable independently of ε , even in the limit ε = 0. In particular, the condition number associated with corresponding linear system should not depend on ε . That is to say, the correct Limit problem (L) has to be ‘‘hardcoded’’ into the equations. Let us briefly recall the Asymptotic-Preserving scheme introduced in [26]. The method relies on the auxiliary variable q introduced by the relation ε∇∥ q = ∇∥ u in Ω . This trick permits to get rid of the terms of order O(1/ε) from the variational formulation. The uniqueness of q is provided by setting q = 0 on a part of the boundary Γin defined by

Γin := {x ∈ Γ /b(x) · n(x) < 0}, i.e. the part of the boundary, where the field lines enter the domain. The following reformulated problem, hereafter called the AsymptoticPreserving reformulation (AP-problem), is proposed: find (u(t , ·), q(t , ·)) ∈ V × Lin , solution of

     ∂u   gN v ds = 0, ,v + (A⊥ ∇⊥ u) · ∇⊥ v dx + A∥ ∇∥ q · ∇∥ v dx −  Ω Ω ΓN V ∗ ,V  (AP )  ∂ t   A∥ ∇∥ q · ∇∥ w dx = 0, ∀w ∈ Lin ,  A∥ ∇∥ u · ∇∥ w dx − ε Ω

∀v ∈ V

(3)

Ω

where

Lin := {q ∈ L2 (Ω )/∇∥ q ∈ L2 (Ω ) and q|Γin = 0}. The AP-problem is equivalent to the original P-problem (2) for fixed ε > 0. Moreover, putting formally ε = 0 in (AP) leads to a well-posed problem

     ∂u   gN v ds = 0, ,v + (A⊥ ∇⊥ u) · ∇⊥ v dx + A∥ ∇∥ q · ∇∥ v dx −  Ω Ω ΓN V ∗ ,V (L′ )  ∂ t    A∥ ∇∥ u · ∇∥ w dx = 0, ∀w ∈ Lin ,

∀v ∈ V

Ω

which is equivalent to the correct Limit problem (L). In this case, the auxiliary variable q acts as a Lagrange multiplier forcing u to be constant along b. The drawback of this method is the choice of the space of the auxiliary variable. Imposing q|Γin = 0 provides uniqueness of a solution but limits the application of the scheme to the case where all field lines are open. Indeed, fixing a value of q on the inflow boundary does not provide uniqueness of q on field lines which do not intersect with the inflow boundary (i.e. on closed field lines). In order to overcome this restriction we propose an approach based on the requirement of zero average of q along the field lines rather than on fixing the value of q on one of the boundaries. Let us first introduce an orthogonal (in L2 sense) complement A of G in V :

A = {q ∈ V /(q, p) = 0 , ∀p ∈ G}.

(4)

It can be easily seen that the average of every function in A is equal to zero. This idea was studied in [27], where the decomposition of the solution into mean and fluctuating parts allowed one to develop an Asymptotic Preserving scheme for diffusion equation in the case of a magnetic field aligned with the coordinate system. A later work [28] explored further this decomposition and proposed a scheme for more general magnetic field geometries (excluding, however, closed field lines) at the cost of including several Lagrange multipliers. Let us now introduce an asymptotic preserving method relying on this kind of decomposition: find (u(t , ·), q(t , ·)) ∈ V × A, solution of

     ∂u   ,v + (A⊥ ∇⊥ u) · ∇⊥ v dx + A∥ ∇∥ q · ∇∥ v dx − gN v ds = 0,  Ω Ω ΓN V ∗ ,V ∂t    A∥ ∇∥ q · ∇∥ w dx = 0, ∀w ∈ A.  A∥ ∇∥ u · ∇∥ w dx − ε Ω

Ω

∀v ∈ V

(5)

J. Narski, M. Ottaviani / Computer Physics Communications 185 (2014) 3189–3203

3193

This reformulation is equivalent to the original singular perturbation problem (1). Indeed, if (u, q) is a solution of (5) then choosing v ∈ G ⊂ V in the first equation yields



∂u ,v ∂t



 + V ∗ ,V

Ω

(A⊥ ∇⊥ u) · ∇⊥ v dx −

 ΓN

gN v ds = 0,

∀v ∈ G

(6)

which is nothing else than the weak formulation of the original problem for test functions in G. Similarly, if a test function of zero average is chosen in the first equation of (5), i.e. v = w ∈ A ⊂ V , then dividing the second equation by ε and adding it to the first leads to the variational formulation of the original problem for test functions in A:

⟨∂t u(t , ·), v⟩V ∗ ,V +

 Ω

A⊥ ∇⊥ u(t , ·) · ∇⊥ v dx +

1

ε

 Ω

A∥ ∇∥ u(t , ·) · ∇∥ v dx −

 ΓN

gN (t , ·)v = 0,

∀v ∈ V .

In the general case the vector space A of functions with zero average along the b lines is difficult to discretize, as the anisotropy direction can be arbitrary, e.g. the magnetic field b can vary in time or be a result of some external simulations. This makes the application of (5) limited only to very simple magnetic field topologies. This issue can be overcome by a penalty stabilization and relaxation of the ∇∥ q = 1ε ∇∥ u relation. This procedure allows one to replace the vector space A by an easy-to-discretize one at the cost of some small additional error. This method shares some analogies with a penalty stabilization used for the Stokes problem [29], where the ∇ · u = 0 requirement is replaced by ∇ · u = h2 ∆p resulting in a coercive bilinear form associated with a finite element formulation of the problem and hence the uniqueness of the solution. The reason of stabilization is quite different in the case of the AP method since the uniqueness of the solution of (5) is already guaranteed. The goal is to develop a scheme which is easy to discretize in a general setting. A new method is obtained by the addition of a small penalization term to the second equation of (5). The APS-scheme reads: find (uα (t , ·), qα (t , ·)) ∈ V × V , solution of

 α       ∂u , v gN v ds = 0, + (A⊥ ∇⊥ uα ) · ∇⊥ v dx + A∥ ∇∥ qα · ∇∥ v dx −  Ω Ω ΓN V ∗ ,V   (APS )  ∂ t   A∥ ∇∥ qα · ∇∥ w dx − α qα w = 0, ∀w ∈ V ,  A∥ ∇∥ uα · ∇∥ w dx − ε Ω

Ω

∀v ∈ V

(7)

Ω

where α is a positive stabilization constant. Note that we are looking for qα in V , a vector space which is not linked to the anisotropy direction and is easy to discretize. Let us observe that the stabilization term ensures that the auxiliary variable qα belongs to the space A for any non-zero α . Indeed, choosing w ∈ G in the second equation of (7) yields

 Ω

qα w = 0 ∀w ∈ G.

(8)

That is to say, qα is unique. However, in the stabilized scheme the relation ∇∥ uα = ε∇∥ qα is no longer fulfilled since a perturbation proportional to α is added. If α is too big, then the stabilization procedure introduces too big an error. On the other hand, in the limit of α → 0 the solution of the APS method converges to the one of (5). The formulation (5) is therefore the limit problem of (7) with respect to α . But if α is set to 0 then the uniqueness of qα is lost. This suggests that in practice one should choose α to be of the order of the truncation error so that the introduced error does not alter the convergence rate of the scheme and the solution remains unique. We postpone a more detailed theoretical analysis of the proposed method to a forthcoming paper. A similar method with stabilization by a diffusion matrix instead of a mass matrix was recently proposed in the context of a a posteriori error indicator and mesh adaptation for strongly anisotropic elliptic equations in [30]. Let us now choose a polygonalization of the domain Ω with polygons of the diameter approximately equal to h and introduce the finite element space Vh ⊂ V . The finite element discretization of (7) writes then: find (uh , qh ) ∈ Vh × Vh , the approximation of (uα , qα ), such that

    ∂ uh   v dx + ( A ∇ u ) · ∇ v dx + A ∇ q · ∇ v dx − gN vh ds = 0, h ⊥ ⊥ h ⊥ h ∥ ∥ h ∥ h  Ω ∂t Ω Ω ΓN   (APS )h     A∥ ∇∥ uh · ∇∥ wh dx − ε A∥ ∇∥ qh · ∇∥ wh dx − hk+1 qh wh = 0, ∀wh ∈ Vh . Ω

Ω

∀vh ∈ Vh

(9)

Ω

Remark that in order to ensure convergence rate in L -norm we have put α = hk+1 , where k is the order of the finite element method. 2

3.2. Semi-discretization in time One should be extremely careful when discretizing in time the (9) scheme as not all numerical schemes conserve the AP property. In fact a chosen method should be L-stable. The notion of L-stability comes from the theory of multistep methods for ordinary differential equations (see [31]). Let us consider a linear test problem y′ = ky subject to y(0) = 1. The solution of this problem is ekt . Many multistep methods (including all Runge–Kutta methods) to solve this test problem can be written as yn+1 = φ(kτ )yn = (φ(kτ ))n y0 , with τ being a 1 time step and φ(z ) a function depending on the scheme (for example φ(z ) = 1− for the implicit Euler scheme). In the case of multistep z

Runge–Kutta methods the function φ(z ) is linked to the Butcher table by the following relation: φ(z ) = 1 + zbT (I − zA)−1 e, where bT represents the bottom line of the table, A is the matrix of coefficients aij and e is the vector of ones [31]. The function φ(z ) is known as a stability function. A method is absolutely stable (A-stable) if |φ(z )| < 1 for all z with negative real part. This implies that the numerical solution approaches 0 as t → ∞, just like the exact solution. The A-stability of the scheme does not however guarantee a good precision of the numerical results, especially when one deals with rapidly decaying functions (|k| is large). Indeed, if |φ(z )| → 1 as z → ∞ then the

3194

J. Narski, M. Ottaviani / Computer Physics Communications 185 (2014) 3189–3203

numerical solution decays very slowly, contrary to the exact solution. This leads to the concept of L-stability. If a method is A-stable and moreover |φ(z )| → 0 as z → ∞ then it is said to be L-stable. In this case the 0 is obtained after just one step for infinitely stiff problems (k → −∞). This is what should be expected from an AP scheme. Indeed, if a method is A-stable and not L-stable, then the AP property is lost for small ε and large time steps τ , as the stiffness ratio is related to 1/ε . If ε is set to zero then the exact solution is immediately constant in the direction of the anisotropy, while the numerical solution is not. Hence, the scheme is not Asymptotic Preserving under time discretization that is not L-stable. A standard first order implicit Euler scheme is L-stable. If however a higher order in time method is desired then certain simple schemes are excluded. For example a Crank–Nicholson discretization is only A-stable and not L-stable, since in this case |φ(z )| → 1 as z → ∞ and the numerical solution approaches zero very slowly and in an oscillatory way (see [26] for numerical examples). The obtained system would not be asymptotic preserving and it would give reliable results only under certain assumptions. Namely, ε should be close to one, the time step should be of the order of ε or the initial value u0 should already be a solution of the stationary equation, i.e. the parallel gradient of u0 should be proportional to ε . This is why we choose in this work to present both a first order implicit Euler scheme and a second order, L-stable Runge–Kutta method. 3.2.1. Implicit Euler scheme Let us introduce  the standard bilinear forms

(Θ , χ ) :=

Θ χ dx,  a∥ (Θ , χ ) := A∥ ∇∥ Θ · ∇∥ χ dx, Ω

Ω

a⊥ (Θ , χ ) :=

 Ω

A⊥ ∇⊥ Θ · ∇⊥ χ dx,

and write the first order, implicit Euler method in a more compact notation: Find (uhn+1 , qhn+1 ) ∈ Vh × Vh , solution of

    (un+1 , v ) + τ a (un+1 , v ) + a (qn+1 , v ) − g ( t )v ds = (unh , vh ) h ⊥ h h ∥ h h N n +1 h h (EAPS ) ΓN  a∥ (unh+1 , wh ) − ε a∥ (qnh+1 , wh ) − hk+1 (qnh+1 , wh ) = 0. 3.2.2. L-stable Runge–Kutta method Any s-stage Runge–Kutta method applied to the following problem ∂u = L(t )u + f (t ), ∂t is conveniently defined using a Butcher diagram: c1

a11

···

cs

as1 b1

··· ···

.. .

a1s

.. .

.. . .

ass bs

The method reads: for a given un , an approximation of u(tn ), un+1 is determined by un + 1 = un + τ

s 

bj (L(t + cj τ )uj + f (t + cj τ )),

j=1

where ui are solutions to the following problems: ui = un + τ

s 

aij (L(t + cj τ )uj + f (t + cj τ )).

j =1

Moreover, if bj = asj for j = 1, . . . , s then un+1 is equal to the last stage of the method: un+1 = us . In order to obtain a scheme which is second order accurate in time, we choose to implement a two stage Diagonally Implicit Runge– Kutta (DIRK) second order scheme. The scheme is developed according to the following Butcher’s diagram:

λ 1

λ

1−λ 1−λ

0

λ . λ

The method is known to be L-stable for λ = 1 − √1 . 2

The second order AP-scheme writes: find (unh+1 , qnh+1 ) ∈ Vh × Vh , solution of

  (un+1 , v ) + τ λ a (un+1 , v ) + a (qn+1 , v ) − gN (tn + λτ )vh ds = (unh , vh ) h ⊥ 1,h h ∥ 1,h h 1,h ΓN    n+1 n +1 k+1 n+1  0 a∥ u1,h , wh − εa∥ (q1,h , wh ) − h (q1,h , wh ) =  1 − λ  n+1 (RKAPS ) (un+1 , v ) + τ λ a (un+1 , v ) + a (qn+1 , v ) − gN (tn + τ )vh ds = (unh , vh ) + u1,h − unh , vh h ⊥ 2,h h ∥ 2,h h 2 ,h λ ΓN    n+1 1 k+1 n+1 a∥ u2,h , wh − ε a∥ (qn2+ , w ) − h ( q , w ) = 0 h h ,h 2 ,h 1 1 unh+1 = un2+ qnh+1 = qn2+ ,h , ,h , 1 n+1 with un1+ ,h and u2,h denoting the solutions of the first and the second stage of the Runge–Kutta method.

J. Narski, M. Ottaviani / Computer Physics Communications 185 (2014) 3189–3203

3195

4. Numerical results We are now ready to test the proposed schemes. To begin, we perform numerical experiments for a b having only open field lines. This setting allows us to compare implicit Euler-APS and DIRK-APS methods with their non penalty stabilized versions (Euler-AP and DIRK-AP) proposed in [26]. The Euler-AP and DIRK-AP schemes are obtained by discretization of the AP-problem and differ from Euler-APS and DIRK-APS in two details. Stabilization terms are absent and the vector space Vh is replaced by Lin,h . The methods are also compared to a standard implicit Euler discretization of the initial singular perturbation problem (PH), given by

(P )hτ

(

unh+1

, vh ) + τ



(

a⊥ unh+1

, vh ) +

1

ε

( ,

a∥ unh

unh+1

, vh ) −

 ΓN

gN (tn+1 )vh ds



= (unh , vh ).

Finally, the DIRK-APS scheme is applied to the case of magnetic islands, where some field lines of b are closed. But let us first introduce a finite element space that we intend to use in all numerical experiments. 4.1. Discretization Let us consider a 2D square computational domain Ω = [0, 1] × [0, 1]. We choose to perform all simulations on structured meshes. Let us introduce the Cartesian, homogeneous grid xi = i/Nx ,

0 ≤ i ≤ Nx ,

yj = j/Ny ,

0 ≤ j ≤ Ny ,

with Nx and Ny being positive even constants, corresponding to the number of discretization intervals in the x- resp. y-direction. The corresponding mesh-sizes are denoted by hx > 0 resp. hy > 0. We choose a standard Q2 finite element method (Q2 -FEM), based on the following quadratic base functions

 (x − xi−2 )(x − xi−1 )     2h2x  θxi = (xi+2 − x)(xi+1 − x)    2h2x  

x ∈ [xi−2 , xi ], x ∈ [xi , xi+2 ], else,

0

 (y − y )(y − y ) j−2 j −1    2  2h  y θyj = (yj+2 − y)(yj+1 − y)    2h2y   0

y ∈ [yj−2 , yj ], y ∈ [yj , yj+2 ], else

for even i, j and

θxi =

  (xi+1 − x)(x − xi−1 ) 

x ∈ [xi−1 , xi+1 ],

h2x

0

else,

θyj =

  (yj+1 − y)(y − yj−1 ) h2y



0

y ∈ [yj−1 , yj+1 ], else

for odd i, j. We define the space



 vh =

Wh :=



vij θxi (x) θyj (y) .

i,j

The spaces Vh and Lin,h are then defined by

Vh = Wh ,

Lin,h = {qh ∈ Vh , such that qh |Γin = 0}.

The matrix elements are computed using the 2D Gauss quadrature formula, with 3 points in the x and y direction:



1

−1



1

f (ξ , η) dξ dη =

1 

ωi ωj f (ξi , ηj ),

i,j=−1

−1



where ξ0 = η0 = 0, ξ±1 = η±1 = ±

3 , 5

ω0 = 8/9 and ω±1 = 5/9, which is exact for polynomials of degree 5. Linear systems obtained

for all methods in these numerical experiments are solved using an LU decomposition, implemented by the MUMPS library. 4.2. Numerical tests 4.2.1. Known analytical solution Let the computational domain Ω be supplied with the boundaries ΓN = {(0, y) ∪ (1, y)|y ∈ [0, 1]} and ΓD = {(x, 0) ∪ (x, 1)|x ∈ [0, 1]} with the boundary conditions gN = gD = 0. Let us now construct a numerical test case with a known solution. Finding an analytical solution for an arbitrary b-field is extremely difficult. We prefer rather to act as in the previous papers [28,26]. We start by selecting a suitable limit solution u0 . The anisotropy direction is defined by isolines of u0 . Let u0 = sin π y + β(y2 − y) cos(π x) e−t ,





where β is an amplitude of variations of b. For β = 0, the limiting solution is constant in the X -direction and represents a solution for b = (1, 0)T . For β ̸= 0 the anisotropy direction is determined by the following implication

∇ ∥ u0 = 0 ⇒ b x

∂ u0 ∂ u0 + by = 0, ∂x ∂y

3196

J. Narski, M. Ottaviani / Computer Physics Communications 185 (2014) 3189–3203

Table 1 The absolute error of u in the L2 -norm for different mesh sizes and ε = 1 or ε = 10−20 , using the singular perturbation scheme (P) and the two proposed AP-schemes for a time step of τ = 10−6 and at instant t = 10−4 . L2 -error ε = 1

h

0.1 0.05 0.025 0.0125 0.00625 0.003125

P

EAP

EAPS

RKAP

RKAPS

5.6 × 10−3 7.1 × 10−4 8.9 × 10−5 1.11 × 10−5 1.39 × 10−6 1.74 × 10−7

5.6 × 10−3 7.1 × 10−4 8.9 × 10−5 1.11 × 10−5 1.39 × 10−6 1.74 × 10−7

5.6 × 10−3 7.1 × 10−4 8.9 × 10−5 1.11 × 10−5 1.39 × 10−6 1.74 × 10−7

5.6 × 10−3 7.1 × 10−4 8.9 × 10−5 1.11 × 10−5 1.39 × 10−6 1.74 × 10−7

5.6 × 10−3 7.1 × 10−4 8.9 × 10−5 1.11 × 10−5 1.39 × 10−6 1.74 × 10−7

1.62 × 10−3 2.20 × 10−4 2.77 × 10−5 3.43 × 10−6 4.2 × 10−7 5.3 × 10−8

1.66 × 10−3 2.37 × 10−4 2.93 × 10−5 3.58 × 10−6 4.4 × 10−7 5.5 × 10−8

1.62 × 10−3 2.20 × 10−4 2.77 × 10−5 3.43 × 10−6 4.2 × 10−7 5.3 × 10−8

1.66 × 10−3 2.37 × 10−4 2.93 × 10−5 3.58 × 10−6 4.4 × 10−7 5.5 × 10−8

L2 -error ε = 10−20 6.9 × 10−1 6.9 × 10−1 6.9 × 10−1 6.9 × 10−1 6.9 × 10−1 6.9 × 10−1

0.1 0.05 0.025 0.0125 0.00625 0.003125

(a) ε = 1.

(b) ε = 10−20 .

Fig. 1. Relative L2 -errors between the exact solution uε and the computed solution for the standard scheme (P), the Euler-AP method (EAP ) the stabilized Euler-APS method (EAPS ), the DIRK-AP scheme (RKAP ) and the stabilized DIRK-APS scheme (RKAPS ) as a function of h, for ε = 1 resp. ε = 10−20 and the time step τ = 10−6 . Observe that for ε = 1 all schemes give the same precision, for ε = 10−20 the standard scheme does not work while all AP schemes give comparable accuracy.

and therefore b can be defined for example as b=

B

|B|

,

 B=

 β(2y − 1) cos(π x) + π . π β(y2 − y) sin(π x)

We set β = 1 in all our simulations so that the direction of the anisotropy is variable in the computational domain. Note that we have B ̸= 0 in the computational domain. Finally, a perturbation proportional to ε is added to u0 and a function u is obtained. In this way we ensure that u converges, as ε → 0, to the limiting solution u0 . For example u = sin π y + β(y2 − y) cos(π x) e−t + ε cos (2π x) sin (π y) e−t .





(10)

We set the initial condition u0 to be equal to u(t = 0), with u defined by (10) and add a suitable force term to the right hand side of the numerical schemes so that u is an exact solution to the problem. In this setting we expect all Asymptotic-Preserving methods to converge at the optimal rate, independently on ε . In order to validate our method and to confirm our expectations we first verify the convergence with respect to the space discretization. We choose a time step small enough that the space discretization error is much bigger than the time discretization error. Then, we perform numerical simulations for 100 time steps for different mesh sizes for all Asymptotic-Preserving schemes and the implicit Euler discretization of the (P) problem. Numerical errors are given in Table 1 and Figs. 1 and 2. The third order space convergence in the L2 -norm for large values of ε is achieved (as expected) by all methods. The stabilization procedure does not alter the accuracy for weak anisotropy. For small values of ε only the Asymptotic Preserving schemes give good numerical solutions. In this case the stabilization term causes a slight decrease in the precision of the results (by 2.9%–4.6% compared to non stabilized schemes), keeping however the order of convergence. Next, the time convergence of the methods is tested. In this case, the mesh size is chosen in such a way that the space discretization does not influence the numerical precision (at least for large time steps). We perform simulations for different time steps and for a fixed final time (t = 0.1). The results are presented in Table 2 and Fig. 3. The time discretization convergence order is confirmed. As expected, the (RKAP ) and (RKAPS ) schemes are of second order in time as long as the error due to the space discretization is smaller than the error induced by the time discretization. The (EAP ) and (EAPS ) schemes are of first order for all values of the anisotropy parameter and the standard (P)-scheme works well and is of first order where the anisotropy is weak (ε close to one). Note that the stabilization

J. Narski, M. Ottaviani / Computer Physics Communications 185 (2014) 3189–3203

(a) h = 0.1.

3197

(b) h = 0.00625.

Fig. 2. Relative L2 -errors between the exact solution uε and the computed solution for the standard scheme (P), the Euler-AP method (EAP ) the stabilized Euler-APS method (EAPS ), the DIRK-AP scheme (RKAP ) and the stabilized DIRK-APS scheme (RKAPS ) as a function of ε , for h = 0.1 resp. h = 0.00625 and the time step τ = 10−6 . Observe that the standard scheme is accurate only for small values of anisotropy (large ε ) while all AP schemes give comparable accuracy in the whole range of ε .

(a) ε = 1.

(b) ε = 10−20 .

Fig. 3. Relative L2 -errors between the exact solution uε and the computed solution with the standard scheme (P), the Euler-AP method (EAP ) the stabilized Euler-APS method (EAPS ), the DIRK-AP scheme (RKAP ) and the stabilized DIRK-APS scheme (RKAPS ) as a function of τ , and for ε = 1 resp. ε = 10−20 and a mesh with 200 × 200 points.

Table 2 The absolute error of u in the L2 -norm for different time step using the singular perturbation scheme (P) and two proposed AP-schemes for mesh size 200 × 200 at time t = 0.1.

τ 0.1 0.05 0.025 0.0125 0.00625 0.003125 0.0015625

L2 -error ε = 1 P

EAP

EAPS

RKAP

RKAPS

1.32 × 10−3 7.2 × 10−4 3.55 × 10−4 1.75 × 10−4 8.8 × 10−5 4.4 × 10−3 2.20 × 10−5

1.32 × 10−3 7.2 × 10−4 3.55 × 10−4 1.75 × 10−4 8.8 × 10−5 4.4 × 10−3 2.20 × 10−5

1.32 × 10−3 7.2 × 10−4 3.55 × 10−4 1.75 × 10−4 8.8 × 10−5 4.4 × 10−3 2.20 × 10−5

8.4 × 10−5 2.43 × 10−5 6.2 × 10−6 1.59 × 10−6 4.8 × 10−7 2.82 × 10−5 2.64 × 10−7

8.4 × 10−5 2.43 × 10−5 6.2 × 10−6 1.59 × 10−6 4.8 × 10−7 2.82 × 10−5 2.64 × 10−7

1.14 × 10−3 6.2 × 10−4 3.07 × 10−4 1.51 × 10−4 7.6 × 10−5 3.81 × 10−5 1.90 × 10−5

1.14 × 10−3 6.2 × 10−4 3.07 × 10−4 1.51 × 10−4 7.6 × 10−5 3.81 × 10−5 1.90 × 10−5

7.4 × 10−5 2.10 × 10−5 5.3 × 10−6 1.33 × 10−6 3.44 × 10−7 1.18 × 10−7 8.3 × 10−8

7.4 × 10−5 2.10 × 10−5 5.3 × 10−6 1.33 × 10−6 3.44 × 10−7 1.18 × 10−7 8.3 × 10−8

L2 -error ε = 10−20 0.1 0.05 0.025 0.0125 0.00625 0.003125 0.0015625

2.28 × 10−1 2.53 × 10−1 2.53 × 10−1 2.50 × 10−1 2.51 × 10−1 2.53 × 10−1 2.53 × 10−1

procedure does not influence the accuracy of the solution if the time discretization error is bigger than the space discretization error. One can also observe that reasonably accurate results are obtained even with relatively large time steps, especially for the Runge–Kutta schemes. This property would allow orders of magnitude gains in integration time for a given physics application with respect to non-AP schemes. Numerical tests have confirmed that all Asymptotic Preserving schemes are accurate independently of ε . The penalty stabilization procedure can introduce a very small amount of additional error in some configurations, but the order of convergence is conserved. In the following experiments, the stabilized schemes are tested on the case of b possessing closed field lines.

3198

J. Narski, M. Ottaviani / Computer Physics Communications 185 (2014) 3189–3203

Fig. 4. Magnetic island for A = 0.01.

Fig. 5. Temperature profiles (top row) and temperature gradients (bottom row) along the X axis for the non perturbed field (A = 0) case, on the left, and for the case of a stationary island (A = 0.01) present in the centre of the domain, on the right, for the Dirichlet BC.

4.3. Temperature balance in the presence of magnetic islands In this section we describe a numerical experiment related to tokamak plasmas. This numerical test case fully demonstrates the novelty of the proposed stabilized AP scheme as it can be applied in more general settings than the previous ones [24,26]. The results of this section show that it is capable of simulating the heat transfer in the presence of closed field lines in the computational domain.

J. Narski, M. Ottaviani / Computer Physics Communications 185 (2014) 3189–3203

3199

Fig. 6. Temperature profiles along the X axis for a moving island (A = 0.01, ω = 10) in the first row, x component of temperature gradient in the second row and a corresponding anisotropy field in the last row for different time steps for the Dirichlet BC.

We consider a square computational domain Ω = [−0.5, 0.5] × [−0.5, 0.5] and a field b with a perturbation consisting of a region with closed lines. This so-called magnetic island is initially localized in the centre of the domain and moving with the velocity ω. The field is given by b=

B

|B|

,

 −A2π sin(2π (y − ωt )) B= , π sin(π x) 

3200

J. Narski, M. Ottaviani / Computer Physics Communications 185 (2014) 3189–3203

Fig. 7. Temperature profiles (top row) and temperature gradients (bottom row) along the X axis in the case of a non perturbed field (A = 0), on the left, and in the case of a stationary island (A = 0.01), present in the centre of the domain, on the right, for the Neumann BC.

where A is a perturbation parameter related to the island’s width w = 4A1/2 /π . This is the largest distance between the two branches of the separatrix, the line that divides the domain into regions of open and closed field lines. The two branches meet at the X-point, the saddle point of the vector potential. The island centre, an extremum of the vector potential, is referred to as the O-point. If A = 0 the obtained field is aligned with the Y axis and points upwards (downwards) for x > 0 (x < 0). For A > 0 the magnetic island consisting of closed field lines appears in the region around x = 0. Referring to the theory of magnetic islands, this magnetic field approximates the final state reached after the tearing instability [4] has gone through its nonlinear phase [5] and eventual saturation [6,7]. The frequency describes the rotation of an island as observed in experiments and numerical simulations. This frequency brings an additional timescale to the problem. In tokamaks, it is typically of the order of the diamagnetic frequency [2], which is usually larger than the transport rate across a typical island. As a consequence, the island cannot be considered as quasi-static, as in the work of Fitzpatrick [9]. It is therefore interesting to study to what extent heat transport is affected by the frequency of an island rotating sufficiently fast. For a static island or one rotating sufficiently slowly, we expect the fast transport along the field lines to flatten the temperature profile in the island region. We choose to impose periodic boundary conditions on {(x, y) ∈ ∂ Ω |y = −0.5} ∪ {(x, y) ∈ ∂ Ω |y = 0.5} so that the domain is topologically equivalent to the surface of a torus. We supply the computational domain with two sets of conditions on remaining boundaries: 1. Dirichlet boundary conditions ΓD = {(x, y) ∈ ∂ Ω |x = −0.5} ∪ {(x, y) ∈ ∂ Ω |x = 0.5} with gD (t , ·) =



1 0

for x = −0.5 for x = 0.5,

where the temperature is exchanged with the exterior only through a ΓD boundary.

(11)

J. Narski, M. Ottaviani / Computer Physics Communications 185 (2014) 3189–3203

3201

Fig. 8. Temperature profiles along the X axis for a moving island (A = 0.01, ω = 10) in the first row and the x component of temperature gradient in the second row for different time steps for the Neumann BC.

2. Neumann and Dirichlet boundary conditions: ΓN = {(x, y) ∈ ∂ Ω |x = −0.5} and ΓD = {(x, y) ∈ ∂ Ω |x = 0.5} with gN (t , ·) = 1,

gD (t , ·) = 0,

(12)

which corresponds to the constant heating of the left side of the computational domain. In the case of boundary conditions of the first type we expect that the presence of the island should increase the gradient of the temperature outside the island region and keep the temperature constant inside the island in such a way that the total energy of the system remains unchanged. If the boundary conditions are of the second type, i.e. in the heating case, the gradient of the temperature in the non-island region should remain constant, leading to a loss of the total energy of the system. We perform simulations with fixed ε = 10−10 and A = 0.01, giving an island of a width w = 0.4/π ≈ 13%. The magnetic field lines are shown on Fig. 4. The initial conditions for both boundary types are the same and correspond to the stationary solution with no island present. That is to say, u0 (x, y) = −x + 12 . We perform our simulations on a fixed grid of 200 × 200 points for 100 time steps τ = 2.5 × 10−3 until the final time 0.25 is reached. We compare the results for a stationary island (ω = 0) and of a moving island (ω = 10) with the no island case. We are interested in the temperature profile along the X axis as well as in the total energy of the system, i.e. the integral of the temperature in the computational domain. 4.3.1. Dirichlet boundary condition on the left edge (11) The numerical results confirm our expectations. The total energy remains constant in the system. The integral of the temperature in the computational domain equals 1/2 in all three cases: a stationary island, a moving island and the unperturbed system. The temperature is constant in the island region leading to a stronger gradient outside the perturbation. The temperature profiles along the X axis are shown in Figs. 5 and 6. In the case of a moving island the width of a constant temperature region in the temperature profile along the X axis is oscillating in time as the island moves in the domain. It is interesting to note that even at the time when there are no closed field lines across the X axis the flat region can be observed. Indeed, the x component of temperature gradient is negative but close to zero in this region. Note also that the temperature gradient increases, approaching zero, near the island at its largest width, i.e. across the O-point. This effect is known as ‘‘profile flattening’’.

3202

J. Narski, M. Ottaviani / Computer Physics Communications 185 (2014) 3189–3203

Fig. 9. Comparison between temperature profiles along the X axis for the stationary and rotating islands. Profile across the island centre (Y = YC ), on the left, and away from the centre (Y = YC ± 0.5), on the right. Note that the temperature profiles for ω = 106 and ω = 105 are superposed.

4.3.2. Neumann boundary condition on the left edge (12) In the case of Neumann boundary condition imposed on the left boundary of the domain the results again are consistent with our expectation. The presence of the island reduces the total energy. The integral of the temperature drops from 0.5 to 0.44 and the maximal temperature from 1 to 0.89 for both stationary and moving islands. The temperature profiles along the X axis are shown on Figs. 7 and 8. The results are very similar to the previous test case. The difference is in the maximal temperature and in the temperature gradient in the X direction, which is now close to that in the regions far from the island. This is of course what is expected and is consistent with the heating imposed on the left boundary by means of the Neumann boundary conditions. 4.3.3. Fast rotating magnetic island with Dirichlet boundary condition on the left edge (11) In the last experiment we study the effect of the island rotation speed on the temperature profile. We expect to see a different temperature profile when the rotation is faster than the transport rate in the parallel direction. To achieve this numerically, we decrease the island width such that A = 0.000625, increase the mesh size to 500 × 500, decrease the size of the computational domain Ω = [−0.125, 0.125] × [−0.5, 0.5] and set ε = 10−3 . We also increase the time resolution and set τ = 0.25 × 10−6 and vary the rotation speed from 103 to 106 . Then we compare the temperature profiles after 10 000 time steps. For the smallest rotation velocity (ω = 103 ) no deviation from the stationary case is observed. For other velocities the temperature profile across the islands centre is only slightly affected by the rotation velocity. The most interesting effect of the velocity is visible on the profile taken away from the island. In the stationary case the profile is a straight line with a constant gradient. If the rotation velocity is big enough, the profile begins to slightly flatten near the centre (X = 0) for ω = 104 and finally looks exactly the same as for the islands centre (ω = 105 and ω = 106 ). In fact, for a sufficiently large rotation speeds the temperature profile becomes homogeneous, i.e. independent of Y . A zoom of the temperature profiles at the islands centre (Y = YC ) and at the maximal distance from the island centre Y (Y = YC + 0.5 if YC < 0 and Y = YC − 0.5 otherwise) is presented on Fig. 9. 5. Conclusion The Asymptotic-Preserving scheme presented here proves to be an efficient, general and easy to implement numerical method for solving strongly anisotropic parabolic problems. As numerical experiments show, the second order two stage Runge–Kutta scheme is of particular interest as it produces accurate results for relatively large time steps allowing a substantial reduction of computational time. The stabilization term introduced here is an important improvement to the previously proposed Asymptotic Preserving schemes. With this new method one can successfully address anisotropy topologies that are fairly common in applications. Finally, the case of magnetic islands studied in the last section brings new and interesting results showing the homogenization of the temperature profile for fast rotating islands. Acknowledgements This work has been partially supported by the ANR project BOOST (Building the future Of numerical methOdS for iTer, 2010–2014), and by the European Communities under the contract of Association between EURATOM and CEA. The views and opinions expressed herein do not necessarily reflect those of the European Commission. References [1] [2] [3] [4] [5] [6] [7]

R. Aymar, P. Barabaschi, Y. Shimomura, The iter design, Plasma Phys. Control. Fusion 44 (5) (2002) 519. J. Wesson, Tokamaks, Clarendon Press, Oxford, 1997. D. Biskamp, Nonlinear Magnetohydrodynamics, Cambridge University Press, 1993. H. Furth, J. Killeen, M.N. Rosenbluth, Finite-resistivity instabilities of a sheet pinch, Phys. Fluids 6 (1963) 459. P.H. Rutherford, Nonlinear growth of the tearing mode, Phys. Fluids 16 (1973) 1903. D.F. Escande, M. Ottaviani, Simple and rigorous solution for the nonlinear tearing mode, Phys. Lett. A 323 (2004) 278–284. F. Militello, F. Porcelli, Simple analysis of the nonlinear saturation of the tearing mode, Phys. Plasmas 11 (2004) L13–L16.

J. Narski, M. Ottaviani / Computer Physics Communications 185 (2014) 3189–3203

3203

[8] Z. Chang, J.D. Callen, E.D. Fredrickson, R.V. Budny, C.C. Hegna, K.M. McGuire, M.C. Zarnstorff, T. group, Observation of nonlinear neoclassical pressure-gradient-driven tearing modes in tftr, Phys. Rev. Lett. 74 (1995) 4663–4666. [9] R. Fitzpatrick, Helical temperature perturbations associated with tearing modes in tokamak plasmas, Phys. Plasmas 2 (3) (1995) 825–838. [10] P. Perona, J. Malik, Scale-space and edge detection using anisotropic diffusion, IEEE Trans. Pattern Anal. Mach. Intell. 12 (7) (1990) 629–639. [11] J. Weickert, Anisotropic diffusion in image processing, in: European Consortium for Mathematics in Industry. B.G. Teubner, Stuttgart, 1998. [12] B. Berkowitz, Characterizing flow and transport in fractured geological media: A review, Adv. Water Resour. 25 (8–12) (2002) 861–884. [13] T. Manku, A. Nathan, Electrical properties of silicon under nonuniform stress, J. Appl. Phys. 74 (3) (1993) 1832–1837. [14] M. Beer, S. Cowley, G. Hammett, Field-aligned coordinates for nonlinear simulations of tokamak turbulence, Phys. Plasmas 2 (7) (1995) 2687. [15] A.H. Boozer, Establishment of magnetic coordinates for a given magnetic field, Phys. Fluids 25 (3) (1982) 520–521. [16] S. Hamada, Hydromagnetic equilibria and their proper coordinates, Nucl. Fusion 2 (1962) 23–37. [17] M. Ottaviani, An alternative approach to field-aligned coordinates for plasma turbulence simulations, Phys. Lett. A 375 (15) (2011) 1677–1685. [18] S. Günter, Q. Yu, J. Krüger, K. Lackner, Modelling of heat transport in magnetised plasmas using non-aligned coordinates, J. Comput. Phys. 209 (1) (2005) 354–370. [19] P. Sharma, G.W. Hammett, Preserving monotonicity in anisotropic diffusion, J. Comput. Phys. 227 (2007) 123–142. [20] P. Sharma, G.W. Hammett, A fast semi-implicit method for anisotropic diffusion, J. Comput. Phys. 230 (12) (2011) 4899–4909. [21] S. Günter, K. Lackner, C. Tichmann, Finite element and higher order difference formulations for modelling heat transport in magnetised plasmas, J. Comput. Phys. 226 (2) (2007) 2306–2316. [22] M.W. Gee, J.J. Hu, R.S. Tuminaro, A new smoothed aggregation multigrid method for anisotropic problems, Numer. Linear Algebra Appl. 16 (1) (2009) 19–37. [23] S. Jin, Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations, SIAM J. Sci. Comput. 21 (2) (1999) 441–454. [24] P. Degond, A. Lozinski, J. Narski, C. Negulescu, An asymptotic-preserving method for highly anisotropic elliptic equations based on a micro-macro decomposition, J. Comput. Phys. 231 (7) (2012) 2724–2740. [25] A. Mentrelli, C. Negulescu, Asymptotic-preserving scheme for highly anisotropic non-linear diffusion equations, J. Comput. Phys. 231 (24) (2012) 8229–8245. [26] A. Lozinski, J. Narski, C. Negulescu, Highly anisotropic temperature balance equation and its asymptotic-preserving resolution. arXiv:1203.6739, 2012. [27] P. Degond, F. Deluzet, C. Negulescu, An asymptotic preserving scheme for strongly anisotropic elliptic problems, Multiscale Model. Simul. 8 (2) (2009/10) 645–666. [28] P. Degond, F. Deluzet, A. Lozinski, J. Narski, C. Negulescu, Duality-based asymptotic-preserving method for highly anisotropic diffusion equations, Commun. Math. Sci. 10 (1) (2012) 1–31. [29] F. Brezzi, J. Douglas Jr., Stabilized mixed methods for the Stokes problem, Numer. Math. 53 (1–2) (1988) 225–235. [30] J. Narski, Anisotropic finite elements with high aspect ratio for an asymptotic preserving method for highly anisotropic elliptic equations, submitted for publication, arXiv preprint arXiv:1302.4269, 2013. [31] E. Hairer, G. Wanner, Solving Ordinary Differential Equations. II, in: Springer Series in Computational Mathematics, vol. 14, Springer-Verlag, Berlin, 2010. Stiff and differential-algebraic problems, second revised edition, paperback.