Synthesizing impurity clustering in the edge plasma of tokamaks using neural networks Open Access

The pursuit of a sustainable and abundant energy source has sparked greater attention toward nuclear fusion, which is the process responsible for powering celestial bodies such as stars. Nuclear fusion is distinct from nuclear fission as it involves the fusion of light atomic nuclei to create heavier ones, resulting in the release of a significant amount of energy. Fusion is highly regarded as the ultimate goal in energy production due to its capacity to generate significant power while causing minimal harm to the environment.

Fusion reactors confine and heat a plasma, consisting of ions and electrons, to enable fusion reactions at high temperatures. However, the task of confining plasma poses a considerable challenge due to the presence of multiple instabilities that can result in energy loss. Impurity accumulation in the plasma core can result in heat loss through radiation, leading to a decrease in confinement quality. Therefore, it is crucial to conduct a thorough investigation of impurity in fusion plasma.

The dynamics of magnetically confined plasma flow is influenced by drift-wave turbulence and zonal flows. The numerical high-resolution simulations for modeling the plasma flow carried out within this work are based on the Hasegawa–Wakatani model, which governs cross field transport by electrostatic drift waves in magnetically confined plasmas. A modified version of the model¹ is likewise investigated to take into account zonal flows. The electric field perpendicular to magnetic field lines is particularly significant because it strongly drives cross field fluxes, impacting edge pressure profiles and stability.^2,3 $E × B$ flows strongly influence the motion of coherent structures, which can account for substantial particle losses in tokamaks.^4–6 The resulting fluxes on walls can cause sputtering, erosion, and impurity injection, further degrading confinement and safe operation.^5,7,8

For the particles, there has been growing interest in Lagrangian perspective in recent years.^9–11 This approach involves the analysis of trajectories of numerous tracers. Numerical simulations involve solving equations to determine the trajectories of test particles within a specific velocity field, such as the $E × B$ velocity field. The Lagrangian approach highlights the significant impact of coherent structures on transport. The combined effect of eddy trapping and zonal shear flows often leads to non-diffusive transport.^12–14 In fusion plasma research, previous studies have primarily focused on passive flow tracers without considering inertial effects in the edge plasma using the Hasegawa–Wakatani model.^15–17 We increase the complexity of these models and add the effect of inertia, which has not been done so far, to the best of our knowledge. Therewith, we can study the behavior of heavy and light atoms and see the impact on the impurity distribution. A crucial aspect of understanding inertial impurities is their self-organization, such as clustering and void formation, which can be quantified mathematically by deviation from Poissonian statistics.¹⁸ Their dynamics can then be analyzed using tesselation based techniques.¹⁹ 

There is a growing interest in studying spatial patterns in turbulence exploring machine learning.^20,21 Our study uses machine learning to estimate mesoscale inertial particle clustering in turbulence, relying on flow field data rather than individual particle information. We apply machine learning tools that have been successfully implemented in the field of hydrodynamics, see Refs. 22 and 23. Three neural networks—Autoencoder, U-Net, and Generative Adversarial Network (GAN)—are trained to predict impurity particle densities from vorticity field snapshots. The results are visually and statistically compared with actual DNS particle data. This approach reduces the cost for predicting particle density distributions.

This paper is structured as follows: Sec. II outlines the theoretical basis used for simulation, introducing the Hasegawa–Wakatani edge plasma model and the motion of inertial impurity particles. Section III details the numerical simulations and results. In Sec. IV, we explore machine learning models such as Autoencoders, U-net, and GAN to synthesize impurity concentration fields. Finally, Sec. V summarizes our findings and suggests potential directions for future investigations.

The numerical simulations are based on Hasegawa–Wakatani (HW) equations for plasma edge turbulence driven by drift-wave instability.²⁴ Our focus is on the two-dimensional, slab geometry of the HW model, see, e.g., Ref. 10. Figure 1 depicts a representation of the flow configuration. The magnetic field lines are assumed to be constant, straight, and perpendicular to the slab. Ions are considered cold, and the effects of the temperature gradient are ignored.

In the context of the classical Hasegawa–Wakatani model (cHW), as detailed earlier, zonal flows are not observed. To specifically obtain zonal flows, a modified Hasegawa–Wakatani model (mHW) can be considered. This modification was introduced in Ref. 1, where the coupling term $c ( ϕ − n )$ is set to zero for modes k_y = 0.

Previous studies^15–17 have assumed that impurity particles perfectly follow the fluid flow without deviating from the flow path. However, this assumption may not always be valid, particularly when the particles have significant mass, resulting in observable inertial effects.

The simulation was conducted within a domain with periodic boundary condition that spans an area of 64², and it was discretized using a resolution of 1024² grid points. The time step was $5 × 10 − 4$⁠, while values for diffusivity of plasma density D and kinematic viscosity ν were both set to $5 × 10 − 3$⁠. The value for κ was set to 1. Simulations for the flow begin with Gaussian random initial conditions and run until achieving a saturated, fully developed turbulent flow (see Ref. 10). After this phase, 1 × 10⁶ uniformly distributed inertial impurity particles are injected.

The systems of HW equations, specifically Eqs. (1) and (2), are solved using a Fourier pseudospectral method. It is particularly well-suited for problems with periodic boundary conditions and provides excellent accuracy for smooth solutions.³⁰ The equations governing the motion of the impurity particles, namely, Eqs. (5) and (6), are solved using a second-order Runge–Kutta (RK2) scheme. This method provides a balance between accuracy and computational efficiency. Simulating a large ensemble of particles presents a considerable computational challenge. Our solution leverages High-Performance Computing (HPC) to overcome this complexity, employing parallel computing with Message Passing Interface (MPI). For further details on the numerical method and its implementation, we refer to Ref. 10.

In this paper, we focus on the quasi-adiabatic regime (c = 0.7), which is the relevant case for edge plasma of tokamaks. For readers interested in zonal flow, relevant discussions can be found in  Appendix A. Figure 4 shows the visualization of vorticity fields and 10⁴ superimposed impurity particles for various Stokes numbers in the quasi-adiabatic regime. From Fig. 4, we can observe clustering of impurity particles for the case of $S t = 0.25$ and St = 1, while for St = 0 and St = 50, the impurity density field homogeneously fills the entire physical domain and does not show significant spatial correlation. When the particles are very light, the response time of the particles is much smaller than the turnover time of the turbulent eddies, i.e., $τ p ≪ τ η$ (⁠ $S t ∼ 0$⁠). Their response time is so quick that they can adjust their velocity to match any changes in the flow immediately. They behave as usual tracers like those in previous studies.^15–17 When $τ p ≫ τ η$ (⁠ $S t ≫ 1$⁠), the inertia of the impurity particles is so high that they are almost unaffected by the eddies. Particles are moving ballistically and are randomly distributed, resulting in less clustering and more dispersion. When the response time of the particles is comparable to the eddy turnover time, i.e., $τ p ∼ τ η$ (⁠ $S t ∼ 1$⁠), particles are subject to centrifugal effect of turbulent eddies. Eddies can accelerate the particles to a degree where they are thrown out of the eddy due to their inertia. This ejection and the subsequent movement toward areas of lower vorticity result in particle clustering. The clustering of impurity particles results in inhomogeneous poloidal distributions (y direction in the domain) of the impurity density, which can modify radial collisional transport.³² 

To analyze and statistically interpret the distribution of impurity particles, we increase the number of particles from 10⁴ to 10⁶. Displaying the position of each particle in the flow is impractical due to this large quantity. Instead, we visualize the particle density. We achieve this by dividing the entire domain into 512 × 512 boxes and defining the impurity density n_imp as the number of particles in each box, i.e., we compute a histogram. Figure 5 shows the impurity density field for the case of c = 0.7 considering different Stokes numbers. When St = 0, the impurity density is distributed uniformly throughout the flow. As St increases, variations in impurity density become apparent, with some areas exhibiting high density and others low. This uneven distribution indicates the formation of voids and clustering of particles. When St = 50, the impurity density returns to a more uniform distribution. The impurity density fields exhibit patterns similar to those seen in the earlier visualization involving 10⁴ impurity particles.

To analyze the impurity density distribution at different scales, we compute the Fourier energy spectra similar to Eq. (7), replacing the flow velocity u by the impurity density distribution n_imp. The spectral characteristics of impurity density field are shown in Fig. 6 in terms of the Fourier energy spectra. The variance of the impurity density is equivalent to the cumulative sum of the corresponding spectral energy. This is visually represented as the area under the energy spectra curve. As observed, the variance initially increases with the Stokes number up to St = 1, indicating stronger fluctuations. From St = 1 to St = 50, there is a decrease, suggesting less fluctuations. For St = 0, representing tracer particles without inertia as in previous studies, the impurity particles are distributed randomly. Random spatial distribution of particles follows a Poisson probability density function (PDF) with mean and variance equal to the mean density $⟨ n imp ⟩$⁠, see, e.g., Ref. 31. The slope of the spectra for St = 0 in 2D is 1, as prescribed by the equidistribution of variance among all wavenumbers. The spectrum for St = 50 is observed to collapse on those of St = 0 at small scales (high wavenumbers), thus resembling the same characteristics as for randomly distributed particles.

Simulating flow without impurity particles to obtain the vorticity field is not computationally expensive, but including and tracking 1 × 10⁶ impurity particles increases costs. Since a specific statistically stationary distribution of impurity density is linked to a specific vorticity field, we propose creating a neural network that can directly predict the distribution of impurity density based on the vorticity field. This network uses as input the vorticity field (easily obtained from DNS) and produces as output the corresponding impurity density field for a specific Stokes number. With this neural network, we can run DNS to get the vorticity field, use it as input for the network, and efficiently predict the impurity density distribution for a specific Stokes number. This approach eliminates the need for costly simulations involving a million particles.

We will evaluate and compare three different neural network architectures that have previously been utilized in particle-laden hydrodynamic turbulence studies:^22,23 the Autoencoder, U-Net, and Generative Adversarial Network (GAN). To determine the most effective neural network architecture for accurately predicting the impurity density fields, a thorough analysis of the statistical attributes of the impurity density fields generated by these neural networks is carried out. In the following, we will consider the cHW model for the quasi-adiabatic case (c = 0.7).

The models used in this study are based on those utilized in Refs. 22 and 23. Figure 7 displays the architectures tailored from the standard structure, which we have implemented using TensorFlow.

• Autoencoder: Autoencoder³³ is the basic model in our study, comprising of an encoder that compresses the data, and a decoder that reconstructs the output from this compressed version.

• U-Net: U-Net³⁴ is essentially an Autoencoder but with added skip connections. These connections facilitate non-sequential connections between layers, which helps in the preservation of information at different scales throughout the network.

• GAN: GAN (Generative Adversarial Network) consists of two neural networks,³⁵ a Generator and a Discriminator, trained simultaneously through adversarial processes. The Generator attempts to produce synthetic data, while the Discriminator tries to distinguish between real and synthetic data. We use U-Net as the generator.³⁶ More details can be found in Refs. 22 and 23.

In this study, we use the data of one long simulation in the statistically stationary regime. This simulation has a temporal span equivalent to 220 eddy turnover times. Within this simulation, we have captured 400 snapshots of the vorticity field. For each snapshot of the vorticity field, there are five corresponding impurity density fields (⁠ $S t = 0.01 , 0.05 , 0.25 , 1 , and 5$⁠). The neural network is then trained to learn the relation between the vorticity field and the corresponding impurity density field. 70% of the dataset is for training the model and 30% for evaluating its performance. Vorticity data are normalized between −1 and 1. Impurity density data are normalized between 0 and 1 to enable the use of binary cross-entropy.

The Adam optimizer³⁷ was used to train the neural networks. Training on a single sample from the training set is termed a “step.” The Autoencoder and U-net models were trained for 28 000 training steps. The GAN model required more steps for convergence and was trained for 56 000 steps. The Python codes were run on an NVIDIA graphics card Tesla V100 32GB.

Figure 8 shows the vorticity field and corresponding impurity density field (⁠ $S t = 1$⁠) obtained from DNS for the case of c = 0.7 (cHW). Figure 9 illustrates impurity density field (St = 1) predicted from Autoencoder, Unet, and GAN for c = 0.7 (cHW) using vorticity field as input. When comparing DNS impurity density fields with those predicted by the Autoencoder and U-Net, we find that both DNS and predicted fields exhibit similar structures, but the clustering in the predicted fields is slightly fuzzy. The GAN model successfully predicts void areas and filament-like structures to some degree of success, but it sometimes struggles with predicting the details of regions with impurities.

To better understand the performance of our models, we analyze the results using the probability density functions (PDFs) and energy spectra of the impurity number density. Figure 10 shows the PDFs of DNS data and predicted data using three different architectures for St = 1. From Fig. 10, we observe that the GAN model excels in predicting the density distribution, as its predictions closely align with DNS. In contrast, U-Net and Autoencoder exhibit shortcomings, particularly in predicting higher densities. The deviations in the PDFs of U-Net and the Autoencoder suggest that these models may struggle to accurately replicate the intricate characteristics present in high-density regions.

Figure 11 shows the energy spectra of DNS data and predicted data using three different architectures for St = 1. In Fig. 11, we observe that the GAN's energy spectra are nearly identical to that of the DNS data. This demonstrates the remarkable accuracy of the GAN model in predicting the energy spectra of the impurity density. On the contrary, both U-Net and the Autoencoder exhibit a significant drop at high wave numbers, which correspond to fine structures in the data. This indicates that U-Net and the Autoencoder are unable to accurately capture these fine-scale structures and may struggle to represent the intricate details present in the impurity density.

Overall, the analysis of the PDFs and energy spectra reaffirms that GAN outperforms the other two for St = 1. Thus, we utilize GAN in the following to investigate more Stokes numbers values. From the PDF of impurity density as shown in Fig. 12, we can see that for $S t = 0.25$⁠, GAN predicts lower impurity density when compared with DNS. For small St (⁠ $S t = 0.01 , 0.05$⁠), GAN predicts more particles than DNS. When dealing with a large Stokes number (St = 5), the PDF is similar to the DNS results. In the case of St = 1, the GAN predictions align perfectly with DNS. For energy spectra, as shown in Fig. 13, it is noticeable that the GAN predictions align very well with the DNS. In the GAN predicted energy spectra, notable spikes were identified at frequencies corresponding to the edges and corners of the domain. This could be because the initial encoding of domain periodicity was inadequate. To correct this, adjustments in the code are needed for true periodic convolution. Unfortunately, as of the time of this writing, TensorFlow does not yet provide native support for periodic padding.

In our study, we train the AI model specifically for c = 0.7, which is a typical value for edge plasma. This model predicts impurity density fields for $S t = 0.01 , 0.05 , 0.25 , 1$⁠, and 5 using the vorticity field at c = 0.7 as input. The model's validity is specific to c = 0.7 and the aforementioned Stokes numbers. Its extrapolability is shown and discussed in  Appendix B.

Now, we know that GAN is the optimal model. Next, we can conduct simulations to generate more vorticity fields without particles, a process which is computationally efficient. We then input these vorticity fields into GAN to predict the impurity density fields. Subsequently, we can perform detailed analyses of transport characteristics and clustering, using methods like multiresolution analysis and Voronoi tessellation. There are several possible ways to enhance the model. One method is to adjust architectural features like the number of layers, filter sizes, and channel counts. Another effective approach is to integrate physical constraints into the model, ensuring that predictions align with established physical laws, such as mass or energy conservation, for more accurate and plausible results.

This study focuses on inertial impurity within fusion plasmas, which is different from previous literature where the focus has typically been on “passive tracers” following the flow. We used high-resolution numerical simulations based on Hasegawa–Wakatani equations with impurity particles having different inertia. We studied the quasi-adiabatic regime (c = 0.7), typical for Tokamak edge plasma. Our findings show that the inertia of particles, quantified by the Stokes number St, plays a vital role in their behavior. At lower St values, impurities move along the flow streamlines. As St increases, impurities tend to cluster at low-vorticity regions. At very high St values, impurities show random movement due to substantial inertia, resulting in much less clustering.

As simulating 1 × 10⁶ or more impurity particles is a computationally intensive task, we utilized machine learning techniques to create a surrogate model. Three neural networks—Autoencoder, U-Net, and GAN—were developed to synthesize impurity densities, using vorticity fields as input. All models produced visually comparable results to simulations, with varying degrees of success. The GAN closely matches DNS values in terms of probability density functions and energy spectra of density. This confirms GAN's effectiveness in modeling both the general distribution and finer structures. Based on the results, future work will involve: improving GAN by adjusting the architectures or involving physical constraints directly into the model; running simulations for more vorticity fields, and then using the GAN to predict impurity fields for analyzing clustering and transport properties; including electromagnetic interactions for more accurate impurity particle models.

Z.L., T.M.O., and K.S. acknowledge the financial support from I2M and the French Federation for Magnetic Fusion Studies (FR-FCM) and the Eurofusion consortium, funded by the Euratom Research and Training Programme under Grant Agreement No. 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. Z.L., T.M.O., B.K., P.K., and K.S. acknowledge partial funding from the Agence Nationale de la Recherche (ANR), project CM2E, Grant No. ANR-20-CE46-0010-01. Centre de Calcul Intensif d'Aix-Marseille is acknowledged for providing access to its high-performance computing resources.

The authors have no conflicts to disclose.

Z. Lin: Writing – original draft (equal); Writing – review & editing (equal). T. Maurel–Oujia: Writing – original draft (equal); Writing – review & editing (equal). B. Kadoch: Writing – original draft (equal); Writing – review & editing (equal). P. Krah: Writing – original draft (equal); Writing – review & editing (equal). N. Saura: Writing – original draft (equal); Writing – review & editing (equal). S. Benkadda: Writing – original draft (equal); Writing – review & editing (equal). K. Schneider: Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

In the modified Hasegawa–Wakatani model at c = 2 with 10⁴ impurity particles as shown in Fig. 14, the zonal flow is very strong, and the vorticity values are weaker compared to cHW. The particles appear to be distributed evenly across the fluid, and in proximity to the zonal flow, it seems that clustering is not evident. To gain deeper insights, we examined the impurity density fields obtained from 1 × 10⁶ impurity particles. As depicted in Fig. 15, there is some preferential concentration of impurities in the shear layer region for St = 0.25, 1, and 5. The impurity particles are ejected by the zonal flow, and for St = 1, the effect seems to be the strongest. For St = 0 and 50, the impurity particles are randomly distributed.

In the following, we test the extrapolability of the GAN model for different c-values. The model has been trained using the data with c = 0.7. As input data, we then use vorticity fields for c-values different from 0.7 and assess the impurity densities obtained with the GAN model. The results are shown from Figs. 16–19. In the hydrodynamic regime (c = 0.01), we observe that GAN predicts the PDF and energy spectra of impurity density very well at $S t = 1$⁠. At other Stokes numbers, discrepancies arise between GAN predictions and DNS results in terms of PDFs and energy spectra. In the adiabatic regime (c = 2), the PDFs and energy spectra of impurity density predicted by GAN tend to have larger values than the DNS results. Note that the results are good from the statistical point of view, but the obtained impurity density fields differ visually from the DNS results.