Investigating the n = 1 and n = 2 error fields in W7-X using the newly accelerated FIELDLINES code

No fusion device can be created without some uncertainty; there is always a slight deviation from the geometric specification. These deviations can add up create a deviation of the magnetic field. This deviation is known as the (magnetic) error field. Correcting these error fields is desired as they cause asymmetries in the divertor loads and can thus cause damage to the device if they grow too large. These error fields can be defined by their toroidal (n) and poloidal number (m). The correction of the n = 1 and n = 2 fields in Wendelstein 7-X (W7-X) is investigated in this work. This investigation focuses on field line diffusion to the divertor, a proxy for divertor heat flux. Such work leverages the 25× speedup obtained through the implementation of a new particle-wall collision model. The n = 1 and n = 2 error fields of the as-built coils model of W7-X are corrected by scanning phase and amplitude of the trim and control coils. Reductions in the divertor load asymmetry by factors of four are demonstrated using error field correction. It is found that the as-built coils model has a significantly lower m/n=1/1 error field than found in experiments (Bozhenkov et al 2018 Nucl. Fusion 59 026004).


Introduction
Steady-state operation is one of the main goals of Wendelstein 7-X (W7-X). A key aspect in this is the divertor performance which ensures proper heat and particle exhaust [1,2]. W7-X is equipped with an island divertor consisting of five upper and lower divertor modules. This island divertor makes use of the magnetic islands at the plasma edge [3][4][5]. The flip-mirror symmetry, about each of the five magnetic field * 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. periods in W7-X (the so-called stellarator symmetry), implies that the divertor modules have similar symmetry. However, E × B drift causes up-down symmetry to not be guaranteed at each toroidal cut [6]. These divertor plates intersect edge magnetic islands which reside at radial positions with low-order rational values of the rotational transform ι 2π = n/m. Edge rotational transforms varies between m/n = 4/5 and 6/5 in W7-X [7]. Configurations with m/n = 5/5 are usually favored as the pumping and divertor shape has been optimized for this configuration [8,9].
The magnetic fields of W7-X have been tailored to spread heat evenly across all ten divertors, with wetted areas of around 1.5 m 2 [10,11]. The target plates in W7-X are designed to withstand a load of 10 MW m −2 in steady-state. However, error fields cause asymmetry between the heat loads on the divertor modules. The head load asymmetries can arise from many factors, such as coil movement under load, component misplacement, and others. Error fields can modify the power distribution over the divertors, resulting in higher power loads in some modules and lower in others. In the worst case, other in-vessel components can even be damaged due to unwanted heat loads [9,12]. Therefore, the W7-X magnetic field error must be kept in check by creating a compensating field.
Error fields as minor as 10 −4 (normalized to average magnetic field strength on the magnetic axis) modify the edge topology enough to add additional undesired islands [9]. In those magnetic configurations with ι 2π = 1 at the edge, this resonance dominates the plasma-wall interaction of the island divertor. Modes that are resonant with this surface (m/n = 1/1, 2/2, 3/3, 4/4) are of the greatest interest [9]. It can be shown that each of these perturbations falls off in the radial direction according to the power of the poloidal mode number. This makes the n = 1 and n = 2 modes the most significant and the focus of this work.
The influence of error fields was taken into account during the construction of W7-X. Efforts were made to minimize the error fields. Coil manufacturing and significant assembly steps were monitored by metrology [13,14], followed by the module adjustment to compensate for construction errors accumulated. The result is a normalized n = 1 error component of below 0.3 × 10 −4 [15].
W7-X started with divertor experiments in the second half of 2017, where investigation of error fields was continued [16,17]. The 1/1 error field was investigated in two ways. The first was by measuring the shift of the helical axis and the second involved canceling the intrinsic 5/5 islands by divertor control coils [9]. These two resulted in a normalized 1/1 error field amplitude of (0.8 ± 0.1) × 10 −4 and (0.5 ± 0.2) × 10 −4 with a phase of around 150 • [9]. It was shown that the 1/1 error field could be compensated using the trim coils in W7-X. The trim coils are designed to be used for this purpose [18]. Only around 100 A are necessary to achieve this, significantly below the 2000 A maximum capability of the trim coils. Having compensated the 1/1 error field, it was also possible to identify the normalized 2/2 error field magnitude at around 0.6 × 10 −4 [9].
Another asymmetry was observed in W7-X during these experimental campaigns involving the strike line location on the divertor modules. The distance from the inboard side (the pumping gap) to the strike line differed between the divertor modules. Symmetric strike line locations ensure similar levels of neutral compression. This will help to achieve the desired particle exhaust targets. The strike line location is dictated by the magnetic island near the divertor module. Thus, an asymmetry in the strike line location indicates that the magnetic configuration is not symmetric. Figure 1 visualizes how the strike line arises from the interaction between the 5/5 island and the divertor plate. The radial position of the magnetic islands can be shifted by changing the rotational transform using the planar superconducting coils. However, this works the same on all divertor modules and thus cannot correct asymmetry [8]. For that, the trim coils need to be used.
A great deal of the effort on error fields in W7-X has focused on experimental analysis. Simulations of divertor loads via field line diffusion (FLD) have been limited given the computational demands of such simulations [19,20]. In this work, we leverage a speedup in the computation of particlewall interactions to perform detailed FLD calculations with error fields in W7-X using the FIELDLINES code.
The following section will discuss the FIELDLINES code and the improvements that have been done to it in detail. Section 3 will discuss the methodology of investigating the error fields with the FIELDLINES code. Section 4 will discuss the results obtained and will compare them to the experimentally obtained observations. The work will end with a conclusion and discussion in section 5.

FIELDLINES
The FIELDLINES code is designed to follow field lines as they pass through any toroidal magnetic confinement fusion device, such as a tokamak or stellarator. The code itself is written in FORTRAN and parallelized using the Message Parsing Interface. The code is interfaced to a wide variety of magnetic field sources, including vacuum (coils and MGRID files), plasma (VMEC [21], HINT [22], or EQDSK [23]), and analytic magnetic field sources. It should be noted that in the case of this work, a W7-X coils geometry data and analytic error field are utilized.
The fundamental purpose of the FIELDLINES code is to follow field lines in a toroidal domain [17]. This is done by solving the following ordinary differential equations: This representation allows the field line to be parameterized by toroidal angle. FIELDLINES offers several ways of solving this ODE. This work will use LSODE to solve the ODE [24]. No plasma model is present in FIELDLINES, as it is not needed to follow the field lines.
FIELDLINES can also perform FLD to model divertor strike locations. New field lines are initialized on the original field lines calculated without diffusion. Then, diffusion is applied to the new field lines at regular intervals along the toroidal direction. The diffusion is done with a displacement related to the magnetic field and an input 'diffusion coefficient'. Diffusion in the toroidal direction is ignored, a similar assumption to assuming only radial diffusion. Since the diffusion is applied at regular intervals along the field line, the unit of the input diffusion coefficient is m 2 m −1 , instead of the usual m 2 s −1 .
In short, at each interval, a displacement is applied where δ R is the radial displacement, n r a normally distributed random value around a mean of one, and D the diffusion coefficient in m 2 m −1 [12]. This implementation of diffusion has been verified, and it has been shown that an input diffusion coefficient of µ = 1 × 10 −6 m 2 m −1 aligns with the experimentally observed diffusion [12]. The conversion factor between m 2 m −1 and m 2 s −1 is the stereotypical velocity of a particle, which in the case of W7-X is ∼1 × 10 6 m s −1 . Each new field line initialized on the original field lines carries an equal fraction of the power flowing through the plasma boundary. As a result, simulation results scale to different possible plasma heating powers. This approach implies that the full heating power is deposited on the targets and the first wall. If a part of the heating power is radiated, only the remaining fraction must be used in the FLD. The result of FLD is that the newly initialized field lines hit something or reach the end of the simulation. FLD is a good proxy for the heat load as it approximates where particles would deposit their energy [12]. A wall model is required to be able to check where FLD hits will be located. The objective of this wall model is to check for collisions between the field lines and the wall shape supplied to FIELDLINES. This is also called the wall collision module. Within the wall collision module, the wall is discretized as a mesh of triangular faces. Whenever a collision has to be checked, a line segment is given to the module. This line segment can be a part of a field line or a particle path. The module returned whether or not it was a hit and, if so, where.
Finding the intersection between a line segment and a triangle is a common geometric problem, and well-defined solutions exist. One common approach is to use barycentric coordinates. In this case, the first step is to find when the line segment intersects the triangle plane. This can be formulated as: with N the normal of the triangle, O the origin of the line segment, R its direction, and D the distance from the origin to the triangle plane. The result is t, the distance along the line segment where it hits the plane, normalized to its length [25]. If t is between 0 and 1, it will hit the plane on this line segment. Only then a secondary check is done whether or not this intersection falls within the triangle itself. This can be done by calculating two parameters: where x ′ is the distance from the origin of the line to the first vertex of the triangle. e 2 and e 1 are the two edges spanning the triangle starting from the first vertex. If both α and β are above zero and their sum is below 1, the line segment will hit the plane of the triangle inside the triangle itself [25]. This problem is also common in codes; for example, it is also used for ray tracing. Several algorithms exist, such as the Badouel algorithm [25], the Möller-Trumbore algorithm [26], and Segura's algorithm [27]. The current implementation within the wall collision module uses an adapted version of the algorithm of Badouel. The adaptation uses as many precalculation steps as possible to reduce computational load during the main loop. All the algorithms mentioned above try to solve the geometric problem using as few operations as possible. The optimal algorithm depends on the application [28].
The previous implementation of the wall collision module checked all triangles to see if a hit occurred. However, the number of triangles grows as the mesh gets physically larger, more complex shapes are added to the mesh, or the resolution of the mesh is increased. Thus, most meshes result in a high number of triangles and thus result in high computational costs for the wall collision module.
This work uses a sped up version of the wall collision module, speeding up FIELDLINES and the other codes within STELLOPT. An acceleration structure has been used to create this speed up. In short, a performance increase of over 25 times has been found. The new acceleration structure scales very well with the number of faces in a wall mesh. In contrast, the old implementation scaled nearly linearly. More detail on how the acceleration structure was implemented and tested can be found in appendix A.

Methodology
This work explores the correction of the n = 1 and n = 2 error fields using FLD simulations and the coil model of W7-X. Self-applied analytic error fields on top of the ideal field and the error fields arising from the 'as-built' electromagnetically loaded coil are considered. The trim and control coils are then used to probe and correct magnetic fields, as is done in experiments. This probing is done by varying the amplitude and phases of the currents in the coils. This is the compass scan technique [29]. The trim coil and control coils can be seen for a single divertor module in figure 2. Divertor strike lines, as modeled with FLD, are examined in terms of their load asymmetry and strike line shape.
Simulations are performed using the same settings as described above for the performance test. However, the magnetic fields are now obtained from three sources: the ideal coils set, the ideal coils set plus an analytic error field, and an as-built electromagnetically loaded coil model. The ideal coil set represents the coil imagined by physicists and engineers, employing a flip-mirror 5-fold symmetry. An analytic error field of the form has been implemented in the code. B err is the error field amplitude in Tesla, n is the mode number, and ϕ is the toroidal angle in radians. B R is the applied error field and B ϕ arises from the requirement that ∇ · ⃗ B = 0. The choice to have B ϕ modulate over B Z , avoids radial shifts of the plasma due to B Z . The as-built coil model is based on detailed metrology of the actual coils of W7-X along with ANSYS modeling of motion of these coils due to pre-load, cooling, and electromagnetic loads. This coil model contains any intrinsic errors related to superconducting coil manufacture and device assembly. In all cases, ideal trim and control coil geometry is assumed. We begin with testing the trim and control coils' ability to compensate for our analytic error field. This compensation is done by applying an n = 1 error field to an ideal model of W7-X. The phase and amplitude of this error field can be chosen at will. A sweep of different trim coil currents and phases is done to compensate for this self-applied error field. The trim coils are energized with a n = 1 pattern, a single cycle of a sine wave.
The success of the compensation is measured by the asymmetries of the heat loads and the strike line locations. For the heat load, the asymmetry is measured by the relative standard deviation of the number of FLD hits on the different divertors of W7-X. The relative standard deviation does not have a unit, as it is the standard deviation divided by the mean FLD hits over all the divertors. Thus, in the case that all divertors get the same load, the asymmetry is zero. The FLD hits on the divertor are also spread out in the radial direction. The center of this spread is determined for each divertor. The standard deviation of these centers for each of the five divertors, as measured in mm, is taken as the measure for asymmetry for the strike line location. These asymmetry values are calculated for the upper and lower divertors separately. As mentioned, up-down symmetry is not guaranteed. This results in four asymmetry values, the heat load and strike line location asymmetry values for the upper and lower divertors, respectively.
After finding the asymmetries for each trim coil current setting with the self-applied error field, a fit in the 2D plane spanned by the amplitude and phase of the applied current can be done. This is done for four asymmetry values separately. These fits will result in the optimal trim coil settings for this error field. The simulation is deemed successful if two checks are passed. First, the individual simulations should result in values for the asymmetries that align with the expectations based on the trim coil current. For example, the asymmetries should be maximized if the trim coils drive an n = 1 field inphase with the applied error field. Second, the optimal coil settings as found by the fit should align with expectations. This process is repeated for a simulation with self-applied n = 1 and n = 2 fields.
The applied error field in this work models a specific toroidal mode number (n) and a broad spectrum of poloidal modes (m). As a result, the applied error field contains both resonant and non-resonant modes. The exact nature of these modes would require decomposition of the applied field, which is a function of the said error field and the underlying magnetic geometry. Such decomposition is beyond the scope of this work, as we sought only to demonstrate compensation of the effect of error fields on the divertor loads.
Similar simulations are performed using the as-built coils model of W7-X. In this case, no analytical error field is applied as the error fields arise from the superconducting coils. With this coil model, another sweep is done using different trim coil currents to balance heat loads and strike line locations on each of the divertors. This result can be compared to the work of Lazerson et al and Bozhenkov et al who investigated the strike line locations and heat load asymmetry respectively experimentally [8,9]. A sweep of the control coils will also be done to reduce the found asymmetries further. The control coils are energized in a n = 2 pattern, two cycles of a sine wave.

Results
The simulation results can be split into two parts. In the first part, we discuss simulations with the self-applied analytic error field. In the second part, we focus on the as-built coil model simulations.

Self-applied error field
The magnitude of the self-applied n = 1 error field is chosen to be 2 × 10 −4 T. This is based upon the 0.8 × 10 −4 normalized value found by Bozhenkov et al [9] combined with the magnetic field strength of W7-X (2.5 T). The phase is set at 36 • .
A simulation is performed with no applied trim or control coil currents to investigate the impact of the error field. As can be seen in figure 3, a distinct asymmetry is found in the divertor heat loads. This aligns with the observed experimental result for a 1/1 error field [9]. The measured heat load asymmetry is 0.31 for both the upper and lower divertor. Minimal asymmetry is observed in the strike line location, as visible in figure 4. The only apparent difference is the number of particles striking the divertors.
A broad sweep of different trim coil currents is done to compensate for the n = 1 error field. The currents used are 50, 100, 150, 200, and 300 A, based upon the experimentally determined optimum of around 100 A. The phases used are 0 • -334 • in steps of 36 • . This results in 50 total data points. Extra data points are gathered around the found optimum after an initial fit was done. The results of this can be seen in figures 5 and 6. Figure 5 shows the effect of the trim coil currents on the heat load asymmetry. The left figure shows the upper divertors, while the right figure shows the asymmetry in the lower   at 210 • . The lower divertor fit results in 113 A at 212 • . These values align with the experimentally found results of around 100 A for an error field magnitude of 0.8 × 10 −4 [9]. As one would expect, the phase found has a 180 • phase shift compared to the input n = 1 error field. Running additional simulations with these values shows that this significantly improves the particle load balance for the divertors. The upper divertor asymmetry is improved from 0.306 to 0.013; the lower divertor asymmetry goes from 0.311 to 0.011; both are reduced by over 95%. Figure 6 shows the asymmetry in the strike line location for varying trim coil currents while applying the n = 1 error field. Low values of asymmetry are found for very different phases at the same current. Also, the convex fit does not seem to align very well with the data. Thus, the sweep of trim coil currents does not result in an optimum for strike line location asymmetry. Also, the strike line location asymmetry is minute; the strike line asymmetry increases by only 0.9 mm after applying the n = 1 error field. During the sweep, a maximum asymmetry of only 12.7 mm is found at 300 A. The strike line movement observed in experiments is significantly higher, peaking at ∼5 cm [8]. This indicates that the n = 1 error field does not cause this asymmetry. Another effect must be present in the device. This is checked by running simulations with 500 A or 1000 A currents through the trim coils. In these simulations, strike line location asymmetry of >12 cm is visible, showing that a trim coil induced n = 1 field can cause strike line location asymmetry. However, the current used in experiments to compensate for the n = 1 field is just over 100 A. This again supports that the n = 1 error field at its measured size does not seem to cause the strike line location asymmetry. Figure 7 depicts Poincaré plots of the outboard magnetic island (ϕ = 0 • ) for the ideal coil set, the ideal coil set plus a self-applied n = 1 error field, and best parameters for trim coil compensation for this error field. This is done to check if the applied error field corrections restore the original configuration. The entire Poincaré plot is investigated, but a focus is put on the 5/5 magnetic island. Figure 7 shows a zoomin of Poincaré plots of the trim coil and error field combinations. The 5/5 island of the ideal coils is presented on the left. The magnetic configuration around the 5/5 island of the ideal coils with the self-applied n = 1 error field is shown in the middle. The center of the 5/5 island has shifted up in zdirection, a clear signal of an n = 1 error field. On the right, the zoomed-in Poincaré plot is shown after applying compensating trim coil currents. The plot clearly shows that the magnetic configuration is restored, as the island has moved back to around z = 0. The compensation has worked. Poincaré plots are made for different angles in the toroidal plane of W7-X, as error fields can cause a difference depending on the phase. Other plots showed similar results to the presented plots at ϕ = 0 • . Three more sets of simulations are done. First, a simulation was done with only an n = 2 error field. The control coils are swept to check the ability of the control coils to compensate for an n = 2 error field. Second, a set of simulations is done with both n = 1 and n = 2 error fields. The previously found optimum value for the trim coils is applied, and the control coils are swept. The third set of simulations is where no trim coil current is applied while the control coils are swept with both n = 1 and n = 2 error fields applied. In all these cases, the n = 2 error field is given a magnitude of 1.5 × 10 −4 T, again based on Bozhenkov et al [9]. Its phase is set to be 144 • to allow separation from the n = 1 field. The result of these three sweeps is summarized table 1. Table 1 shows that compensating the n = 2 error field can be done in FIELDLINES. It is required to have compensated the n = 1 field with the trim coils; otherwise, this field dominates the heat load asymmetry. This is visible in the last row, where a sweep of the control coils is done without trim coil current. This does not successfully compensate for the applied error. The three cases show that it is possible to compensate for selfapplied error fields.
The phase needed in the control coils to compensate is not 180 • opposite of the applied phase to the n = 2 error field. The optimal phase found for both sets of divertors differs significantly. Only one optimum is found to reduce the effect of the applied n = 2 field, not two with opposite phases, as might be expected. This hints at other harmonics influencing the compensation of the n = 2 error field. It could be that control coils do not drive the clean n = 2 field and thus induce other error fields. The coupling to the B ϕ due to the requirement that ∇ · ⃗ B = 0 could also cause a single optimum to appear instead of two optima. Regardless, the current amplitude for the n = 2 field compensating is within the expectations.
The strike line asymmetry is left out of table 1 as none of the self-applied error fields induced a significant asymmetry in the strike line location.
To conclude, the FIELDLINES code can simulate n = 1 and n = 2 error fields, and trim coil and control coil current sweeps can be used to compensate for self-applied error fields. Optima that align with the expected phase and amplitude of a compensating field were found. Finally, it was demonstrated that the n = 1 error field must first be neutralized to compensate the n = 2 field adequately. A similar approach was taken for the experimental error fields [9].  Goodness of fit is indicated by given R 2 .

As-built coils
Having demonstrated that we can analyze and compensate known n = 1 and n = 2 error fields, we apply our analysis to the 'as-built' coils model. This model incorporates error fields arising from manufacturing, assembly, and electromagnetic deformations of the superconducting coils. A broad sweep of trim coil currents is performed to probe for the present but unquantified n = 1 error field. After finding an optimum, a sweep of the control coils reduces the asymmetry even more. The results of the trim coil sweep can be seen in figures 8 and 9. A complete overview of both sweeps can be seen in table 2.
The trim coils can compensate for the present n = 1 field and reduce the asymmetry in the heat flux by around 60%. As can be seen in figure 8, the optimum is different for the lower and upper divertor. The upper divertor results in a current of 38 A at 22 • . The lower divertor results in 30 A at 55 • . Good agreement with the fit is also found, with a R 2 of around 0.9 for both divertors. These values for the current are lower than expected based on the results found by Bozhenkov et al where over 100 A was found [9]. The 22 • -55 • phase matches the result found by Bozhenkov et al who found a phase of around 150 • for the error field, as this is opposite to the phase found for the compensation in this work. However, the significantly lower current still indicates that the as-built coils model does not reflect the full extent of the error fields in W7-X or has inaccuracies. The significant difference between the optima for the upper and lower divertor can be explained by the lack of up-down symmetry in each toroidal cut. Higher-order error fields could also affect this difference.
No clear optimum is found for the strike line location asymmetry during the trim coil sweep as can be seen in figure 9. The best value did reduce the asymmetry. However, no clear pattern was visible to determine which trim coil current should be used to reduce strike line location asymmetry. This can also be seen in the R 2 of the fit, which has dropped to around 0.5, showing that the pattern cannot be explained with an elliptical fit. Just as before, the strike line location asymmetry is minimal, smaller than observed in experiments. This concludes that the as-built coils model does not include the effect that causes the strike line location asymmetry. The strike line location asymmetry observed could be caused by other factors Figure 9. Strike line location asymmetry for different trim coil currents energized in a n = 1 fashion for the 'as-built' coils model. Asymmetry measured as standard deviation (mm). Color of dots indicates results simulation. Background is fit. Cross indicates optimum of fit. Goodness of fit is indicated by given R 2 . such as divertor misalignment or non-ideal errors in other systems. Another option is that the as-built coil does include the effect that causes the strike line location asymmetry, but that this effect is underestimated. This would align with the underestimation of the error fields causing the heat load asymmetry. A current of 35 A at 36 • is applied to the trim coils during the control coils sweep. This value is an average of the previously found optima. The sweep of the control coils again reduces the found heat load asymmetry, resulting in a further 40% decrease for the upper divertor and a 24% decrease for the lower divertor. However, the optima for the lower and upper divertor are at entirely different locations. The upper divertor requires nearly twice the current, and there is a 130 • phase difference between the two phases. This is again likely caused by higher-order effects or the control coils inducing other harmonics.
The optimal current found for the control coils is higher than the current for the trim coils. However, no conclusion can be drawn on the size of the n = 2 field compared to the n = 1. The trim and control coils are not placed in the same position, and thus their currents do not relate to the same magnetic field strength in the plasma. To conclude, the impact of the coils on the n = 2 field requires further investigation, and indications of higher-order fields are found. The results also indicate that other factors than the coils mainly cause the n = 1 error field.

Conclusion
In this work, we present a newly improved version of FIELDLINES, achieving a speedup of between 25 times and 1000 times. This speedup was put to use by investigating the n = 1 and n = 2 error fields in W7-X. The ability to correct error fields with the trim and control coils was demonstrated by applying an analytic error field. The heat load asymmetry was reduced by compensating for the n = 1 error field using the trim coils. No significant asymmetry in the divertor strike line location could be detected for the applied n = 1 error field. An n = 2 field was also applied and compensated for successfully using the control coils. An unexpected phase difference between the applied n = 2 and the required control coil current was found. A potential cause for this could be coupling between B R and B ϕ .
The heat load asymmetry in the 'as-built' coils model was successfully compensated using applied n = 1 trim coil phasing and n = 2 control coil phasing. However, the optimal value for the trim coil current does not match the experimentally determined value and is significantly lower. This indicates that the n = 1 error field in the 'as-built' coils model is significantly smaller than the experimentally observed 1/1 field. Also, a significantly different control coil current is needed for the upper divertor and the lower divertor, indicating higher-order error fields. Additionally, the 'as-built' coil model does not show the significant strike line position asymmetry measured in experiments. Possible causes for these asymmetries in experiments could be divertor misalignment or the non-ideal performance of other magnetic components. Drift effects, while not included in the model, are not considered an optional cause as they cannot induce toroidal asymmetry [6].
Several options for future work exist. First, higher-order and poloidal self-applied error fields could be investigated in FIELDLINES. This could help the understanding of the impact of the error fields on the divertor loads. It could also explain the unexpected control coil phase obtained. Second, higherorder and poloidal error fields could also be investigated for the as-built coils model to apply this knowledge to practice. Lastly, an option would be to investigate other non-ideal parts of W7-X to find the cause of the n = 1 error field and the strike line asymmetry. Creating models for these and simulating these with FIELDILNES could allow them to be investigated similarly.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors. collision module As mentioned, an acceleration structure has been used to achieve the obtained speed up. Acceleration structures have been around in the field of line-triangle intersection detection since the 1980s. The concept itself is simple: instead of checking every single triangle for a possible intersection, design a system that allows most triangles to be skipped [30,31].
This simple concept has resulted in many different types of acceleration structures. Concepts include uniform grids, bounding hierarchy volumes, kd-trees, and octrees [32][33][34]. All concepts work in a similar method. They introduce a spatial decomposition (called blocks in this work) and first check in which block the line is present. In that case, only the triangles of that block have to be checked for line-triangle intersections. Another option is that the line spans several blocks, resulting in multiple blocks having to be checked. Either way, this allows for a significant reduction in the number of triangles to be checked. A uniform grid acceleration structure is implemented. Its advantages are its simplicity, easy mesh creation, and easy traversal between the blocks. Also, a simple method is chosen as all acceleration structures show similar performance levels compared to using no acceleration structure.
The resulting code is tested by running a typical FIELDLINES run. These runs are performed on 16 nodes of the MPCDF Draco cluster. Each node features 32 Intel 'Haswell' Xeon E5-2698v3 cores. During this run, the focus is on vacuum magnetic fields. The simulation domain has a 128 major radius, 176 toroidal, and 128 vertical gridpoints over which C1 splines are constructed for the quantities RB R /B ϕ and RB Z /B ϕ . This choice of toroidal gridding ensures that a toroidal grid point is placed at each of the ten stellarator up/down symmetric planes in the simulation. This preserves the underlying magnetic field structure. The radial grid runs from R = [4.25, 7.25], the vertical grid from Z = [−1.75, 1.75] and the toroidal grid from ϕ = [0, 2π]. The simulations must be performed over the entire toroidal domain as error fields break the five-fold field periodicity inherent in the stellarator. A diffusion coefficient µ = 1 × 10 −6 m 2 m −1 is chosen. As mentioned, this value is consistent with previous simulations and has been shown to have excellent agreement with the diffusion in W7-X [12]. FLD simulations begin with ∼256 000 particle initialized along a vacuum flux surface near the edge of the confined plasma volume. Such data is obtained using the ideal unperturbed coil set. These particles are followed once parallel to the field lines and once anti-parallel. The 'standard' configuration of W7-X is used for all simulations in this work.
In this test, the improved code shows a performance increase of over 25 times. It has also shown that the new acceleration structure scales very well with the number of faces in a wall mesh. In contrast, the old implementation scaled nearly linearly. For more complex meshes, the speedup can be up to a hundredfold. This is due to the uniform grid ensuring fewer face checks have to be done for each line segment. Figure A1 shows how the computational time drops as the block of the uniform grid is made smaller. The different lines represent increasingly complex meshes with more faces in them for every refinement. The geometry is the same; the faces are made smaller. This method ensures that all other aspects are equal. The x-axis shows an increasing number of blocks (and thus smaller blocks) from left to right. The figure clearly shows that increasing the number of blocks reduces the computational time for every mesh refinement. The improved code is also validated by comparing the physical results to the results of the old version. Figure A2 shows how the scaling with the number of faces in the wall mesh is no longer linear for the same simulation. The scaling improves as the number of blocks is increased. The x-axis shows the increasing number of faces in the mesh. This again is done by making the mesh faces smaller with the same geometry. Each line is a different block size, resulting in a different number of blocks. The figure again shows that having a small block size gives better performance. This holds for any mesh, regardless of the number of faces in that mesh. It is observed that keeping the ratio of faces to blocks below 100 ensures the best performance while preventing an excessive number of blocks. A limit on the number of blocks has not been found in the tests done.