Anisotropic diffusion as a proxy model for the estimation of heat-loads on plasma-facing components

To facilitate the estimation of heat loads on plasma-facing components in fusion devices in various different magnetic geometries, a heat load proxy model was developed based on anisotropic diffusion. In this work, this model is compared to the so-called field-line diffusion approach. To facilitate the evaluation of these models, a novel synthetic camera-based approach for obtaining heat load distributions from Monte Carlo samples was also developed and implemented. With the assistance of this synthetic camera, heat load predictions for the Wendelstein 7-X divertor were obtained and compared with infrared camera observations. It was found that the anisotropic diffusion-based model achieved a closer match to infrared camera observations, while still being suitable in computational effort for large magnetic configuration database scans.


Introduction
Monte Carlo models are a popular approach to simulate the heat-and particle-transport in magnetically confined plasmas. This can mainly be attributed to their excellent ability to handle the differences between transport parallel and perpendicular to the magnetic field lines, which can be multiple orders of magnitude apart, as well as their excellent parallelizability. Furthermore, the parallel and perpendicular transport processes can be modeled in separate subroutines of model implementations, which increases the quality and maintainability of transport codes significantly. An important application of these models is the simulation of heat fluxes onto plasmafacing components.
Of particular relevance are two basic linear transport processes, the diffusion and convection process. These two provide straight-forward mappings between fluid-and Monte-Carlo-views of the same physics. On the one hand, the convection process, which is connected to a first-order fluid equation, can be implemented using single-directional random walks (e.g. exponentially distributed or fixed-distance movement). The diffusion process, on the other hand, is related to a parabolic second-order fluid equation and, on the microscopic level, maps to bi-directional random walks, such as Brownian motions.
As mentioned before, diffusive and convective Monte Carlo models form important basic building blocks for the modeling of heat-and particle transport in fusion plasmas. Of particular importance in the plasma edge are the Braginskii equations, which are, at their heart, non-linear diffusion equations. However, the nonlinearity in these equations is numerically challenging to treat. Common Monte Carlo implementations trace test particles individually without consideration of the direct interaction with other tracer particles, which is a concession made to limit code complexity (exceptions to this tendency include e.g. particle-in-cell methods). The non-linear interactions need to then be modeled by treating the particle statistics (densities, velocities etc) as fields that the tracer particles can interact with. This approach requires an alternating loop between test particle tracing and field accumulation, which substantially increases the required execution time.
However, the full physics picture provided by a nonlinear model is not always necessary to derive the information requested. Sometimes, it is sufficient to have a reasonably well working model which can be used to perform studies on physics parameters not directly related to the transport itself. The discussion presented here is aimed at such an attempt of proxy modeling.

Anisotropic diffusion as proxy model
The modeling approaches presented in this discussion primarily aim to support the near-and mid-term future divertor operation and design of Wendelstein 7-X [1,2]. In the upcoming campaigns, this stellarator device will aim to demonstrate quasi-steady-state operation in high-performance regimes. However, its island divertor configuration can be sensitive to changes in the magnetic topology related to the toroidal plasma current and the plasma beta [3,4]. To assess these changes (and to prepare future divertor adjustment research), the underlying heat-load model needs to combine an acceptable execution time, flexibility in the magnetic background field and reasonable accuracy of the heat-load predictions.
Prior to the developments presented here, there were two main models in use for edge physics modeling at Wendelstein 7-X. The first one is the widely-known and -employed edge simulation code EMC3-EIRENE [5,6]. This code implements a full Monte Carlo model for the Braginskii equations. To compensate for the high computational demands of such a nonlinear simulation, EMC3-EIRENE uses a precomputed interpolation grid to accelerate field line tracing. Unfortunately, the process to generate these grids has not yet been fully automated, which precluded application of this code to our large equilibrium database waiting to be analyzed.
The second pre-existing heat-load model is provided as part of a field line tracing webservice [7], which can optionally enable a perpendicular diffusion process inside its tracing routine. This implements a hybrid transport model, combining a convective parallel transport with a diffusive perpendicular transport process. The convective parallel transport makes this code extremely fast, as particles quickly hit the divertor once they leave the confined plasma. This model was initially used to study the effects of plasma-beta and plasmacurrent related magnetic topology changes on the divertor heat loads. However, its predictions were found to lack some of the heat load patterns observed in experiment [8] (see section 3). Particularly, field-line diffusion failed to reproduce heat load patterns created by counter-streaming heat flows. In such a counter-streaming flow scenario, the heat first flows into one direction parallel to the magnetic field (either with or against the magnetic field vector), and then reverses direction at least once.
As a response, we decided to improve upon the second approach by replacing the convective parallel transport process with a diffusive one. Internally, this was realized by replacing the single-direction field-line tracing with a bidirectional Gaussian random walk. The sign of the sampled tracing length would then be used to determine the tracing direction. The resulting model represents a middle ground between the highfidelity physics modeling of EMC3 and the particularly fast computation of field line diffusion. Furthermore, such a model is still based on simple field line tracing. Therefore, it does not require a field-aligned grid, which requires substantial manual labor to set up.

Field line diffusion
The field line diffusion model is an extension of a field line tracing process. Instead of simply following the field line until hitting a plasma-facing component, the test particle follows the field line for a random distance ∆l, after which it is displaced by a random displacement ∆⃗ x ⊥ perpendicular to the magnetic field (with a randomly chosen direction). These quantities are sampled from their respective distributions as: with a convective streaming velocity v, a perpendicular diffusion coefficient D ⊥ and a mean free path length λ. For small step sizes, this model implements a Monte Carlo process for a hybrid convective-diffusive transport. In a short mean-freepath limit (continuous fluid equation), the tracer particle density follows the process: where ⃗ B is the magnetic field,P ⊥ is an operator representing projection perpendicular to the magnetic field, ⃗ J is the tracer particle flux, n the tracer particle density and f is an inhomogenous source term representing particle seeding or heating. This type of model is very fast to solve (requiring about 4 core hours for 10 5 test particles) and simple to implement, as it requires only a small extension of a field line tracer and scales linearly with the connection lengths in the scrape-off layer. Unfortunately, this type of model is not able to capture counter-streaming flow geometries. In such geometries, particles stream in the direction of the magnetic field in some regions, and against the magnetic field direction in others. These effects are uncommon in axisymmetric fields, but can lead to additional strike patterns in 3D fields as they are found in Stellarators.

Anisotropic diffusion
Anisotropic diffusion can itself be again viewed as a simple extension of the field line tracing process. The key difference to field-line diffusion, however, is that the parallel transport model is also chosen to be a diffusive one. For a random walk-based implementation, this implies that the parallel random walk must be a bidirectional random walk that can step in both forward and backward direction. The natural choice here is to implement both the parallel and the perpendicular random walks as Gaussian random walks approximating Brownian motions. For simplicity, the model presented here simply considers a parallel and an isotropic (omni-directional) random walk. The movement distance ∆l and the perpendicular displacements ∆⃗ x follow now the distributions: with a parallel diffusion coefficient D ∥ , an isotropic diffusion coefficient D and a time step size δt. For infinitesimal time steps, such a model can be interpreted as solving an equation similar to equation (3): This gives rise to two equivalent interpretations, a parallel diffusion with coefficient D ∥ and an isotropic diffusion with coefficient D, or a parallel diffusion with coefficient D ∥ + D, a perpendicular diffusion with coefficient D, and the source term f.

Parameter dependency of flux distributions and computing time
As the fluid transport equations represent the small-step-size limits of the Monte Carlo processes, both models contain scaling parameters that do not appear in their respective fluid interpretation. For the field-line diffusion, this parameter is the mean-free path λ, while the time step size δt is its equivalent for the anisotropic diffusion model. These parameters represent computational accuracy tradeoffs. Small values lead to more physically accurate solutions of the fluid equations, but also to increased computational time. This is particularly important for the anisotropic diffusion model, as the effective distance traversed by a diffusing particle does not scale with the time, but with its square root σ ∝ √ Dt. This in turn implies, that such a model's execution time scales with the square of the connection length in the edge. To combine upstream accuracy with acceptable computation times, we chose to linearly increase the time step size with the step number in the random walk. This curbs the computational cost at the expense of potentially missing some flow reversal points in the downstream region.
In addition to the above parameters, there is another important simplification to be considered. The purpose of the presented models is to calculate the distribution of attached heat loads in steady-state conditions. Therefore, only the fluxes ⃗ J in steady state cases dn dt = 0 are of practical relevance. In such a case, it is possible to linearly rescale all velocities and diffusion coefficients by an identical factor c, which will modify only the tracer particle densities n by a factor of c −1 without modifying the flux distributions. Therefore, the heat loads do not depend strictly on the transport parameters, but only on their ratios. The characteristic quantities for steady state heat load solutions are therefore: • D ⊥ /v for the field line diffusion model, empirically fixed to the value 3 · 10 −7 m by matching to infrared camera observations; • D /D ∥ for the anisotropic diffusion process, empirically fixed to the value 10 −7 by matching to infrared camera observations.
To limit the sampling noise of the Monte Carlo model to an acceptable level while attaining high resolution, we chose to trace 10 5 test particles for each model. For the anisotropic diffusion model, this calculation was done in parallel on the JURECA supercomputer [9], while the field-line diffusion model was traced in a serial implementation [7]. As should be expected from a tradeoff for increased accuracy (due to worse computation time scaling with connection length), the total CPU consumption of the anisotropic diffusion model already substantially eclipses its simpler counterpart. Whereas the field-line diffusion calculation takes roughly 4 h on a single core, the anisotropic diffusion calculation requires a wall clock time of about 1 h on a 128 CPU setup (leaving the anisotropic diffusion model at a substantially more expensive 128 corehours). For our purpose, this is tolerable, but for applications without access to a parallel machine, this might already present a substantial hurdle.

Heat-load density estimation from Monte-Carlo samples
A common challenge that Monte-Carlo methods face is the back-conversion of test particle distributions into fluid parameter distributions, a problem commonly known in statistics and machine-learning as 'kernel density estimation'. Volumetric quantities such as temperature and density distributions are usually inferred by convolution or using grid cells. Heatload estimation, however, requires the estimation of a surface density on a (usually irregular) 2D geometry embedded in 3D space. This challenge is usually tackled by one of the two following approaches: • The first approach is to directly use the geometry mesh (which is also used during field line tracing) for heat load inversion, essentially by binning particles over the elements of this mesh. The accuracy of this process is tightly linked to the surface resolution of this mesh. If the mesh is too coarse, small-scale hot spots will be missed during evaluation, while a too fine mesh will result in overly large noise from Monte Carlo sampling. While dynamic subdivision techniques exist to remedy this problem, the irregularity of triangle meshes makes this a challenging problem to address. • The second approach is to selectively map parts of the plasma-facing components onto a 2D grid. Once such a mapping is constructed, it is straight-forward to convert particle samples into 2D distributions. Since the target region of the mapping is a regularly spaced rectangular region, changing the sampling density (even with local adaptivity) becomes a relatively simple process. The main drawback of this approach is the large amount of manual work involved, because such transformations are usually manually specified, and are cumbersome to derive for complex shapes which arise outside of the divertor target elements.
To retain the advantages of the second approach while removing most of the manual setup work involved, we decided to automate the transform from 3D to 2D space following the ideas of a synthetic camera. Since an imaging transform is any way necessary for visualization, this removes any duplicate effort related to heat flux inversion. For this purpose, we make the following assumptions about our synthetic camera formalism: • The camera transform is given as T , where the numerator calculates the 2D position relative to the camera origin and viewpoint and the denominator calculates the depth.
• The camera viewport is given by the interval [0, M] × [0, N] ⊂ R 2 and is subdivided into M × N pixels.
The resulting density estimation method consists of the following steps: (a) Apply the perspective transform T to all test particle strike points, transforming them from 3D space into the synthetic camera's image plane. (b) Remove all invisible strike point samples: 1. Remove all occluded strike points based on a maximum per-pixel depth ⃗ c ·⃗ x + d, 2. Remove all strike points outside the camera viewport. (c) Count the number of strike points inside each pixels and divide by the total number of strike points (including hidden / invisible ones). As the pixels are of size 1 × 1, this gives the viewport density of strike points. (d) Multiply the strike point density by a per-pixel correction factor to obtain the heat load density (assuming an even energy distribution over test particles).
In order to carry out this procedure, one needs to obtain two pixel-related pieces of information ahead of time: • The depth of the closest geometry piece along the sight line related to the pixel, • The correction factor required to calculate the heat load density from the pixel-space sample density.
The first quantity can be obtained using well-known rasterization procedures. The second quantity can be obtained using information derived from the transform T. The value one needs to compute here (or at least approximate to first order) inverse of the pixel's area if it were to be back-projected onto the visible geometry mesh. This is challenging to compute, as the transform from 3D to 2D space is, by itself, degenerate. We therefore reuse the information provided during triangle rasterization to compute the required correction factor.
During rasterization, the following pieces of information (besides the imaging transform T) are available to us: • An affine mapping H from the reference triangle T R = {x, y ∈ [0, 1] |x + y ⩽ 1} onto the mesh triangle currently undergoing rasterization, given by the triangle's three spanning points h 0 , h 1 , The current pixel (n, m) ⊂ N 2 viewing the rasterized triangle, • A point x R ∈ T R on the reference triangle such that it gets mapped onto the current pixel T (H (x R )) = n + 1 2 , m + 1 2 , Based on this information, given an infinitesimal rectangular area dA R around x R , we can compute: • the corresponding area on the mesh dA M is the Jacobian of the combined transform from the reference triangle to the camera viewport around x R (Note: H is affine, T is fully nonlinear).
Based on this information, we can estimate (to first order) the inverse back-projected area of the pixel as the ratio of these two: This calculation works with arbitrary camera transforms, captures all transform effects (including zoom, distortion, distance, etc) and only requires forward pass information about the camera transform. Furthermore, it turns out that the inverse pixel can be stably computed, whereas the pixel area could potentially become extremely large (due to a surface being extremely close or the camera looking at it at a steep angle). The generality of this approach has already proven to be highly useful in heat load studies related to equilibrium effects on the magnetic field topology, where it was possible to quickly investigate the complex geometries outside of the divertor plate. This investigation revealed previously unnoticed areas, which were at risk of facing potential overloads in finite-beta configurations [4].

Comparison between models in W7-X reference configurations
In order to compare the validity of both models, heat load distributions are presented exemplarily for two magnetic reference configurations (the 'standard' and 'low-iota' configurations) of the Wendelstein 7-X Stellarator in finite-beta fields (see [4] for a discussion of the calculation procedure and the magnetic configurations), and then compared to infrared camera observations in experimental discharges. As mentioned in section 2.3, we used 10 5 Monte Carlo samples or the simulations (which fixes the attainable resolution of the simulated heat load densities). To eliminate the influence of drift effects (which are not simulated by either model) onto the IR observations, IR camera images from an upper and lower divertor in a single module (module 1) were combined. Please see [1,2] for a more thorough discussion of the Wendelstein 7-X island divertor and [10] for an overview of the magnetic reference configurations.
The heat-flux observations are obtained by running the 2D diffusion code THEODOR [11] on individual fingers of the divertor tiles, constraining the plasma-facing surfaces temperatures to values measured by the divertor infrared cameras [12]. Please see [13,14] for a more detailed description of this process. The heat fluxes are then normalized to P div = P heating − P rad , where P heating is the total power injected into the device (in our case ECRH heating power), while P rad is the radiated power estimated by the bolometer diagnostic. Possible dependencies of these derived heat fluxes include the surface emissivity of the divertor (which is not only material-, but also roughness-dependent), the non-uniformity of the IR camera, reflections of infrared light emitted by other components, as well as the alignment of the cameras themselves. All of the mentioned dependencies, with the exception of possible reflections, are included in the calibration of the heat-load computation process. However, possible deviations beyond the calibration still represent potential error sources, which are difficult to quantify.
• In the magnetic standard configuration (figure 1), both models correctly predict the locations of both the main strike lines on the horizontal and vertical targets (lower and upper components in the pictures) in the left part of the toroidally extended divertor. However, infrared measurements show an additional secondary strike pattern on the right hand side of the horizontal plate, which can not be reproduced by the field line diffusion model without increasing the perpendicular diffusion to levels at which the main strike line is distributed over nearly the entire divertor. This difference indicates that this secondary pattern is related to a counter streaming flow distribution discussed above. • In the low-iota configuration (figure 2) at finite beta, both models show little deviation in the qualitative structure of the divertor heat loads. Like above, there exists an additional heat flux pattern on the right hand side of the horizontal target, as well as an extension of the main strike line into the divertor's middle section. However, these effects appear to be less prominent than in the magnetic standard configuration discussed above. Furthermore, there are disagreements concerning the toroidal distribution of heat loads inside the wetted regions. The anisotropic diffusion model appears to slightly better capture the heat loads onto the secondary spot-like structure adjacent to the main strike line on the horizontal target. On the other hand, it also predicts a concentration of heat flux near this spot inside the main strike line. This contradicts experimental observation, which show the heat load to be concentrated on the center of the strike line. The original field line diffusion model appears to capture the heat load distribution inside the wetted spots better for this configuration. It should be mentioned that this configuration is subject to substantial up-down heat load asymmetries [15], which were eliminated by averaging between observations from an upper and a lower divertor infrared camera.
It can be seen from these observations that the relative behavior of both models and the experiment is highly configuration-dependent. In particular, the characteristics of the main strike line are well captured by classical field-line diffusion. However, the anisotropic diffusion model appears to achieve a closer match to the experimental observations concerning the structure of secondary heat loads (specifically, which areas appear to receive heat loads and which do not). For a design study encompassing many parameters (such as divertor shape parameters), the suggested approach would be to first use field-line diffusion to prune the parameter space, and then apply more expensive models like the anisotropic diffusion model on a reduced set of cases to look for secondary heat load patterns. We would like to stress that we do not believe the anisotropic diffusion model to be strictly more accurate. Both models are proxy models that skip over part of the underlying physics. Therefore, final checks should be performed with a high-fidelity physics-based code where possible.

Conclusions & outlook
Based on the comparison in section 3, it can be concluded that the improved ability to capture counter-streaming flows represents a relevant edge of the anisotropic diffusion model over its simpler and faster field line diffusion counterpart when modeling divertor heat loads, particularly when trying to identify location of potential overload events. Unless the existence of counter-streaming flows can be ruled out based on prior knowledge, care should be taken to include them in the physical modeling process. This work did not investigate alternative approaches to embed such effects into the simpler field-line diffusion model. One such approach, based on repeatedly tracing in alternating directions (by reflecting particles at their impact points) a fixed number of times are currently being tested on Wendelstein 7-X. It could also be concluded that in a direct comparison, the anisotropic diffusion model is not necessarily more accurate, particularly when examining the heat distribution inside wetted regions.
It should also be mentioned that the proxy models discussed here completely neglect the effects of electric fields inside the plasma. In particular, no sheath models or ambipolarity conditions were included for the plasma edge. The large electrostatic potential differences across such an electrostatic sheath can potentially change the flow structure of electrons and ions in the plasma edge. However, to implement such a model for parameter studies, one would have to compromise between accuracy of the (nonlinear) sheath simulation and computation time constraints.
The flexibility of the synthetic camera-based heat load interpretation was found to be of substantial practical value during the modeling studies. Currently, the synthetic camera approach is implemented as a python code accelerated by the NUMBA compiler [16]. It should be possible to implement it directly as a GPU pixel shader. Such an implementation could potentially offer massive speedups over a CPU implementation, as graphics cards are specifically tailored to accelerate this type of computation.
Currently, the relative comparison of the presented models was only performed on a qualitative basis. In the future, it would be desirable to also evaluate the models quantitatively. Such a comparison metric would have to take into account the design constraints of divertor scenarios (different heat load limits on different subcomponents, position deviation vs maximum local heat load deviation), but could potentially provide a valuable proving ground for models to compete in.

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