Treatment of Advanced Divertor Configurations in the Flux-Coordinate Independent turbulence code GRILLIX

Advanced divertor configurations modify the magnetic geometry of the diverter to achieve a combination of strong magnetic flux expansion, increased connection length and higher divertor volume - to improve detachment stability, neutral/impurity confinement and heat-channel broadening. In this paper, we discuss the modification of the Flux-Coordinate Independent (FCI) turbulence code GRILLIX to treat generalised magnetic geometry, to allow for the investigation of the effect of magnetic geometry on turbulent structures in the edge and SOL. The development of grids and parallel operators from numerically-defined magnetic equilibria is discussed, as is the application of boundary conditions via penalisation, with the finite-width method generalised to treat complex non-conformal boundaries. Initial testing of hyperbolic (advection) and parabolic (diffusion) test cases is presented for the Snowflake scenario.


INTRODUCTION
Limiting the heat-and particle-exhaust to prevent damage to the divertor & first wall is identified as one of the principle challenges in the development of fusion power reactors [4] . For ITER high-power operations, the 'ITER baseline divertor solution' (a single-null 'conventional divertor' (CD) operated at partially detached conditions [4,5] ) will require operation close to the material limits for the divertor targets [6] . Operating conditions for DEMO and future fusion power reactors will be significantly outside the parameter space tested in existing experimental devices and it is uncertain whether the ITER baseline solution will be sufficient for these devices. As a contingency, several 'advanced divertor configurations' (ADCs) are currently the focus of WP-DTT1-ADC project -including the 'Snowflake divertor' (SF, [7,8] ), the 'X-divertor' (XD, [8] ), the 'Super-X divertor' (SXD, [9] ) and a 'Double-null' (DN) variant of a conventional divertor. These configurations are designed to have (compared to the ITER baseline scenario) [4,10] higher poloidal flux expansion -which increases the width of the directed heat-flux channel, to increase the area over which the heat is deposited in the divertor targets -, increased divertor parallel connection lengthwhich gives more field-line-parallel distance over which heat flux can be radiatively dissipated -, and improved neutral baffling -which allows for higher neutral and seeded-impurity densities near the divertor targets, increasing SOL/divertor radiation while keeping core radiation at acceptable levels.

FIGURE 1
Preprocessed Snowflake magnetic equilibria. Contours show the given poloidal flux Ψ (in Wb), and the red polygon indicates the first-wall plus divertor. Additionally, the blue & orange lines gives the upper edge of polygons which provide region-localised flux limiting to eliminate the 'dome' parallel boundary, the blue-shaded region gives the end-ofdomain boundary, and the blue and red crosses give the magnetic axis and primary x-point respectively. In this figure values of Ψ which are outside of the identified flux-limits are excluded, so the contour region is approximately equal to the resulting numerical grid.

TURBULENCE MODELLING WITH GRILLIX
GRILLIX is a drift-reduced Braginskii 3D fluid turbulence code, which is particularly notable for its use of the flux-coordinate independent approach (FCI). For a complete description of the code, see Stegmeir et al., 2019 [1] regarding the use of FCI (including the support operator method), extension to a global 1 electromagnetic model, and the penalisation method, and Zholobenko et al., 2019 [11] , for extensions including ion-thermal effects and implicit heat conductivity. This paper focuses on modifications to the grid generation method and the establishment of parallel operators for numerical equilibria, and a modification which allows for the use of a minimal width penalisation function. Although these developments are made in the context of the GRILLIX code, they are generally applicable to any turbulence code employing FCI.
The complex magnetic geometry of advanced divertor configurations can introduce challenges for codes employing field-or flux-aligned coordinates. Firstly, each of the (or, high order) X-points require extraordinary treatment, to avoid the introduction of coordinate singularities. Secondly, since turbulence occurs on roughly isotropic scales, resolving turbulent dynamics in regions of high flux-expansion (i.e. near poloidal-field nulls, or near flux-expanded strike points) requires very high resolution around the midplane. This can introduce a severe CFL criterion, which requires a greatly reduced timestep (particularly for explicit numerical schemes).
In contrast, the flux-coordinate independent method uses a grid which is independent of the magnetic field and instead encodes the magnetic field structure in the parallel operators. This approach has both advantages and disadvantages compared to the more typical field-or flux-aligned coordinates. Coordinate singularities are avoided, and X-points, O-points and the separatrix can be treated without changing the numerical method. The grid resolution may be set freely, allowing for good resolution in regions of high-flux expansion and with a minimally-restrictive CFL criterion. Additionally, the method allows for the treatment of 3D equilibria such as RMPs and stellarators [12] -although currently GRILLIX assumes an axisymmetric field.
Conversely, the interpolation routines incur a (slight) computational cost outright, but more significantly this also means that efficient parallelisation is more complex since the concept of up and downwind parallel 'neighbours' is not well defined. The field-line tracing and interpolation unavoidable adds a certain degree of numerical dissipation -and since the parallel dynamics are much faster than perpendicular dynamics (flute mode behaviour), this could overwhelm the actual perpendicular dynamics unless special care is taken in the construction of parallel operators (i.e. support operator method). Finally, because the grid is not conformal to the boundary, the application of boundary conditions is complicated significantly.

Preprocessing
The advanced divertor configurations were supplied as eqdsk-standard files, which give the poloidal flux function Ψ as a function of the radial and vertical directions. Additionally, the first-wall and divertor were given as a set of ( , ) polygon-points.
Before importing this data into GRILLIX, a preprocessing step was used to identify features to make the equilibria compatible with the set of equations and boundary conditions in GRILLIX -for the Snowflake scenario these are indicated in figure 1 , including a) a lower limit of Ψ, below which grid points are not included in the computational grid -excluding the core where the fluid model is not valid; b) an upper limit of Ψ which eliminates points which are on a flux surface which impacts the first-wall -giving parallel boundary conditions only at the divertor targets; c) region-specific Ψ limits & corresponding region polygons, which eliminate points which are on a flux-surface which impacts the 'dome' region between two divertor targets; d) a polygon identifying the end of the domain -where simple boundary conditions are applied, to provide an end-of-domain for the penalisation stencil; and e) positions of a primary X-point at ( , ) and the magnetic axis at ( , ) -used for normalisation.

Poloidal flux Ψ and first-wall/divertor polygon  div
The equilibrium data and results from preprocessing are passed to GRILLIX and used to build a bicubic spline interpolator [13] . This allowed for the poloidal flux-function to be determined at any point within the computational grid, as well as the derivatives of Ψ at those points. Due to the axisymmetry of the magnetic field, the field components are defined entirely by the poloidal flux function via ) which automatically fulfils ∇ ⋅ = 0 [14] .
The first-wall and divertor polygon, the flux-limiting regions and the end-of-domain polygon are treated via a custom polygon object, which stores the ( , ) array of polygon points and uses the 'winding algorithm' [15] to identify whether a query point lies within the interior or exterior of the polygon. Together, the Ψ interpolator and the polygons are used to determine the computation grid and ghost points from a full Cartesian grid.

Field-line tracing
Once the grid has been established, the parallel operators are constructed via field-line tracing. For a given start point ( , , ) (where is the toroidal angle coordinate), the coordinate variation along the fieldline is given by = and = and so the fieldline trace from a toroidal angle to = + Δ can be expressed as The integration is performed using the dop853 integration library [16] . The integrator uses an 8th-order Runga-Kutta method with a 5th order error estimator (plus a 3rd order 'dense output' interpolation). The error estimator can be used to automatically adapt the step-size of the integrator to achieve a specified tolerance. By adding a modified 'trace' toroidal angle * to the state vector, the adaptive step-size routine (combined with a stop condition) may be used to perform an accurate conditional trace. A interior trace of * may be constructed by checking whether the state vector coordinates are inside or outside of the first-wall/divertor polygon  div , and 'integrating' * with * = 1 if ( ′ , ′ ∈  div ) 0 otherwise . Once the point crosses out of the firstwall/divertor polygon, the integration is aborted. A corresponding exterior trace can be constructed from the inverse integration condition. This method allows the determination of the trace angle * to the divertor targets in the direction with or against the toroidal field, marked as +, > 0 and −, < 0 respectively. For points outside of the vessel, only one of these angles will The plots are annotated with the separatrix (white) for and divertor (red) for Ξ.

FIGURE 3
Parallel diffusion case, for a blob initiated at ∕ 0 = 0.89, ∕ 0 = −0.495 on a single plane with an initial density of 1. The left image shows the density distribution on a single poloidal plane after = 27.5, as well as the magnetic field-line corresponding to the centre of the blob (black-dashed line). The initial circular blob can be seen in its original position, as well as strongly sheared blobs corresponding to (−3 ≤ ≤ 2) × 2 fullrotation projections of the blob. The ( < −3) and ( > 3) projections are in the penalisation area, and so are not visible. The right plot shows the density profile of the blob in the parallel direction (with ( = 0.5) set equal to 1), as well as the penalisation function (black dashed line). It is seen that the initial density evolution is a symmetric decay. Once the forward density diffusion front reaches the boundary condition (where = 1), the Neumann boundary condition modifies the dynamics. be defined. The external trace angles are defined as ±, = ±, (for Σ > 0 the interior trace angle between the divertor targets) which gives continuous functions + and − . At the targets, ± = ± Σ and ∓ = 0 if the field line is directed towards the target, and vice versa.

Penalisation functions
is the characteristic function, ≪ 1 is the relaxation parameter, and is a penalisation value which is set such that the boundary conditions are satisfied exactly as → 0. Within GRILLIX, the toroidally-projected width of the function must be greater than or equal to the toroidal angle between poloidal planes. This prevents an issue where the staggered and full-grids used for the finite-volume discretisation of parallel dynamics decouple at the boundary. By using a finite width and gradually transitioning from interior to exterior (penalised) dynamics, this decoupling can be avoided. The smoothed penalisation characteristic function can be determined directly from the trace angles + and − via Here  3 is a 3rd order smooth-step function -a sigmoid function which gives a smooth (Hermitepolynomial) transition from 0 to 1, over a region of width . Additionally, a function Ξ = 1 − 2 ⋅ ( + + − ) (where  is a Hermite step-function) is developed to indicate, for a given point, whether the nearest plate is in the direction of (Ξ = +1), or against (Ξ = −1), the toroidal field direction. This is required to indicate the sign for the penalisation values of the parallel velocity and parallel current, as well as to determine which direction is 'upstream' (towards the interior) for boundary conditions which set ghost values in terms of upstream neighbours (i.e. ≥ 1 st order Dirichlet, Neumann and higher derivatives). These functions are shown for the Snowflake divertor in figure 2 .
The use of a smooth function is required for numerical stability in GRILLIX. However, this has the disadvantage that the boundary condition is not applied exactly at the plate, but rather the system of equations smoothly transitions to the boundary condition. This complicates the determination of boundary fluxes, and -in cases where the forcing terms of each variable are highly disparate -may result in the 'boundary' being located at different points for different quantities. As such, choosing a penalisation width as small as possible while still avoiding poloidal plane decoupling is desirable. Future work aims to consider whether other parallel stencils such as WENO could alleviate this issue, or alternatively whether an alternative set of nondecoupling boundary conditions [17] could be made compatible with the implemented turbulence model. Of particular interest is the set of 'strong-sink boundary conditions' developed by Parades et al. [18] , which could allow for a unified set of boundary conditions for both the divertor and first wall.

TESTING OF NUMERICAL EQUILIBRIA
The physical model in GRILLIX is expressed as a coupled set of non-linear partial differential equations (see [1] ), which may be broadly defined as either parabolic, hyperbolic or elliptic PDEs, or a mix thereof. We consider prototypical parabolic and hyperbolic equations -the parallel diffusion and parallel advection equations -to test the implementation of the numerical equilibria in GRILLIX.
The parallel advection equation is given by + ∇ ⋅ ( ∥̂ ) = 0 and ∥ + ∥ ∇ ∥ ∥ = − ∇ ∥ . The density is located on the full-grid , while the parallel velocity is located on the staggered grid * . These grids are coupled via the parallel operators ∇ ∥ (...) ∶ → * which maps from full to staggered grid, and ∇ ⋅ (̂ ...) ∶ * → which maps from staggered to full. The boundary conditions applied are Bohm for the parallel velocity ( ∥ → ± = ±1), and Neumann for the density (∇ = 0). The time normalisation is equal to 0 for this system of equations. The results of the advection test are shown in figure 4 , and discussed in the corresponding caption.
The parallel diffusion equation is = ∇ ⋅ (̂ ∇ ∥ ), with Neumann density boundary conditions. The time normalisation for this equation is 2 0 . The results from the diffusion test are shown in figure 3 , and discussed in the corresponding caption.

CONCLUSIONS AND FUTURE WORK
Advanced divertor configurations aim to reduce the heat flux impacting the divertor targets via increased flux expansion & parallel connection lengths, and improved control of detachment. To assist with the interpretation of experimental results, modelling with turbulence codes is desirable. The flux-coordinate independent approach is expected to be well-suited to studying complex magnetic geometries due to the ability to set the grid independent of the magnetic field structure. However, there are also challenges associated with this method -particularly the application of boundary conditions for boundaries which are non-conformal to the grid. In this paper, we discuss an extension of the penalisation method discussed in [1] , which allows for minimisation of the penalisation-width required for coupling the staggered and full-grids at the boundary. Additionally, an automatic method for determining the direction to the nearest plate based on a conditional field-line tracing is introduced. The modifications to the code are applied to prototypical parabolic (parallel-advection) and hyperbolic (paralleldiffusion) test cases, which represent the two classes of differential equation which are affected by the magnetic field structure through the parallel operators. The parallel diffusion case is seen to decay as expected initially and is modified by the Neumann boundary condition for late-stage evolution. Blob projections are visible at (2 ) projections along the field-line. The separatrix and X-points formation is clearly seen in the advection test case.
Work to apply the full turbulence model on an ADC case is ongoing, with no saturated results available yet.