A Flux-Coordinate Independent (FCI) approach for anisotropic systems, not based on magnetic flux coordinates, has been introduced in Hariri and Ottaviani [Comput. Phys. Commun. 184, 2419 (2013)]. In this paper, we show that the approach can tackle magnetic configurations including X-points. Using the code FENICIA, an equilibrium with a magnetic island has been used to show the robustness of the FCI approach to cases in which a magnetic separatrix is present in the system, either by design or as a consequence of instabilities. Numerical results are in good agreement with the analytic solutions of the sound-wave propagation problem. Conservation properties are verified. Finally, the critical gain of the FCI approach in situations including the magnetic separatrix with an X-point is demonstrated by a fast convergence of the code with the numerical resolution in the direction of symmetry. The results highlighted in this paper show that the FCI approach can efficiently deal with X-point geometries.

The presence of a magnetic field in plasmas is known to introduce a strong anisotropy in the system. This is commonly met in fusion and most astrophysical plasmas. Developing numerical codes that take into account the strong plasma anisotropy and make efficient use of computational resources is crucial, for instance, to simulate ITER-like tokamak plasmas. Efforts have been made in this direction that lead to introducing optimized coordinate systems, so-called field-aligned coordinates, that take advantage of this anisotropy. Two-dimensional field-aligned coordinates, based on magnetic flux variables, were first presented in Refs. 1–6. They have a fundamental drawback of not handling situations including the magnetic separatrix, due to the singularity of the field-aligned metric. The approach proposed later in Ref. 7 avoids this shortcoming, but still relies on flux coordinates, for instance (r,θ). A generalized 3-dimensional coordinate system, referred to as the Flux-Coordinate Independent (FCI) system, not based on magnetic flux coordinates, is developed in Ref. 8 and elaborated in Ref. 9. The FCI coordinate system is constructed in a way that avoids the use of flux variables to discretize the fields in the poloidal plane. As discussed in Refs. 8 and 9, it allows the parallel derivative to be constructed by tracing the magnetic field lines from one perpendicular plane to the next, and interpolating to find the desired quantity. This frees us from using flux coordinates in the perpendicular plane and allows for the simulation of complex magnetic geometries. The aim of the present work is to highlight the advantage underlying the use of the FCI field-aligned system for X-point magnetic configurations. Examples of such configurations are magnetic islands (driven by magnetohydrodynamics modes or embedded in the magnetic equilibrium such as in stellerators) and axi-symmetric X-points on the Last Closed Flux Surface (LCFS) of tokamak plasmas.

These configurations, found in fusion as well as astrophysical plasmas, have a singularity in the metric of a coordinate system that uses the magnetic flux as a coordinate. This difficulty originates in the nature of the flux function, which has a saddle point (known as X-point) in one or more locations. Using the FCI approach avoids this problem and opens the way for the numerical simulations of plasma turbulence in a tokamak geometry encompassing both closed and open magnetic field line regions. First, this is a necessary feature of codes addressing the question of the transition from the low (L) to the high (H) confinement regime in X-point geometry. It is worth reminding that ITER would operate in the H confinement mode to achieve the expected thermonuclear fusion power. Understandably, research in this direction is a hot topic and investigating turbulence in X-point geometries is crucial for predicting the L-H power threshold. Second, studies of the interaction between micro-turbulence and magnetic islands exhibit a significant interplay affecting their dynamics. The most recent investigations on this topic are reported in Refs. 17–19.

The ultimate result of this paper is showing that the FCI approach allows, in particular, not only a more natural treatment of the operations in the poloidal plane, but it also easily deals with X-point configurations and with O-points such as the magnetic axis. In fact, the expression of the parallel gradient operator, which the approach developed in Ref. 8 relies on, is valid everywhere including the X- and O-points. For this purpose, we use the code FENICIA presented in Refs. 8 and 9, where the Flux-Coordinate Independent field-aligned approach has been implemented. FENICIA is based on a modular numerical scheme specially designed to address the anisotropic transport of any set of equations belonging to the general class of models (1) in Ref. 8. Since this code is easily extensible to different meshes and coordinate systems, we use it here to demonstrate the application of the FCI coordinate system to a magnetic island in a slab version of the code. The outline of this paper is as follows. A brief review of the FCI approach is given in Sec. II. Analytic solutions of a sound-wave propagation model for an island magnetic equilibrium are presented in Sec. III. Tests are performed at the exterior and interior of the island to show the convergence of the numerical solution to the analytic one in Sec. IV. Across the separatrix, where there are not yet known analytic solutions, two conclusive numerical convergence tests are presented in Sec. V, allowing us to show the robustness of the FCI approach in X-point magnetic configurations.

In order to fully describe a given field in a tokamak plasma, one needs high resolution in any given poloidal plane, but only a small number of these planes. Field-aligned coordinates allow for a sparse mesh along the field lines. They can be, however, constructed in a way that avoids the use of flux coordinates to discretize the fields in the poloidal plane. A new FCI approach, not based on flux variables to discretize the fields in the poloidal plane, has been developed in Ref. 8. The method exploits the elongated nature of microinstability driven turbulence whose scales perpendicular to the guiding magnetic field are typically of the order of the ion gyroradius, while parallel scales are of the order of the machine size. This paper extends the results of Ref. 8 to validate the FCI approach in an X-point configuration. The approach can be easily generalized to non-axisymmetric and stochastic magnetic field configurations. Details are beyond the scope of this paper and shall be discussed in an upcoming work. In the following, we consider a class of static low-β equilibria, such that the suitably normalized axisymmetric magnetic field can be written in the form

B=b(x)+ẑ,
(1)

where we employ a three dimensional Cartesian reference system (x, y, z), such that ẑ is the direction of the magnetic axis, the main magnetic field along z is constant and normalized to unity, and b(x) is the poloidal magnetic field in the poloidal plane (x, y). The two-dimensional vector x indicates the position in this plane. The poloidal field can be written in terms of a flux function ψ(x) such that

b=×(ψẑ).
(2)

It follows that bx=yψ and by=xψ. Magnetic surfaces can be labeled by the value of ψ. A point where ψ has a saddle point is an X-point and a point where ψ has a local extremum is an O-point. Both X- and O-points are such that ψ=0. Passing through the saddle are two separatrix legs, on which ψ is constant. Both closed and open field lines can be treated. The parallel derivative operator is given by

=b·+/z=[ψ,·]+z,
(3)

where the Poisson bracket operation is defined by [ψ,·]=xψyyψx.

The central idea is to work with a grid independent of magnetic flux surfaces, not using the magnetic coordinates, but computing the parallel derivatives directly along the magnetic field lines. One then looks for a change of coordinates from the original (x, y, z) to a new set (ξα,s), such that s can be treated as a slowly varying coordinate and only the two ξα(α=1,2) carry the information on the small scales. Taking advantage of what was learned in Ref. 6, we divide the domain in a certain number of sectors centered around zk, and extending to the boundary in the (x, y) directions, with k labeling a given sector. This leads to a set of transformations of the form

{ξα=Vα(x)+Cα(x)(zzk)s=zzk,
(4)

where Vα(x) and Cα(x) are yet unknown functions. In terms of the new variables the parallel derivative is given by

=bαVβxαξβ+(zzk)bαCβxαξβ+Cβξβ+s.
(5)

The FCI approach relies on expressing the parallel gradient operator as

=As,
(6)

thus eliminating the fast-varying derivatives. Here, A is a local function. One can show that it is always possible to construct such a coordinate system (ξα,s). The sufficient condition for the existence of the above transformation breaks at the extrema and saddle points of ψ, that is, at X- and O-points, respectively, where relation (6) is not valid. However, the key point for the FCI approach is that, precisely at the X- and O- points, the parallel gradient can still be expressed as a partial derivative along a single coordinate, namely z: |Xpointz. This, in turn, allows to construct the FCI approach everywhere, including the X and O-points, based on a field line map in the original Cartesian coordinates. In practice, considering a second-order finite difference representation, the expression of the discrete parallel gradient operator in (x, y, z) can be written as

ϕϕ(x+Δx,y+Δy,z+Δz)ϕ(xΔx,yΔy,zΔz)2Δz.
(7)

So, a same point in space can be labeled by (ξ1,ξ2,s+Δs)(x+Δx,y+Δy,z+Δz), where (Δxy) are obtained by integrating along the field lines

dxdz=yψ,dydz=xψ,
(8)

from (x, y, z) = (x0, y0, zk) to (x,y,z)=(x0+Δx,y0+Δy,zk+Δz) for a given displacement Δz = Δs.

More generally, the direct calculation of at a given point is done by considering an arc of field-line passing through a point (x, y, z). One can show that there always exists a parametrization x(τ), y(τ), z(τ) to get the following system:

dxdτ=yψ=bx,dydτ=xψ=by,dzdτ=1.
(9)

For any differentiable function f(x, y, z), the derivative along the field-line with respect to τ is computed as follows:

ddτf[x(τ),y(τ),z(τ)]=dx(τ)dτxf+dy(τ)dτyf+dz(τ)dτzf=yψxfxψyf+zf,=(f)τ.
(10)

Right at the X-point, since bx= 0 and by= 0, it trivially follows from Eq. (10) that f reduces to zf. As a result, Eq. (10) is valid everywhere (provided that z can be used to parametrize unambiguously the field line, which is the case if Bz does not change sign). Its finite difference form writes

f=f(τ+Δτ)f(τΔτ)2Δτ,
(11)

where Δτ = Δz. Similarly, if toroidal coordinates are employed in the poloidal plane, (R,φ,Z), the infinitesimal displacements become: (dR,Rdφ,dZ). These must be proportional to the components of B. Thus,

dRBRdτ,RdφBφdτ,dZBZdτ,
(12)

where τ is a parameter for the position along the field line. It is convenient to fix τ such that dφ/dτ=1 (case for toroidal sectors). The following field line equations result

dRdτ=RBRBφ,dφdτ=1,dZdτ=RBZBφ.
(13)

Then, for a function f(R,φ,Z),

ddτf(R,φ,Z)=dRdτRf+dφdτφf+dZdτZf=RBφ[BRRf+BφRφf+BZZf]=RBφ(B·f).
(14)

Again, right at the X-point, B·f=Bf trivially reduces to Bφφf. At a given point φ, this leads to the following expression in finite difference form:

B·f=BφR[f(φ+Δφ)f(φΔφ)2Δφ],
(15)

where f(φ±Δφ) corresponds to the value of f(R,Z,φ) at points (R±,Z±,φ±Δφ). Here, R± and Z± are obtained by integrating the field-line Eq. (13).

We conclude that the construction of the above scheme based on (6) is well-behaved everywhere, including the X-point. The construction of such a system allows decoupling of the geometric description in the poloidal plane from the treatment of the parallel derivatives. Any coordinate system can, therefore, be used in the poloidal plane. The parallel derivative is then computed in one go, by using values at the end points of arcs of magnetic field lines. This generally requires interpolating in the poloidal plane since, usually, end points are not grid nodes. Accuracy of the interpolation operation would constrain the resolution in the poloidal planes to be adequately high, but high resolution is needed anyway, in order to keep the necessary information on the fine structure and to carry out the operations in the poloidal plane to a satisfactory accuracy. Thus, the need to use interpolation for the parallel operations does not introduce substantial constraints as shown by using the Hermite cubic interpolation in Ref. 9. In light of these findings, a new modular code named FENICIA was developed.8,9 It uses the FCI approach and solves the general class of plasma models

tL·S=E(S)+I·S,
(16)

where S is a structure of vectors representing the state of the system. E(S) is a nonlinear operator that can be treated explicitly without much penalty. L and I represent linear operators such that the reduced problem obtained by setting E(S) = 0 could be treated implicitly. The splitting of the right hand side (r.h.s.) between E and I is by no means unique and depends on the physics to be studied. As a general rule, only the physics occurring at the timescale of interest for the specific problem are treated explicitly. The main constraint with respect to a generic r.h.s. is that the implicit problem be linear. We further assume that L and I are time-independent. And boundary conditions are built into L beforehand.

As extensively discussed in Ref. 10, slowly growing tearing instabilities reconnect magnetic flux surfaces to form magnetic islands. Analytical methods were derived to examine the linear stability and radial distribution of tearing modes.11–14 The nonlinear growth of these modes was also investigated in Ref. 15 and extended to cases where current nonlinearities are important as shown in Ref. 16.

The results presented in this paper are obtained using the code FENICIA. Consider a magnetic equilibrium characterized by a magnetic island whose half radial width is equal to δx=2A, with A being a parameter. Such an equilibrium corresponds to a magnetic field given by Eqs. (1) and (2), with

ψ(x,y)=(x1)22+Acos(y).
(17)

The considered domain, for an island centered at x = 1, is xminxxmax (in practice, xmax = 1 + Δ and xmin = 1 − Δ, with 0 < Δ < 1) and πyπ, with 0z2π.

Let us consider the following sound-wave model, belonging to Ref. 16, that propagates parallel to the magnetic field

tϕ+Cu=0,tu+C(1+τ)τϕ=0.
(18)

ϕ is the electrostatic potential and u is the wave's parallel velocity. We define two dimensionless parameters: C = a/R× 1/ρ* where a is the tokamak minor radius, R is the tokamak major radius, ρ* = ρs/a is the reduced gyro-radius with ρs being the ion sound Larmor radius; and τ = 1.

The aim of the present section is to construct analytic solutions of model (18) in the magnetic equilibrium defined above, i.e., in the presence of a magnetic island. To do so, it is convenient to employ flux coordinates. It readily appears that ψ is a label of magnetic surfaces, since ψ=0.

Let us introduce the following new set of coordinates:

ρ=ψ/A=(x1)22Acosy,
(19)
η=g(ρ)0ydy[cosy+ρ]1/2.
(20)

The coordinate η is a straight-field-line angle-like variable and has still to be calculated. As detailed in Appendix, for η to be an angle ranging from −π to π, g(ρ) has to be given different expressions depending on whether ρ is bigger or smaller than 1, ρ = 1 being the radius at the separatrix of the island. There are indeed three regions, each with its own specific expression of g and η. The interior region of the island corresponds to values of ρ < 1, while ρ > 1 characterizes the exterior region of the island. We can actually distinguish two exterior regions: the left side and the right side of the island. g(ρ) in both of the exterior regions has the following expression:

g(ρ)|ρ>1=±π2(1+ρ)1/2[K(21+ρ)]1,
(21)

where the sign ± depends on whether the right side (+) or the left side (−) of the exterior region is considered (cf. the Appendix). K is the elliptic integral of the first kind

K(x)0π/2dθ(1xsin2θ)1/2.
(22)

Conversely, g(ρ) can be shown to be as follows for the interior region ρ < 1:

g(ρ)|ρ<1=π22[K(1+ρ2)]1.
(23)

To illustrate this, contour plots of ρ and η are shown in Figs. 1(a) and 1(b), where ρ defines the island geometry and η defines the inner and outer regions of the island. It can be shown that the parallel gradient takes the following expression in terms of this new set of coordinates:

=g(ρ)2Aη+z.
(24)

We consider wave-like solutions of model (18) of the form

(ϕ(ρ,η,t)u(ρ,η,t))=(ϕ0(ρ)u0(ρ))cos[mηnzω(ρ)t],
(25)

with (m, n) standing for the wave numbers in η and z, respectively, and ω(ρ) being the mode frequency. Injecting these expressions in Eq. (18) leads to the following system:

(ωC[g(ρ)2Amn]C(1+τ)τ[g(ρ)2Amn]ω)(ϕ0(ρ)u0(ρ))=0.
(26)

The eigenfrequency solution of the dispersion relation depends on ρ

ωmn±(ρ)=±C(1+ττ)1/2[g(ρ)2Amn].
(27)

The eigenvectors are then

u0(ρ)=C(1+τ)τωmn±(ρ)ϕ0(ρ).
(28)

For a given expression of ϕ0(ρ), u0 can be calculated by using Eq. (28), with ωmn±(ρ) given by Eq. (27). Analytic solutions of the model (18), provided by Eq. (25), will be compared to their numerical counterparts in Sec. IV. While the above solutions can be derived for ρ < 1 and ρ > 1, there is no obvious analytic solution to this problem at ρ=1. This issue is left for future investigation.

FIG. 1.

The new coordinates ρ (a) and η (b) as a function of the grid mesh coordinates (x, y).

FIG. 1.

The new coordinates ρ (a) and η (b) as a function of the grid mesh coordinates (x, y).

Close modal

The first test aims at verifying that the numerical simulation of sound-wave propagation outside the X-point gives results in good agreement with the analytic solutions of this problem up to the order of accuracy of the chosen numerical scheme. FENICIA is used to solve the system of Eq. (18) with a version embedding an island magnetic equilibrium. With respect to the model presented in Refs. 8, Dirichlet boundary conditions are still used in the x direction, but periodic boundary conditions are now used in the y and z directions. The results of our investigation are presented into the influence of a static axi-symmetric magnetic island, including the X- and O-points, on sound-waves. We consider a box of size 2Δ = 0.2 (normalized to the radial position of the separatrix), an island of size 4A with A = 10−3 and a mode (m, n) resonant at ρ = ρmn, i.e., m/n=1/(g(ρmn)2A). This means that k=0 at ρmn. For this value of A, we show 1/(g(ρmn)2A) as a function of ρ for both the exterior and the interior of the island as plotted in Fig. 2. For rational values of 1/(g(ρ)2A), there exists a resonant (m, n). Accordingly, the position around ρ = 1 is directly constrained by high m or high m/n values, especially when approaching the separatrix starting from the interior of the island. We initialize a perturbation of the form

ϕ(t=0)=ϕ0(ρ)cos(mηnz),
(29)

where ϕ0(ρ) is a Gaussian structure centered around ρmn and having the following form:

ϕ0(ρ)=exp[(ρρmn)2Δρ2]×(ρρmn)m×(ρρbdρmnρbd)2,
(30)

with Δρ=ρmn/m and ρbd is the value of ρ at the boundaries x={xmin,xmax} for y = 0. The small radial width of the envelope ϕ0, around the resonant surface ρmn, ensures that k1. Two cases are to be examined: the exterior of the island (ρ > 1), and the interior (ρ<1). The only difference between the two cases comes from the two different expressions of g(ρ) given by Eqs. (21) and (23). For the exterior region of the island, we consider a perturbation ϕ0(ρ) centered around ρmn = 1.25 for (m,n) = (24,1). The initial condition for a simulation of size (nx,ny,nz)=(800×800×20) is shown in Fig. 3(a). Note that time t is normalized to the Bohm timescale (which is long by a factor 1/ρ* compared to the significant turbulence time) and the time step considered for the following simulations is Δt=103. In Fig. 3(b), we show the solution at time t = 1. For clarity, a zoomed view of the initial condition and the final solution are also shown in Figs. 3(c) and 3(d). The emphasis is now on showing that the exact slab solution given by (25) is recovered at the exterior of the island, where ρ>1. In Fig. 4, we plot the relative error

(ϕnumϕexact)21/2(ϕexact)21/2,
(31)

between the exact and the numerical solution (ϕexact and ϕnum, respectively) as a function of time for increasing spatial resolutions (nx,ny,nz)={(400,400,20);(600,600,20);(800,800,20)}, where ·=02πdy02πdz1xmaxdx. The volume integral is defined over half of the domain (right side of the island). The calculations give similar results if considering the left side of the island. The numerical results show faster than linear convergence to the analytic solution, as seen in Fig. 4. The same tests are performed at the interior of the island where we consider a perturbation ϕ0(ρ) centered around ρmn = 0.58 for (m,n) = (45,1). The initial condition for a simulation of size (800 × 800 × 20) is shown in Fig. 5(a), and the final solution at t = 1 is given in Fig. 5(b). A zoomed view of the initial condition and the final solution is also shown in Figs. 5(c) and 5(d). The numerical solution is again compared to the analytic solution given by (25) for ρ < 1. This is illustrated in Fig. 6, where we plot the relative error as a function of time. We observe that the numerical solution to the sound-wave problem converges to the one derived analytically as the grid size Δx decreases. The relative error for the resolutions considered here is yet high (∼30%). We interpret this to indicate that the poloidal resolution is not sufficient to be able to resolve the m = 45 mode at the interior of the island.

FIG. 2.

Plot of m/n as a function of ρ for both the exterior and the interior of the island.

FIG. 2.

Plot of m/n as a function of ρ for both the exterior and the interior of the island.

Close modal
FIG. 3.

(a) Initial condition ϕ(t=0); at the exterior region of the island (b) numerical solution ϕ at the final time step t = 1; (c) zoomed view of the initial condition ϕ(t=0); (d) zoomed view of the solution ϕ at the final time step t = 1.

FIG. 3.

(a) Initial condition ϕ(t=0); at the exterior region of the island (b) numerical solution ϕ at the final time step t = 1; (c) zoomed view of the initial condition ϕ(t=0); (d) zoomed view of the solution ϕ at the final time step t = 1.

Close modal
FIG. 4.

Normalized RMS of the difference between the analytic solution ϕexact and the numerical solution ϕnum at the exterior of the island as a function of time given by different spatial resolutions.

FIG. 4.

Normalized RMS of the difference between the analytic solution ϕexact and the numerical solution ϕnum at the exterior of the island as a function of time given by different spatial resolutions.

Close modal
FIG. 5.

(a) Initial condition ϕ(t=0); for the interior region of the island (b) numerical solution ϕ at the final time step t = 1; (c) zoomed view of the initial condition ϕ(t=0); (d) zoomed view of the solution ϕ at the final time step t = 1.

FIG. 5.

(a) Initial condition ϕ(t=0); for the interior region of the island (b) numerical solution ϕ at the final time step t = 1; (c) zoomed view of the initial condition ϕ(t=0); (d) zoomed view of the solution ϕ at the final time step t = 1.

Close modal
FIG. 6.

Normalized RMS of the difference between the analytic solution ϕexact and the numerical solution ϕnum at the interior of the island as a function of time given by different spatial resolutions.

FIG. 6.

Normalized RMS of the difference between the analytic solution ϕexact and the numerical solution ϕnum at the interior of the island as a function of time given by different spatial resolutions.

Close modal

At the separatrix (ρ = 1), the analytic solutions (25) are no longer valid. As stated in Sec. III, we do not have an obvious analytic solution at the separatrix. Given this limitation, rather than considering an initial condition as a function of the variable η, which diverges at rho = 1, we instead consider a perturbation crossing the separatrix of the form

ϕ(t)=N0(ρ)cos(mynz),
(32)

where N0(ρ) is a Gaussian structure given by

N0(ρ)=e(ρ1)2/Δρ2,
(33)

where Δρ is a parameter set to Δρ = 0.5 here for the perturbation to cover both the left and right-hand sides of the separatrix. A series of simulations were run with Δ = 0.1, A = 10−3 and (m,n) = (5,1). Figures 7 and 8 show the evolution of the electrostatic potential as a function of time (normalized to Bohm time). The time step for this simulation is Δt = 10−3, and the box size is (200 × 200 × 20). As it is apparent from Fig. 7, modes across the separatrix gradually evolve and shear. Eddies leaving one side of the box reenter on the other side due to the periodicity of the simulation domain. The elongated nature of structures along the magnetic field lines is identified through 3D illustrations of the same simulation, shown in Fig. 8. It becomes evident that the direction parallel to the magnetic field does not coincide with the direction of symmetry (z). Because our ultimate goal was to prove the validity of the FCI approach for X-point geometries, three main tests are to be considered hereafter. The first test consists in showing the conservation of an energy-like quadratic quantity of the sound-wave model (18). For this purpose, we perform a scan in nx and in nz, the number of points in the x and z directions, respectively (in all simulations we considered nx = ny). For a fixed number of points in the z direction, nz = 20, a scan over the following pairs is done in the x and y directions: (nx,ny)={(100,100);(200,200);(400,400);(600,600)} until t = 1. The energy-like quantity E(ϕ2+(1+1/τ)1u2)/2dxdydz is conserved by Eq. (18). It is plotted in Fig. 9(a) as a function of time. Notice that t = 1 is a Bohm time that is long compared to turbulence time scales, which typically range between 1/ρ and 10/ρ with ρ being the ratio between the ion larmor radius and the system size. Figure 9(b) shows that the conservation of E improves as the grid spacing Δx decreases. For a resolution (nx, ny) = (100,100), the relative change in Energy is equal to 4.5 × 10−3 and in the well-resolved case, where (nx, ny)= (600,600), the relative change decreases to 0.5 × 10−3. When simulating turbulence, errors of the same order of magnitude are expected. For the special case where the box size is (nx, ny, nz) = (200,200,20) shown in the 2D and 3D snapshots above, we performed the run up to t = 12. Though there are fine structures of the grid size scale appearing after t = 6, the plot in Fig. 10 shows a relatively good conservation of the energy. The relative variation of energy as a function of time reaches at most 5%. Similarly, the conservation of the energy improves with increasing resolution in the z direction. A box of size (nx, ny, nz) = (400, 400, nz) is considered where nz={20,40,60,80,100}. In Fig. 11(a), the energy variation converges towards zero as Δz decreases. Furthermore, its relative change is in the worst considered case equal to 5 × 10−4 as shown in Fig. 11(b).

FIG. 7.

Time evolution of a perturbation that straddles the separatrix including the X-point.

FIG. 7.

Time evolution of a perturbation that straddles the separatrix including the X-point.

Close modal
FIG. 8.

3D time evolution of a perturbation that straddles the magnetic separatrix including the X-point.

FIG. 8.

3D time evolution of a perturbation that straddles the magnetic separatrix including the X-point.

Close modal
FIG. 9.

(a) Conservation of energy with respect to time for different spatial resolutions. (b) The relative change in energy with respect to Δx.

FIG. 9.

(a) Conservation of energy with respect to time for different spatial resolutions. (b) The relative change in energy with respect to Δx.

Close modal
FIG. 10.

The relative change in energy with respect to time.

FIG. 10.

The relative change in energy with respect to time.

Close modal
FIG. 11.

(a) Plot of the energy as a function of time with increasing resolution in the z direction and (b) The relative change of energy as a function of Δz.

FIG. 11.

(a) Plot of the energy as a function of time with increasing resolution in the z direction and (b) The relative change of energy as a function of Δz.

Close modal

The second target test intends for the demonstration of the convergence of the numerical solution when crossing the separatrix region. For this purpose, a series of simulations were performed with a box of size (nx, ny) = (400, 400) and a scan over nz = {5, 10, 20, 40, 60, 80, 100}. We start by considering the two solutions given by the simulations with nz = 5 and nz = 10. The idea is to calculate the difference between them and repeat the process over each pair of the entire set of solutions ϕ given by all the simulations. This is called the moving difference. In Fig. 12, the moving difference of the solution ϕ for nz = {5, 10, 20, 40} at the first poloidal plane iz = 1 is plotted in 2D at time t = 1. At this stage, for the three pairs given by {ϕ(nz=5),ϕ(nz=10)},{ϕ(nz=10),ϕ(nz=20)},{ϕ(nz=20),ϕ(nz=40)}, we qualitatively obtain an error that is of the order of the solution. This means that convergence is not yet reached. In Fig. 13, the moving difference is calculated for pairs of solutions given by {ϕ(nz=40),ϕ(nz=60)},{ϕ(nz=60),ϕ(nz=80)},{ϕ(nz=80),ϕ(nz=100)}. We observe that beyond nz = 40, the difference between each pair of solutions is nearly zero. This is an indication that convergence is reached, and a resolution of nz = 40 is sufficient to resolve the physics problem considered. This constitutes the main strength of the field-aligned FCI system of coordinates and validates its application to X-point configurations with an exponential convergence rate shown in Fig. 14. The average is calculated here by dividing the sum of the difference between each pair of solutions by the total number of solutions. It is, thus, a quantitative mean that allows us to conclude that the numerical solution converges as the resolution in the z direction decreases. In fact, at nz = 40, the norm is approximately 10−2 which is consistent with the fact that we used second order centered finite differences to compute the parallel gradient operator . With this result, we conclude that the numerical solution of the sound-wave propagation problem given an equilibrium with a magnetic island has converged with the expected convergence order. To finish, we do a last test showing the order of convergence of the numerical solution. In fact, for the same set of simulations, i.e., (nx, ny) = (400, 400) and nz = {5, 10, 20, 40, 60, 80, 100}, the test consists of choosing a reference case supposed to be the closest possible case to the analytic solution of the sound-wave problem. With this hypothesis, it is legitimate to calculate the average of the difference between all the simulations and the reference simulation chosen to be that with nz = 100 (Root Mean Squared of the difference). The resulting Fig. 15(b) shows that as nz tends to 100, the error tends to 0 with an estimation of the order of convergence given by the loglog plot in Fig. 15(b). The convergence is indeed fast, of the order of a = 2.6, the corresponding slope of the loglog plot. This is once again in agreement with the use of second order finite differences to compute the parallel gradient operator. In summary, even with a highly elongated island in the perpendicular plane, which make the requirements for resolution particularly difficult, the latter results show the robustness of the FCI system. At this stage, we can finally conclude that the FCI system of coordinates is quantitatively validated for X-point geometries.

FIG. 12.

Moving difference {ϕ(nz=5),ϕ(nz=10)},{ϕ(nz=10),ϕ(nz=20)},{ϕ(nz=20),ϕ(nz=40)}: the first two columns show 2D plots of the solution ϕ for two different values of nz. The third column shows a 2D plot of the difference between the two chosen solutions. Note that in this figure, the same color-bar scale is used for all the plots.

FIG. 12.

Moving difference {ϕ(nz=5),ϕ(nz=10)},{ϕ(nz=10),ϕ(nz=20)},{ϕ(nz=20),ϕ(nz=40)}: the first two columns show 2D plots of the solution ϕ for two different values of nz. The third column shows a 2D plot of the difference between the two chosen solutions. Note that in this figure, the same color-bar scale is used for all the plots.

Close modal
FIG. 13.

Moving difference {ϕ(nz=40),ϕ(nz=60)},{ϕ(nz=60),ϕ(nz=80)},{ϕ(nz=80),ϕ(nz=100)}: The first two columns show 2D plots of the solution ϕ for two different values of nz. The third column shows a 2D plot of the difference between the two chosen solutions. Note that in this figure the same color-bar scale is used for all the plots.

FIG. 13.

Moving difference {ϕ(nz=40),ϕ(nz=60)},{ϕ(nz=60),ϕ(nz=80)},{ϕ(nz=80),ϕ(nz=100)}: The first two columns show 2D plots of the solution ϕ for two different values of nz. The third column shows a 2D plot of the difference between the two chosen solutions. Note that in this figure the same color-bar scale is used for all the plots.

Close modal
FIG. 14.

Average of the norm of the moving difference showing convergence of the numerical solution in nz at an exponential rate.

FIG. 14.

Average of the norm of the moving difference showing convergence of the numerical solution in nz at an exponential rate.

Close modal
FIG. 15.

(a) Average of the difference between the numerical solution at different nz values and the solution at the reference case nz = 100 and (b) loglog plot (natural logarithm) showing the convergence rate having a slope a = 2.6.

FIG. 15.

(a) Average of the difference between the numerical solution at different nz values and the solution at the reference case nz = 100 and (b) loglog plot (natural logarithm) showing the convergence rate having a slope a = 2.6.

Close modal

Employing coordinate systems which take advantage of the physical characteristics of magnetized plasma, namely the strong anisotropy of spatial scales between the parallel and the transverse directions (with respect to the equilibrium guiding magnetic field B) allow the simulation of more complex systems with today's computational resources. A new three-dimensional Flux-Coordinate Independent system, referred to as FCI in Refs. 8 and 9, is shown to have a number of advantages over earlier approaches. It first permits more flexible coding by decoupling the magnetic geometry from the computational grid. Second, it allows for coarse grids since it is field-aligned (the idea is to follow the field lines to calculate the parallel gradient operator). Third, the FCI system copes with X-point magnetic geometries, a situation where previous approaches based on magnetic flux coordinates have a fundamental problem due to the singularity of the field-aligned metric. In the present work, we show that the FCI approach can be extended to deal with X-point configurations such as magnetic islands. For verification of the approach, a study of sound-wave propagation across the X-point of a magnetic island was carried out using the code FENICIA (cf. Refs. 8 and 9). Analytic solutions to the sound-wave problem were constructed both inside and outside the island. The numerical results show that the FCI approach converges to the derived analytic solution. A test case is considered with an initial perturbation that straddles the magnetic separatrix including the X-point. In the absence of an analytic solution, conservation properties were verified and conclusive tests performed to show the convergence of the numerical solution with respect to the number of points in the z direction (direction of symmetry). This demonstrates the capability of the approach to study X-point physics with reduced numerical resolution in the parallel direction. If the magnetic island is highly elongated in the perpendicular plane, the requirements for the numerical resolution increase. Results show the robustness of the FCI approach also in this stressed test. The FCI approach allows for a coarser grid by exploiting the anisotropy of the system and for flux-coordinate independent operations in the poloidal plane. This opens the way to the implementation of complex magnetic geometries, without difficulties encountered with X- and O-points configurations. The flexible nature of FENICIA presented in Refs. 8 and 9 allowed us to demonstrate the application of this coordinate system to a magnetic island in a slab, thus the validity of its application to X-point configurations.

F.H. acknowledges discussions with L. Villard, Y. Peysson, V. Naulin, and S. Brunner. F.H. would, in particular, like to express her appreciation to S. Cowley for his encouragement and valuable feedback on the FCI approach. Further thanks to X. Litaudon, E. Joffrin, S.-I. Itoh, and K. Itoh for advice and support benefiting this work.

This work, supported by the European Communities under the contract of Association between EURATOM and CEA, was carried out within the framework of the European Fusion Development Agreement. The views and opinions expressed herein do not necessarily reflect those of the European Commission. Part of this work was also supported by the ANR Contract No. E2T2.

We seek to define proper coordinates associated to the field lines in the magnetic island equilibrium. For that purpose, let us consider a magnetic field given by Eqs. (1) and (2), with

ψ(x,y)=(x1)22+Acos(y).
(A1)

The parallel gradient operators is expressed as

=[ψ,·]+z.
(A2)

In order to construct analytic solutions to the above model in the presence of a magnetic island, the idea is to define a coordinate system: (x,y)(ρ,y) such that ρ is a magnetic surface label, i.e., ρ=0. One would write

ρ=ψ/A=(x1)22Acosy,
(A3)
y=y,
(A4)

where ρ′ is a normalization of ψ so that ρ′= 1 refers to the separatrix. Note that in the neighborhood of y = 0 (O-point), cosy1y2/2ρ(x1)22A1+y2/2, which defines the equation of an ellipse. With this set of coordinates, the spatial derivatives in the x and y directions write

x=xρρ+xyy=xρ/A,
(A5)
y=yρρ+yyy=sinyρ+y,
(A6)

where x′= x−1. This leads to a parallel operator of the form

=xy+Asinyx+z=x(sinyρ+y)xsinyρ+z=xy+z.
(A7)

The parallel gradient operator can then be written as

ϕ=x(ρ,y)yϕ+zϕ,
(A8)

where x(ρ,y)=±(2A)(ρ+cosy)1/2. At this point, to ensure the separation of variables, we define a new system: (ρ,y)(ρ,η) such that ρ = ρ′ and

xy(2A)g(ρ)η=±(2A)(ρ+cosy)1/2y.
(A9)

We conclude from (A9) that

dη=±g(ρ)dy(ρ+cosy)1/2.
(A10)

So, the new variable η can be expressed as

η=[±g(ρ)0ydy(ρ+cosy)1/2]ρ=x2/(2A)cosy+η0±,
(A11)

with g(ρ) chosen such that for all x′:

η(x,y=π)=πη(x,y=π)=π.
(A12)

This gives

π=±g(1+x22A)0πdy[cosy+(1+x22A)]1/2x.
(A13)

So

g(ρ)=±π[0πdy(cosy+ρ)1/2]1x.
(A14)

The function η is regular everywhere, but delicate near the X-point where ρ = 1 and 0ydy/(ρ+cosy)1/2 diverges. When ρ < 1, the argument of the integral involves the term (cosy+ρ)1/2 which becomes negative at cosy=ρ. We choose to limit the integration to the upper bound yM such that cosyM=ρ. The situation is described through a schematic of the island in Fig. 16. The explicit expression of g(ρ) at the exterior of the island (i.e., ρ > 1) can be derived as follows:

0πdy(ρ+cosy)1/2=0π/22dθ(ρ+cos2θ)1/2,=0π/22dθ[ρ+(12sin2θ)]1/2,=0π/22dθ(ρ+1)1/2[12sin2θ/(1+ρ)]1/2.
(A15)

Use definitions of elliptic integrals of first kind

K(m)0π/2dθ(1msin2θ)1/2,
(A16)

to get

0πdy(ρ+cosy)1/2=2(1+ρ)1/2K(21+ρ).
(A17)

Note that K(2/(1+ρ)) diverges logarithmically as ρ1+. Finally, for ρ>1,g(ρ) can be written as

g(ρ)=π2(1+ρ)1/2[K(21+ρ)]1.
(A18)
FIG. 16.

A schematic of an island with the new variables.

FIG. 16.

A schematic of an island with the new variables.

Close modal

At the interior of the island (i.e., ρ < 1), however, the limit of integration is yM such that cosyM=ρ. Thus

η=π+[±g(ρ)0ydy(ρ+cosy)1/2]ρ=x2/(2A)cosy,
(A19)

but now

g(ρ)=π2[0yMdy(cosy+ρ)1/2]1x.
(A20)

In order to express g(ρ) as a function of K, the elliptic integral, we proceed by first rewriting the expression I=0yMdy(cosy+ρ)1/2 by change of variables. We indeed have

cosy=1(1+ρ)sin2α,
(A21)

where 0απ/2. Thus

(cosy+ρ)1/2=(1+ρ)1/2(1sin2α)1=(1+ρ)1/2(cosα)1.
(A22)

Moreover, by differentiation, we get

sinydy=2(1+ρ)sinαcosαdα.
(A23)

Also,

sin2y=1cos2y=2(1+ρ)sin2α[11+ρ2sin2α],
(A24)

yielding

siny=[2(1+ρ)]1/2sinα[11+ρ2sin2α]1/2
(A25)

and

dy=[2(1+ρ)]1/2cosα[11+ρ2sin2α]1/2dα.
(A26)

Hence,

I=20π/2dα[11+ρ2sin2α]1/2,
(A27)
=2K(1+ρ2).
(A28)

Finally, for ρ < 1, we have

g(ρ)=π23/2[K(1+ρ2)]1.
(A29)

For x′< 0, as shown in Fig. 16, η′ needs to be completed. We define

η=π[g(ρ)0ydy(ρ+cosy)1/2]ρ=x2/(2A)cosy,
(A30)

so that π/2η3π/2 when x′ < 0. And

η=2π+[g(ρ)0ydy(ρ+cosy)1/2]ρ=x2/(2A)cosy,
(A31)

for x′> 0, but y < 0 if one wants 3π/2η2π in the fourth quadrant.

1.
S.
Cowley
,
R.
Kulsrud
, and
R.
Sudan
, “
Considerations of ion-temperature-gradient-driven turbulence
,”
Phys. Fluids B
3
,
2767
(
1991
).
2.
G.
Hammett
,
M.
Beer
,
W.
Dorland
,
S.
Cowley
, and
S.
Smith
, “
Developments in the gyrofluid approach to tokamak turbulence simulations
,”
Plasma Phys. Controlled Fusion
35
(
8
),
973
(
1993
).
3.
M.
Beer
,
S.
Cowley
, and
G.
Hammett
, “
Field-aligned coordinates for nonlinear simulations of tokamak turbulence
,”
Phys. Plasmas
2
,
2687
(
1995
).
4.
A. M.
Dimits
, “
Fluid simulations of tokamak turbulence in quasiballooning coordinates
,”
Phys. Rev. E
48
,
4070
4079
(
1993
).
5.
B.
Scott
, “
Global consistency for thin flux tube treatments of toroidal geometry
,”
Phys. Plasmas
5
,
2334
(
1998
).
6.
B.
Scott
, “
Shifted metric procedure for flux tube treatments of toroidal geometry: Avoiding grid deformation
,”
Phys. Plasmas
8
,
447
(
2001
).
7.
M.
Ottaviani
, “
An alternative approach to field-aligned coordinates for plasma turbulence simulations
,”
Phys. Lett. A
375
(
15
),
1677
1685
(
2011
).
8.
F.
Hariri
and
M.
Ottaviani
, “
A flux-coordinate independent field-aligned approach to plasma turbulence simulations
,”
Comput. Phys. Commun.
184
(
11
),
2419
2429
(
2013
).
9.
F.
Hariri
, “
FENICIA: A generic plasma simulation code using a flux-independent field-aligned coordinate approach
,” Ph.D. thesis (
Aix-Marseille University
, AMU,
2013
).
10.
D.
Biskamp
,
Nonlinear Magnetohydrodynamics
(
Cambridge University Press
,
1997
), Vol. 1.
11.
H.
Furth
,
P.
Rutherford
, and
H.
Selberg
, “
Tearing mode in the cylindrical tokamak
,”
Phys. Fluids
16
(
7
),
1054
1063
(
1973
).
12.
S. C.
Cowley
and
R.
Hastie
, “
Electron diamagnetism and toroidal coupling of tearing modes
,”
Phys. Fluids
31
(
3
),
426
(
1988
).
13.
S. C.
Cowley
, “
Some aspects of anomalous transport in tokamaks: Stochastic magnetic fields, tearing modes and nonlinear ballooning instabilities
,” Ph.D. thesis (
Princeton University
, NJ,
1985
).
14.
J.
Connor
,
S. C.
Cowley
, and
R.
Hastie
, “
Micro-tearing stability in tokamaks
,”
Plasma Phys. Controlled Fusion
32
(
10
),
799
(
1990
).
15.
P. H.
Rutherford
, “
Nonlinear growth of the tearing mode
,”
Phys. Fluids
16
(
11
),
1903
(
1973
).
16.
F. L.
Waelbroeck
, “
Nonlinear growth of strongly unstable tearing modes
,”
J. Plasma Phy.
50
,
477
494
(
1993
).
17.
A.
Ishizawa
and
F. L.
Waelbroeck
, “
Magnetic island evolution in the presence of ion-temperature gradient-driven turbulence
,”
Phys. Plasmas
20
(
12
),
122301
(
2013
).
18.
W.
Hornsby
,
A.
Peeters
,
M.
Siccinio
, and
E.
Poli
, “
On the dynamics of vortex modes within magnetic islands
,”
Phys. Plasmas
19
,
032308
(
2012
).
19.
F.
Militello
,
F. L.
Waelbroeck
,
R.
Fitzpatrick
, and
W.
Horton
, “
Interaction between turbulence and a nonlinear tearing mode in the low β regime
,”
Phys. Plasmas
15
(
5
),
050701
(
2008
).