Fluid simulations of plasma filaments in stellarator geometries with BSTING

Here we present first results simulating plasma filaments in non-axisymmetric geometries, using a fluid turbulence extension of the \boutxx~framework. This is made possible by the implementation of the Flux Coordinate Independent scheme for parallel derivatives, an extension of the metric tensor components which allows them to vary in three dimensions, and development of grid generation. Tests have been performed to confirm that the extension to three dimensional metric tensors does not compromise the accuracy and stability of the associated numerical operators. Recent changes to the FCI grid generator in \boutxx, including a curvilinear grid system which allows for potentially more efficient computation, are also presented. Initial simulations of seeded plasma filaments in a non-axisymmetric geometry are reported. We characterize filaments propagating in the closed-field-line region of a low-field-period, rotating ellipse equilibrium as inertially-limited by examining the velocity scaling and currents associated with the filament propagation. Finally, it is shown that filaments in a non-axisymmetric rotating ellipse equilibrium propagate in a toroidally nonuniform fashion, and it is determined that the long connection lengths in the scrape-off-layer enable parallel gradients to establish, which has consequences for interpretation of experimental data.


I. INTRODUCTION
Neoclassical transport is the dominant loss mechanism in sufficiently hot stellarator plasmas and can dominate in the plasma core [1]. In the outer, colder parts of the plasma, however, turbulence becomes more important and therefore dominates the plasma edge region [2]. Since the Wendelstein 7-X stellarator [3] has been optimized to have low neoclassical transport, turbulent transport could become comparable to neoclassical losses even in the center of the plasma. Wendelstein 7-X has already demonstrated novel edge physics; poloidally rotating filaments as measured by visible cameras [4], and a high-frequency variation of limiter heat fluxes [5] merit numerical investigation. Furthermore, the edge of Wendelstein 7-X in the island divertor configuration exhibits long connection lengths, such that cross field transport can become comparable to parallel transport. Predicting this crossfield transport in high density, collisional, detached plasmas without an ad-hoc assumption for diffusion is a motivation of this work. It is becoming increasingly important to simulate turbulence in non-axisymmetric configurations.
In stellarator core plasmas, the most common method for simulating plasma turbulence is with gyrokinetic codes such as GENE [6], which is feasible due to the closed flux surfaces and the low collisionality. However, the simulations are computationally expensive for long (on the order of confinement time) temporal and global spatial scales. Additionally, GENE simulations are currently limited to flux-tube and flux-tube-ensemble geometries.
The high collisionality of tokamak and stellarator edge plasmas facilitates a fluid approach to turbulence simulations. While there are several fluid turbulence simulation codes for tokamak geometries [7][8][9], previous attempts to develop such a simulation framework for stellarators have been unsuccessful.
The recent implementation of the Flux Coordinate Independent (FCI) [10] method for parallel derivatives in BOUT++ has allowed for simulations in non-axisymmetric geometries [11,12]. Instead of aligning the computational grid to magnetic field lines, the FCI method uses interpolation of field line mapping on poloidal (or, in the case of linear geometries, azimuthal) planes to obtain values for finite-difference differentiation parallel to the magnetic field. In BOUT++, a cubic Hermite spline is utilized, although other methods have been implemented [12]. The FCI method removes the inherent singularities in flux or field aligned coordinates around magnetic null points. Additionally, since the computa-tional grid is no longer aligned to the magnetic field, the simulation of complex geometries including X-points is possible. For a more complete discussion of the FCI method, see References [10][11][12].
Here, we present the first results simulating plasma fluid turbulence in non-axisymmetric geometries, made possible by extensive modifications to the BOUT++ framework [13,14].
Section I A describes the recent modifications to the BOUT++ framework which are relevant for this work. Initial testing of the modified framework is described in Section II, where Sections II A and II B test the accuracy FCI parallel gradient operators and their associated boundary conditions, and Section II C reports the modifications to the Laplacian inversion algorithms. Section III introduces a new curvilinear coordinate system for FCI simulations in BOUT++ which is used in Section IV to simulate plasma filaments in non-axisymmetric geometries; filaments in the closed-field-line region of a rotating ellipse geometry are determined to be inertially-limited and exhibit a toroidally non-uniform propagation, a result which has implications for interpretation of experimental data. Finally, Section V describes how the curvilinear FCI grids can be used for simulation of realistic geometries, namely Wendelstein 7-X.

A. Modifications to the BOUT++ framework
The BOUT++ framework is a modular, object oriented and open source framework for fluid simulations with an international team of developers [13]. This paper presents recent progress in modifying BOUT++ to Simulate Turbulence In Non-axisymmetric Geometries under the "BSTING" project.
Previous work in simulating non-axisymmetric geometries has focused on the conventional BOUT++ framework, which is a 3D code but was written with metric tensor components which vary in two dimensions due to an assumption of toroidal symmetry. For an accurate simulation of plasma dynamics in stellarators, BSTING must include metric components which are fully three dimensional. This extension to three dimensions is simple in principle (and was in fact mentioned in the introduction of the original BOUT++ paper [13]), but unfortunately the geometrical components are integral to many different parts of the code, and the work presented here has required extensive modifications to the framework.
The majority of modifications are primarily focused on the numerical methods of spatial operators and do not affect file handling, parallelization, post processing, and many other functions in BOUT++. Development has focused on implementing operators relevant to edge transport and turbulence simulations: spatial derivatives of scalar fields which vary in three dimensions, and Laplacian inversion. Here we address the most relevant issues: the accuracy of spatial gradient operators, boundary condition implementation, and Laplacian inversion which allows plasma potential to be calculated from vorticity. The following section provides initial tests for the implementation of these methods.

II. TESTING
The development of BSTING is an extensive modification to the BOUT++ framework, and therefore careful testing of numerical accuracy is required. In this section, we concentrate on ensuring the accuracy of spatial derivatives, boundary conditions, and Laplacian inversion.
All tests in this section use a geometry where the poloidal planes are described by the radial x-coordinate and vertical z-coordinate while the y-coordinate describes the toroidal (or longitudinal in linear geometries) direction. The FCI operators therefore interpolate the relevant values based on field line mapping in the x-z planes. In Section III we will discuss an alternative coordinate system for complex geometries.
A. Flux surface mapping using heat diffusion A potential issue with the implementation of the FCI scheme as discussed in section I A is that since the poloidal planes are not orthogonal to the magnetic field lines, there could be a considerable pollution of perpendicular dynamics due to the projection of parallel effects [15].
A simple and common test to ensure the proper calculation of parallel dynamics using the FCI method in complex geometries is to implement a parallel diffusion model such as that shown in Equation 1.
where b is the magnetic field vector. Here the diffusion model in Equation 1 is used to test the numerical diffusion in a rotating ellipse equilibrium as done in References [11,12].
Specifically, we will simulate this model on a rotating ellipse geometry, the flux surfaces of  [11,12] -however this result differs in that it uses fully three dimensional metric tensor components, whereas the previous results utilized a metric tensor that varied in only two dimensions. This added flexibility also allows for non-axisymmetric toroidal geometries. Figure 3 indicates the flux surfaces as calculated by BSTING in a toroidal rotating ellipse geometry. The red surfaces indicate the 2D projection on each poloidal plane, and the blue/green cloud is the interpolated function between the poloidal planes.
This heat flux mapping indicates that the FCI operators are capable of simulating nonaxisymmetric geometries after the transition to three dimensional metric tensors in BSTING.
The following section will use a more quantitative method to ensure the numerical operators and implementation of boundary conditions with three dimensional metric tensors have sufficiently small numerical error. lines can leave the domain before reaching the next toroidal plane, therefore leading to non-uniform grid point spacing for interpolation and complicating the correct calculation of derivatives. There have been a few recent advances in boundary condition calculation for FCI operators; BOUT++ utilizes the Leg-Value-Fill (LVF) method detailed in Reference [12].
In this section we extend previous testing [12] using the Method of Manufactured Solutions [18,19] to ensure that the extension to three dimensional metric tensors has not diminished the accuracy and stability of the framework. Two coupled differential equations were therefore simulated for a single time step: where parameters are identical to those in Reference [12]; namely, D=10, and the domain measures 0.1 x 10 x 1 (x,y,z) meters. The magnetic geometry is a sheared slab, such that (B x ,B y ,B z ) = (0,1,0.05 + (x-0.05)/10). The manufactured solutions are also those from Reference [12]: whereȳ andz are normalized between 0 and 2π. The diffusion terms in Equations 2 and 3 scale with y-spacing, and do not affect the convergence of ∇ . Therefore the grid is scaled in y and z simultaneously. Figure 4 indicates the convergence of FCI operators in BSTING, including LVF boundary conditions.  Having established the accuracy and stability of the FCI operators and the associated LVF boundary conditions in BSTING, the following section describes the implementation of Laplacian inversion routines which allow for the calculation of plasma potential from vorticity.

C. Laplacian inversion with complete poloidal metrics
One of the advantages of BOUT++ is its modular nature; numerical methods can be modified without compromising the stability or accuracy of the rest of the framework. For this reason, several different methods for Laplacian inversion have been implemented in BOUT++.
Unfortunately for BSTING, many of these routines assume a periodicity in one direction (the z coordinate, usually the toroidal angle in tokamak simulations), since BOUT++ was originally designed to simulate turbulence in tokamak scrape-off-layers. Recent work on implementing the Hermes model [20] in BOUT++ has included several new numerical methods. One of these is the implementation of a Laplacian inversion routine in three dimensions, which inverts an inhomogeneous Helmholtz equation in the conservative form: where where J is the Jacobian, g ij are the metric tensor components, and A, B and b are variables which are specific to each situation -for instance b is often vorticity in plasma turbulence simulations. The current implementation of this solver utilizes the PETSc suite of data routines [21], which is available with several features including preconditioners for efficient computation. This implementation differs from conventional BOUT++ since it includes the off-diagonal metric terms (g xz ). By setting the metric tensor components, g ij , to non-zero values and comparing the implemented inversion routine using PETSc to explicit calculation of Equation 7 indicated a difference of less than 10 −15 . Testing with zero-value diagonal metric tensor components indicated similar errors relative to the implementation without off-diagonal metrics in BOUT++, suggesting proper convergence of the inversion routines.
Having implemented the FCI operators and Laplacian inversion with Cartesian poloidal grids, the BSTING project is now capable of simulating turbulence in non-axisymmetric geometries. A significant challenge for this method, however, is to handle the entire plasma cross section in a Cartesian poloidal grid while neglecting the plasma core and far edge.
One solution to this issue is to use a penalization function to mask the areas where the variables should not be evolved. This method has been used previously in BOUT++ [16] to remove solid-density magnetic coils in the simulation domain and is currently used in with FCI operators in GRILLIX [9] to mask the plasma core and far scrape-off-layer. The disadvantage of this method is that it requires a large poloidal grid for a relatively small computational area. In the following section we present a new method for generating FCI grids in BOUT++ and BSTING which does not use a grid over the entire plasma cross section, potentially providing faster computation.

A. Implementation of Elliptic Grids
While all previous simulations using the FCI method have used poloidal planes with Cartesian coordinates [9][10][11][12]22], this is not required. The method is independent of the poloidal grid system as long as interpolation in these planes is correctly calculated and communicated. Here we present recent results using structured, non-Cartesian poloidal grids which are still logically rectangular [23,24]. As an illustration of this method, Figure 5 illustrates a sample grid with independent inner and outer surfaces.
These new grids have been added to the BOUT++ FCI grid generator, Zoidberg, and are included in a recent release of BOUT++ (version 4.1). These grids are particularly advantageous as they include a periodic direction which could potentially increase computational efficiency. A grid is generated by prescribing an inner and outer surface, and then inverting an elliptic equation to connect the inner and outer points. Both the inner and outer surface shapes are independently prescribed, and can be described using various methods: Zoidberg includes an flux surface shape generator, which will describe a shape based on elongation, triangularity and indentation. Alternatively, one can use the Zoidberg field line tracer to construct flux surfaces from a given magnetic field (i.e. from VMEC, a vacuum field solver, or an analytic magnetic field description), and generate a shape based on this flux surface mapping.
These grids provide an additional degree of flexibility and avoid some potential prob-  This two-field period, rotating ellipse geometry has a major radius of 2.5m, The inner surface is described by a flux surface, but the rest of the grid is not aligned to flux surfaces; the outer surface is a circle centered around the magnetic axis with a radius of 50cm.
Therefore, this geometry incorporates both open and closed field lines. The following section utilizes a finite-β electromagnetic isothermal reduced magnetohydrodynamic model similar to that used in the isothermal version of TOKAM3X [26] which evolves vorticity ω, electromagnetic potential A , electron density n, and parallel momen-tum Γ = m i nv . Electron and ion temperatures T e and T i are assumed constant, though independently specified. The magnetic field is described by a constant equilibrium field B 0 and a time-evolving poloidal field such that: where A is the parallel component of the vector potential and a large-aspect ratio approximation has been utilized such that ψ = RA .
The equations are described as follows in SI units: Here ∂ ≡ b·∇ and ∇ f ≡ ∇·(bf ) = B∂ f B . The pressure is p = p e +p i = n(T e +T i ). The vector b 0 = e φ is the "toroidal" magnetic field unit vector, and b = B/B 0 is the unit vector along the total magnetic field, assuming that the poloidal magnetic field is small relative to the toroidal field. Gradients in the poloidal plane, which is not necessarily perpendicular to the magnetic field (in the case using FCI derivatives, as is used here), are defined by In this model, the magnetic drift term is treated generally (in comparison to, for instance, Equation A.14) and is written as: which uses ∇ × B · ∇p = 0 which is valid in equilibrium since J · ∇p = 0. The curvature operator is then defined as: which has a similar form as that derived in the appendix (Equation A.8), meaning that we can use the bracket coefficient to calculate the curvature effects in curvilinear grids. This is especially convenient as the magnetic field does not, in general, vary solely with the major radius in stellarators -an approximation which is often used in fluid turbulence simulations [16,27,28]. In the simulations presented here, all cross-field drifts are implemented with the 2 nd order Arakawa brackets [29].
B. A weakly-non-axisymmetric, rotating-ellipse geometry As an initial investigation of turbulence in non-axisymmetric geometries, a seeded plasma filament in a rotating ellipse geometry was considered. While there have been experimental investigations of turbulent filaments in stellarators [30], this study will serve as the first example of fluid turbulence simulations in non-axisymmetric geometries. A seeded filament test offers a somewhat straightforward approach to studying important phenomena in plasma transport. Previous studies in BOUT++ have investigated filaments in slab [31], toroidal pinch [27], and X-point geometries [16,32].
For the studies presented here, an analytically calculated, low-field-period rotating ellipse geometry was chosen due to the relatively straightforward implementation and analysis.
These analytic equilibria are a necessary step before geometries like W7-X. Wendelstein 7-X grids for use in BSTING are described in Section V, but turbulence studies in these more complex geometries will be a subject of further study. Furthermore, low-field-period rotating ellipse geometries exhibit a magnetic field which generally varies as 1/R (see Figure   1 from [33]), allowing for a more straightforward analysis since this configuration is most similar to axisymmetric configurations. Figure 7 illustrates the degree of non-axisymmetry by plotting the variation of the magnetic field multiplied by the major radius, since a plot of the magnetic field strength would be dominated by the predominantly 1/R variation.

FIG. 7: Variation of the non-toroidal magnetic field at three different toroidal locations -
obtained by multiplying the total field by the major radius R, and calculating the difference with respect to the mean value.
From Figure 7 it can be deduced that the magnitude magnetic field which does not vary like 1/R only changes toroidally by less than a percent, indicating a small degree of nonaxisymmetry in the magnetic field strength, which can affect the drive term for filament propagation (Equation 19).

C. Filament characterization
To characterize filament propagation in this non-axisymmetric geometry, a field-aligned plasma filament is first initialized; an approximately circular density perturbation at (R,Z,φ) = (2.5m,-0.3m, 0.0) is prescribed and a simple parallel diffusion model as in Equation 1 is first simulated to achieve an initial condition of a field-aligned filament. As this is a lowshear geometry, the filament approximately becomes field aligned once the initial distribution diffuses once toroidally. The initial field-alignment is determined when the maximum value of the density on a plane varies by less than 5% in a timestep (100/ω Ci , where ω Ci is the ion cyclotron frequency). This condition is satisfied after 100 timesteps, or ten thousand ion cyclotron times. This field-aligned density distribution, where the peak density perturbation is n = 1.05 × 10 19 m −3 , is then used as an initial condition for the seeded filament simulation using the model described in Section IV A. All other plasma fields are not initialized and, once the field-aligned filament is achieved, are allowed to develop independently. The ion and electron temperature is set to 100eV and the background density is n 0 = 1 × 10 19 m −3 .
Plasma filaments (or blobs) are often characterized by the method by which the charge separation is resolved; if charge is carried via parallel currents through the sheath, filaments are considered "sheath limited". If the connection length to the sheath is large, however, this charge separation can be short-circuited via perpendicular currents and the filaments propagate in a so-called "inertially-limited" regime [34]. Filament propagation is also characterized by the scaling of the propagation speed as a function of its poloidal cross section, δ ⊥ ; inertially limited filaments scale proportional to δ 1/2 ⊥ , whereas sheath-limited filaments scale as δ −2 ⊥ . For more complete discussions of filaments, see References [34,35]. Therefore, one can determine the filament propagation regime by plotting the scaling of the maximum speed as a function of filament diameter δ ⊥ . The edge and scrape-offlayer of stellarators such as Wendelstein 7-X can exhibit large connection lengths [5]. As an initial insight into filament behavior in a non-axisymmetric field with long connection lengths, filaments were seeded in the closed-field-line region in the weakly non-axisymmetric geometry discussed in the previous sections. The scaling of these filaments is shown in Figure 8, where δ ⊥ = 1 is normalized to 7cm, the initial filament diameter for the filaments in the following section (IV D).
Similar to the tokamak (axisymmetric) case, the scaling of filaments initialized in the closed-field-line region propagate in an inertially-limited regime, as indicated by the δ 1/2 ⊥ scaling in Figure 8. As a confirmation of the inertially-limited propagation, Figure 9 illustrates the currents which dictate the propagation of the filament at t ≈ 4µs.
Since the divergence of the parallel current is much smaller than the perpendicular currents, the potential difference is resolved via short-circuiting perpendicular currents, instead of traveling along field lines to the sheath. This again supports the characterization an inertially-limited regime. As this is only a weakly-non-axisymmetric field, it is reasonable to find similarities to filaments in an axisymmetric field, for instance in Reference [32], where inertially-limited filaments were characterized in a MAST (tokamak) geometry. For a more strongly-non-axisymmetric geometry such as Wendelstein 7-X, the filament propagation may exhibit different behavior, since the filament drive changes directions relative to the major radius within a field period. While filament simulations in Wendelstein 7-X await a future publication, the following section discusses how even a weakly-non-axisymmetric field can alter the toroidal uniformity of the filament propagation.
If the magnetic geometry is not axisymmetric, the filament drive due to the magnetic field curvature can vary along the length of a filament. If the drive is toroidally nonuniform, one would expect the propagation to also vary toroidally. It is often assumed, however, that filaments propagate uniformly along field lines, for instance in [30]. To test the effects of a non-axisymmetric magnetic field, we can investigate the propagation of a filament at different toroidal locations. Figure 10 illustrates the filament velocity (solid) and displacement (dotted) of a 100eV plasma filament at various toroidal angles.  Figure 11 which indicates how the magnetic drive term (black, also fitted), and the resulting filament velocity vary as a function of toroidal angle. Here, the filament velocity is normalized to the average toroidal field at 100 timesteps.
The blue squares in Figure 11 indicate the normalized velocity at each toroidal position, averaged over the 100 timesteps. The fill cloud indicates the standard deviation of the toroidally-normalized velocity for these sample timesteps.
The non-axisymmetric propagation of filaments can be clarified by considering the timescales associated with filament propagation. First, we approximate the timescale for parallel propagation along a filament to follow the relation t ∼ l cs where l is the length along the filament and c s is the ion sound speed. In the simulations presented here, c s ≈ 6.9 × 10 4 m/s, which indicates that information takes about 14µs to propagate one meter. Therefore, if the filament is driven non-uniformly, the time which the filament needs to restore the symmetry is longer than the propagation timescale t ⊥ , which can be approximated by assuming L ⊥ ≈ δ ⊥ ≈ 7cm and v ⊥ ≈ 13km/s -indicating therefore that This assertion can be tested by increasing the speed at which this restoration is performed, for instance by increasing the sound speed. When simulations were performed with hotter (1keV), smaller filaments -thus keeping the pressure constant -the standard deviation of the position of the filaments averaged 79% of that for the colder simulation, indicating that a hotter filament propagates more uniformly. This can also be seen in the resulting speed of the hotter filament, shown as red triangles in Figure 12, which does not vary as strongly with toroidal location. It is also possible, however, that the filament is restored to uniform propagation toroidally at the Alfvén velocity. This would also explain the more uniform propagation for a hotter filament, since the density perturbation was reduced to provide an equal drive (from pressure), and the Alfvén velocity is a function of the plasma β. To determine the extent to which this non-axisymmetric nature is affected by the Alfvénic effects, one can simulate a filament in an electrostatic case. In an electrostatic case, all terms in the model described in Section IV A which are dependent on the plasma β are neglected, which in essence provides an infinite Alfvén speed. Figure 13 illustrates how the propagation of a filament in an electrostatic and electromagnetic filament compare as a function of toroidal angle. Figure 13 indicates that the non-uniform propagation is not an electromagnetic effect and thus cannot be adequately mitigated by parallel transport at infinite Alfvénic speeds, since the electrostatic and electromagnetic case exhibit very similar characteristics.

V. WENDELSTEIN 7-X CURVILINEAR GRIDS
BSTING is designed to provide numerical support for experimental measurements. The curvilinear grid system presented in Section III has therefore been applied to Wendelstein 7-X geometries using various descriptions of the magnetic field. As this geometry is considerably more complicated than the analytically-prescribed rotating-ellipse equilibria presented earlier, the following sections extend the flux surface mapping tests to the W7-X grids.

A. Inherent Perpendicular Diffusion in W7-X Curvilinear Grids
Here we present the development of curvilinear poloidal grids for Wendelstein 7-X geometries using outputs from the VMEC code [36]. To test the implementation and limitations of grids in this complicated geometry, the parallel diffusion model in Section II A, Equation 1 was modified to include a perpendicular diffusion, as shown in Equation 21. . This corresponds to a resolution of approximately 0.3mm -although this obviously is not uniform -which is a relatively coarse resolution for a Wendelstein 7-X turbulence study (ρ s ≈ 0.1mm). Figure 15 indicates that the inherent numerical perpendicular diffusion caused by pollution from parallel dynamics is less than a factor of 10 −9 times smaller than the parallel diffusion, as this is where the points begin to diverge significantly from the zero-diffusion case (as indicated by the dashed line at f 150 /f 0 = 1.0). This inherent perpendicular diffusion is sufficiently less than transport due to plasma drifts and turbulence [15]. This is encouraging as this result is for a moderate-resolution grid, and higher-resolution grids will most likely be necessary for future turbulence simulations in Wendelstein 7-X.
B. W7-X curvilinear poloidal grid for the edge and scrape-off-layer The grids described in the previous section are generated from VMEC [36,37] equilibria, which assume closed flux surfaces. The edge of Wendelstein 7-X is much more complex as it includes magnetic islands and stochastic magnetic field lines. As such, another tool must be developed to trace field lines for grids which can accurately describe this region. To this and 16c illustrate one such grid, which uses the Wendelstein 7-X web services vacuum field solver and components database [38]. The inner surface is generated by tracing flux surfaces using a vacuum field solver, which simplifies core boundary conditions and potential coupling to core profiles and sources, and the outer surface is generated based on a description of the Wendelstein 7-X divertor and first wall developed by Michael Drevlak for fast particle calculations, and is also available on the Wendelstein 7-X webservices. Zoidberg has also been modified to use EXTENDER [39], allowing both plasma-generated magnetic fields and a smooth vacuum solution outside of the last closed flux surface.

VI. CONCLUSIONS
The BOUT++ framework has been extended to allow the metric tensors to vary in three dimensions. This provides greater flexibility to the framework. One major advancement is the implementation of curvilinear grids for use with the Flux Coordinate Independent (FCI) method. Initial simulations of filament propagation in non-axisymmetric geometry have been performed, and the filaments have been characterized to propagate in the inertiallylimited regime. Furthermore, simulations indicate that even a weakly-non-axisymmetric field can significantly alter the propagation of filaments. The long connection lengths of the scrape-off-layer in non-axisymmetric geometries facilitates the establishment of parallel nonuniformity, an effect which must be considered when interpreting experimental data.
Since three dimensional effects are becoming increasingly important -for instance the application of edge magnetic perturbations -the results presented here are applicable to both tokamak and stellarator configurations.
Future work will include simulations of filaments in the Wendelstein 7-X stellarator, where the non-uniform drive of a filament can be more pronounced. The curvature drive in Wendelstein 7-X reverses direction relative to the major radius within a single field period, which could lead to highly non-uniform propagation of filaments, or perhaps even prohibit the radial propagation of coherent filament structures.

VII. ACKNOWLEDGMENTS
The authors would like to acknowledge the work of the BOUT++ development team. The The poloidally curvilinear coordinate system used in this work dictates that numerical operators in the perpendicular (x-z) plane must carefully incorporate the geometry into the calculation. Here we will concentrate on two operators in particular -Poisson brackets and an example of a curvature operator.
The operator 1 B b × ∇g · ∇f appears often in plasma models and represents phenomena such as E×B advection. It often appears in equations in the form of Poisson brackets, and is what is referred to here as the bracket operator. To determine the modifications for the bracket operator in BOUT++, we start by defining real space coordinates R(x, z) and Z(x, z) which depend on the radial coordinate x and the poloidal coordinate z. In the current formulation, x ranges from 0 to 1, and z from 0 to 2π. From here, we determine the coordinate vectors by taking derivatives along the real-space coordinates: where x i is either the x or z coordinate. We can now define the metric components as: g xx = e x · e x g xz = e x · e z g zz = e z · e z (A. 2) The y-direction is considered to be orthogonal to the x − z plane, and is defined as the toroidal angle spanning 0 to 2π. The nonzero metric components are therefore simply: where R is the major radius. The unit vector b is considered to be perpendicular to the x − z plane and is defined as: b = e y √ e y · e y = e y √ g yy = ∇y √ g yy (A.4) We can then begin to construct the bracket operator by taking: The terms in the square brackets is defined as the Poisson bracket, which is what is conventionally described in BOUT++ by the bracket operator. Noting this, we arrive finally at: where we see that a coefficient of √ gyy JB is required for proper calculation of E×B advection in curvilinear grids. In Clebsch coordinates, however, it is worth noting that ∇z × ∇x = 1 J e y = B and therefore √ g yy /J = B and this coefficient becomes 1.
Curvature effects are one of the most important aspects of turbulence simulations, as this can drive drifts and ballooning behavior which contributes to radial transport. The introduction of curvilinear poloidal grids has necessitated careful implementation of curvature operators. To determine the effects of curvature on a quantity f , we must determine how to calculate (b × κ) · ∇f . As a simple example to illustrate this, we begin by assuming that the curvature vector is of the form: which, when dotted with ∇f , then allows the inclusion of curvature effects in curvilinear poloidal grids. This form of the curvature operator can then be used for large-aspect ratio simulations where the magnetic field varies inversely with major radius, an approximation which is often used in plasma fluid turbulence simulations [16,32]. A more general curvature operator is derived in Section IV A.