Novel SOLPS-ITER simulations of X-point target and snowflake divertors

The design and understanding of alternative divertor configurations may be crucial for achieving acceptable steady-state heat and particle material loads for magnetic confinement fusion reactors. Multiple X-point alternative divertor geometries such as snowflakes and X-point targets have great potential in reducing power loads, but have not yet been simulated widely in codes with kinetic neutrals. This paper discusses recent changes made to the SOLPS-ITER code to allow for the simulation of X-point target and low-field side snowflake divertor geometries. Snowflake simulations using this method are presented, in addition to the first SOLPS-ITER simulation of the X-point target. Analysis of these results show reasonable consistency with the simple modelling and theoretical predictions, supporting the validity of the methodology implemented.


Introduction
Maintaining manageable heat and particle loads on materials presents a significant challenge to overcome for high power magnetic confinement fusion reactors. For reactor-like machines such as ITER and SPARC [1], empirical scalings predict deposited target heat loads in excess of 25 M Wm −2 [2,3]; significantly higher than the ≈15 M Wm −2 steady-state * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. surface heat load limit on tungsten monoblock tiles [4]. Additionally, the high loads of energetic particles in these machines can lead to erosion rates of more than 15 nm s −1 [2] on tungsten tiles, limiting the effective lifetime of these tiles and requiring more frequent, costly replacement of the divertor components. As a consequence, methods of reducing heat and particle loads in the divertor will be important for operating future devices; methods such as injecting radiating impurities and operating with alternative divertors.
Alternative divertors are divertors which leverage novel magnetic or geometric features differing from the 'standard' divertor; features which are predicted to benefit divertor exhaust in some way. One such alternative divertor is the Snowflake, which implements additional nulls near primary X-point, creating a higher order X-point, which has been shown to significantly increase divertor connection lengths [5] and may enhance cross-field transport through churning modes [6]. The X-point target (XPT) is another alternative divertor, which has a secondary X-point near the divertor targets [7]. The low poloidal field near this secondary X-point can lead to longer connection lengths, and there can be significant power sharing between the multiple outer divertor targets-depending on the separation of the two separatrixes and strength of cross-field transport [8,9]. Furthermore, the secondary X-point can be located at higher major radius than the primary, creating a magnetic flux expansion and lowering heat flux densities at the targets (similar to the Super-X divertor configuration [10]). Finally, the low poloidal field and local minimum in total field created by the secondary X-point could aid in detachment control [11,12].
Despite the potential of Snowflake and XPT divertors, neither geometry has been widely simulated in codes with kinetic neutrals, particularly in the simulation code SOLPS-ITER [13,14]. In fact, though snowflakes have been simulated in codes such as EMC3-EIRENE [15][16][17], they have only been simulated once before in SOLPS-ITER [18], and XPTs have never been simulated previously using the code due to the rigidity of the physics modules and grid generators [14]. The benefits of using SOLPS-ITER include the complex kinetic treatment of neutrals through its coupling to the EIRENE Monte-Carlo solver, the ability to activate drifts and currents, and the extensive history and benchmarking of the code [14]. In this work we propose a new generalised method for simulating both X-point target and low-field side snowflake divertors within SOLPS-ITER, which may allow more widespread simulation of these alternative geometries; furthering the understanding of such divertors. Building on work performed in [18], the modifications to SOLPS-ITER have been generalized and are integrated with the rest of the code (version 3.0.8) with no loss of functionality.

Methodology
When running a SOLPS-ITER simulation of a new geometry, it is first necessary to generate relevant input files. This is typically done using a combination of the DivGeo software for basic input specification, the standard grid builder CARRE [19], and other preprocessing routines. After generating the relevant input files, the physics simulations can be run, which for SOLPS-ITER is typically a coupled simulation of the fluid solver B2.5 and the kinetic Monte Carlo code EIRENE [20].
However, previous versions of SOLPS-ITER only work for an expected list of geometries, including single null (SN) divertors, double null divertors, slab geometries, annular geometries, limited geometries and stellarator islands. For more unconventional grids such as the XPT, there are two main issues that prevent using this standard workflow. First, the physics modules of SOLPS-ITER expect a double null geometry if a grid is given that contains two X-points. The second issue is that the grid builder CARRE is unable to build grids using equilibria with two X-points not connected to the same O-point.
Consequently, to run a SOLPS-ITER simulation of a Snowflake or XPT requires solutions to these two issues. Solutions to these challenges have been developed throughout this work and will be presented in the following section.

Modifying SOLPS-ITER source
The first challenge to be solved before XPT or Snowflakes can be simulated in SOLPS-ITER is the fact that the fundamental physics modules are not expecting such geometries. Instead, when two X-points are detected in the geometry, SOLPS-ITER assumes the grid is a double null. These assumptions are present because the B2.5 plasma grid must be rectangular when mapped into poloidal and radial index coordinates.The introduction of X-points or more than two targets introduces 'cuts' into the grid, where the poloidal neighbour of a cell is not the cell that is adjacent in physical space. To prepare the geometry correctly, SOLPS-ITER assumes a certain order of Xpoint and target cuts. If two disconnected X-points are detected, the code assumes the ordering is that of a disconnected double null (DDN), which is shown in figure 1. Figure 2 on the other hand shows the geometry of a lower XPT and its rectangular mapping, which is similar to a snowflake minus with the secondary X-point on the low-field side. A low-field side snowflake plus is poloidally similar as well, but with different radial region definitions. From figures 1 and 2, one can see that the DDN and XPT geometries are similar, with three radial regions, four cuts due to X-points, and one cut due to intermediate targets. However, the figure also shows that the ordering of these cuts is subtly but crucially different. In the DDN case the target cut is the 3rd cut to be encountered poloidally, whereas in the XPT case the intermediate targets comprise the 4th poloidal cut. Moreover, if we label the primary X-point as 1, and the secondary as 2, then in order of increasing in poloidal index the DDN encounters 1,2,2,1; contrarily, the XPT encounters 1,1,2,2.
Because the ordering of XPT cuts are different from what SOLPS-ITER expects, a standard run of the code will fail if an X-point target or Snowflake grid is provided. To overcome this, several B2.5 preprocessing routines, and B2.5-EIRENE interfacing routines were edited to detect whether a grid has a DDN or an XPT (equivalently a low-field side snowflake) ordering of primary and secondary X-point cuts. If a DDN is detected, the code will run in its unedited form, expecting the cut order shown in figure 1. If an XPT or snowflake is detected, however, new code has been written to unpack grid data according to an ordering of cuts shown in figure 2. With this modified code, XPTs and low-field side snowflakes (both plus and minus) can be simulated without taking away from any previous SOLPS-ITER functionality. The code modifications have not yet been implemented to allow for a high-field side snowflake, though this could be achieved in much the same way.

Input file generation
The second issue to overcome in simulating XPTs or snowflakes in SOLPS-ITER is the fact that relevant input files  cannot be generated through the standard grid builder CARRE. As such, a custom workflow was developed to prepare the required input files for a SOLPS-ITER run without the use of CARRE. This workflow consists of three unique steps, followed by a continuation of the regular SOLPS-ITER workflow.
The first step in the custom XPT simulation workflow is similar to that of a regular workflow, in that it revolves around generating a basic input structure with the graphical user interface DivGeo. Plasma facing components (PFCs) are loaded into DivGeo and adjusted manually if necessary. The plasma species, wall material, and recycling coefficients are specified. The second step of this process is generating the XPT grid mesh from an equilibrium. To do this, we have used INGRID [21], an open-source grid generator than can use snowflake and X-point target equilibria, among others.
After the grid data is generated from INGRID and the species and surface data is obtained from DivGeo, the third step is to manually convert the surface data obtained from DivGeo and grid data into SOLPS-ready input files using custom python routines. After this third step, all geometrydependent input files should be generated, and the remaining geometry-independent input files can be written in the usual way. The standard preprocessing routines can then be run, followed by submission of the simulation itself until convergence is achieved. Converged SOLPS-ITER simulations have been obtained for a snowflake-minus, snowflake-plus, and for the first time for an XPT geometry. The results of these geometries will be discussed over the following sections.

TCV-like snowflake plus simulations
As part of the validation of this novel workflow, a snowflake plus geometry has been simulated in SOLPS-ITER. The equilibrium used is similar to that of Tokamakà configuration variable (TCV), though the PFCs are unique to this study. An input power and core density of 20 KW and 1.1×10 19 m −3 were used, and the electron temperature profile of the converged solution is shown in figure 3(a). This configuration has been simulated to ensure that in the limit of large separation between the two separatrixes, that the plasma behaves much like a SN. Indeed, the profiles from this figure are well behaved, and in general this configuration is not greatly effected by the presence of the secondary X-point. Moreover, inspecting the plasma fluid velocity profile in figure 3(b), it is clear there is no unexpected or irregular transport near the additional X-point boundaries.

SPARC-like X-point target results
For this study a SPARC lower XPT geometry was simulated to verify the modified code and workflow. In this geometry the two separatrixes are 0.4 mm apart mapped to the outer midplane, and the XPT simulation was compared to a lower SN with a similar plasma current. The XPT grid files were generated using the novel workflow in section 2, whereas the SN input files were generated using the standard CARRE workflow as a benchmark for comparison. PFC contours for both simulations are based upon those of SPARC, but are modified to allow a wider simulation domain. When the relevant input files were obtained, both simulations were run with the same coupled SOLPS-ITER routines. Both grids were roughly matched in terms of radial boundaries; the SN extending from −7.4 mm to +1.7 mm from the separatrix at the outer midplane, and the XPT extending from −7.8 mm to +2.0 mm.
The SN and XPT simulations were pure deuterium simulations without drifts. The same input power of 200 kW and approximately the same outer midplane density of 1.15 × 10 −19 m −3 were used for both simulations. Note that these values are not representative of SPARC, but have been used to access a regime of modest heat and particle flux (where extensive validation of SOLPS-ITER has been performed) without the need for impurities. The transport coefficients were also held constant for both runs, to achieve a ≈0.4 mm heat flux fall off width at the X-point mapped to the midplane and a ≈4 mm density fall off width at the outer midplane. Crucially, because no impurities are used and because both grids have similar total flux expansions, there should be little heat flux dissipation in both SN and XPT simulations. Moreover, the separation of the separatrixes is relatively large at ≈1λ q , so there should not be significant modification of peak heat fluxes due to the secondary X-point. As such, the agreement of these simulations can serve as an effective benchmark for the new code changes and workflow. Run times before convergence was reached was similar for the two configurations, with the SN taking roughly 7 h and the XPT taking 10 h. The higher run time for the XPT is expected due to the higher resolution of grid cells in the outer divertor.
The 2D parallel electron heat flux density profiles of the SN and XPT simulations are shown in figures 4(a) and (b) respectively. From these figures, it is clear to see good overall qualitative agreement between the two configurations, both in profile shape further upstream and the maximum and minimum heat flux values throughout the divertor. Furthermore, the radial profiles of the parallel heat flux density at the divertor entrance can be seen in figure 5(a), plotted as a function of r u , which is the radial distance from the separatrix mapped to the midplane. These upstream profiles show good quantitative and qualitative agreement, as was also found in [15]. In particular, both profiles show a peak parallel heat flux density of 80 M Wm −2 , and a fall-off width of ≈0.4 mm.
Given that the upstream profiles of the two configurations are similar, and that there are no strong heat dissipation mechanisms such as impurities, then the target parallel heat flux density profiles in the two cases should be similar. Indeed, the target parallel heat flux densities are shown in figure 5(b), with the fluxes for both XPT primary targets (defined in figure 4(b)) plotted. The figure shows a similar qualitative profile between the SN and XPT configurations. Moreover, when the target profiles for the two primary targets for the X-point target are overlaid, they form one relatively continuous profile apart from at the secondary X-point at ≈0.4 mm from the primary separatrix. This tends to support the idea that these physics simulations are running with no issues, and that the baseline grid is mapped correctly.  In this attached, low power regime with no impurities, one would expect the target heat flux profiles to be dominated by conduction and radial diffusion. To test this, a simple conduction-diffusion model for target parallel heat flux has been compared to the simulation data for the outer targets. The model assumes the target parallel heat flux density profiles q ||,t are a convolution of an exponential decay at the upstream, with a Gaussian representing heat flux broadening due to radial diffusion [22]: Here q ||,0 is the peak upstream heat flux density, λ q is the heat flux decay width, and S is the width of Gaussian broadening, which are all fitting parameters. R u and R t are the major radii at the target and upstream respectively, and θ is 1 inside the SOL of interest and 0 outside. When describing  The target electron temperature profiles as a function of radial distance from the separatrix (mapped to the outer midplane) for a SPARC single null and XPT divertor. (b) The profiles of connection length, product of connection length and upstream electron heat flux density, and the difference between upstream and downstream temperatures to the power 2/7, plotted as a function of radial distance from the separatrix (mapped to the outer midplane) for a SPARC XPT divertor. SOL splitting in the XPT, the primary and secondary SOL are modelled separately. This model has been fit to the data, and is shown in figure 5, with the two XPT outer targets fitted separately. Indeed, very good agreement is observed with the SOLPS data. In particular, the goodness of fit for the two XPT outer targets indicates the impact of this alternative divertor in this regime can be effectively described by simple models of radial transport.
Finally, the target electron temperature profiles for both simulations are plotted in figure 6(a). Unlike the heat flux profiles, these profiles differ significantly between the two grids. However, this is not unexpected given that ultimately the grids are not identical, and in fact the connection length profiles vary between the two. What is crucial to determine, therefore, is whether these differences in target parameters can be explained by physics analysis. Using simple arguments of the two point model, the target temperature should vary [23]: where L || is the connection length and κ 1 is the Spitzer conductivity divided by T 5/2 , which is approximated as a constant in simple two point model analysis. Consequently, given a constant κ 1 and assuming the heat flux density does not change significantly over the divertor volume, the difference in upstream and downstream temperatures to the power 7/2 should vary with q ||,u L || . The profiles of connection length, the product of connection length and heat flux density, and the difference in upstream and downstream temperatures to the power 7/2 are plotted in figure 6(b). Indeed, this figure shows the XPT simulation temperature profiles are in very good agreement with simple two point model physics. There is small deviation from the trend at R-R sep = 0.25 mm, but this can be attributed to flux limited heat conduction and increased heat flux losses near the secondary separatrix.

MAST-U-like snowflake minus results
The final simulations used for this study are of three MAST-U-like snowflake geometries. They are lower low-field side snowflake minus configurations, with various degrees of X-point separation, as can be seen in figure 7. All three equilibria have been run with identical input powers and core boundary densities of 400 KW and 1 × 10 19 m −3 , along with identical transport coefficients. These equilibria have been used to verify the behavior of snowflakes as X-point separation is varied. One would expect that as the X-point separation increases, the ratio of power incident on the first and third outer target (shown in figure 7) to change, with more power incident on the third. To test these predictions, the fractional power load on each outer target has been plotted as a function of X-point separation for these grids, and is shown in figure 8. These are total loads, taking into account plasma, neutral, and radiation contributions. Here the X-point separation is defined as the radial separation of the primary and secondary separatrix when mapped to the midplane. The figure indeed shows consistency with expectations, as the fractional power to target 3 increases with X-point separation, whilst the target 1 power loading decreases. The second target receives very little power, which is expected given that it receives no parallel transport from the SOL, and any power loads must come from cross-field transport, neutral loading, or radiation.
Finally, when X-point separations are varied in a snowflake, it is important to understand how the entire plasma state varies, not simply target loading. One important phenomenon to study (and one to which SOLPS-ITER lends itself well) is how neutral transport can change with snowflake configuration. To this end, the EIRENE neutral densities (including atomic and molecular) have been plotted for each simulation in figure 7, and are shown in figure 9. In these plots, the neutral density distribution changes as the peak in target ion flux changes from the first to the third target. In particular, the grid with the smallest X-point separation shows lower neutral density in the lower divertor chamber, and higher density in the SOL than the equilibria with larger separations. Though these simulations are not run with sputtered impurities, it is important to note that the separation of the X-points can also affect sputtered impurity content, as having too much sputtered impurities such as carbon in the upper SOL may negatively impact core performance and the pedestal.

Conclusions
Throughout this work, a new way of simulating X-point target and low-field side snowflake divertors in SOLPS-ITER is presented. The SOLPS-ITER source has been edited to accept multiple X-point alternative divertors, whilst still maintaining functionality for standard divertors. INGRID, an interactive Python grid generator has been used to generate the grid data, which is then converted into SOLPS-ITER ready input files using custom python routines. Results and analysis from SOLPS-ITER simulation of a TCV-like snowflake plus, a SPARC-like X-point target, and MAST-U-Like snowflake minus equilibria have been presented, and demonstrate the robustness of the new workflow. Consistency of these simulations with predictions from the two-point model supports the validity of the results presented. Moreover, the results show expected behavior when the X-point separation for a snowflake grid is varied. In future work this workflow should be used to leverage the full complexity of SOLPS-ITER, including drift-activated runs, and full multi-species transport.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.