Gyrokinetic flux-driven simulations in mixed TEM/ITG regime using a delta-f PIC scheme with evolving background

Magnetic fusion research relies heavily on its accurate modeling by computer simulations. In the most promising reactor configuration, the tokamak, the plasma is confined by magnetic fields in a toroidal vacuum chamber. A complete description of the plasma involves simulating regions of the core, edge, the scrape-off-layer (SOL), and plasma–wall interaction.

To simulate fusion plasmas, many methods exist, which can be categorized by the physical assumptions made, dictated by the physical process of the plasma volume considered. This work focuses on the gyrokinetic particle-in-cell (PIC) method.^1–5 The gyrokinetic formalism^6–8 reduces the number of phase space variables from six to five, approximating the dynamics of plasma particle trajectories by gyrorings bound to evolving gyrocenters. The reduction in dynamics implies a timescale separation between the fast cyclotron motion and the typical fluctuation time scales involved in turbulent processes. The PIC scheme relies on representing the evolution of particle distributions in phase space in terms of Lagrangian “markers,” whose characteristics are dictated by the equations of motion. In this approach, the sources of the fields, i.e., the charge and current densities, are obtained via Monte Carlo integration, thus giving the distribution a statistical interpretation.⁹ This also means that it inherits the main problems of Monte Carlo sampling, notably the statistical sampling error problem referred to as “noise” in the following. Unless appropriate noise control measures are implemented, this problem severely limits the physically relevant simulation time. This is especially true for cases that exhibit significant deviation from initial profiles and/or large relative fluctuation amplitudes.

The aim of this work is to improve the statistical sampling noise problem of the PIC scheme to apply it to long transport timescale simulations, i.e., targeting cases with significant deviations from initial profiles. For plasma simulations where the distribution function f of each species (ions, electrons, possible impurities) does not deviate more than a few percent from its initial state $finit$⁠, one usually adopts the delta-f splitting.^4,10 The approach decomposes f into a stationary (and often analytic) distribution f₀ and a time-dependent perturbation part δf, where only the latter is represented by numerical markers. This so-called delta-f PIC method is to be contrasted with the full-f PIC scheme, which represents the whole f in terms of markers. The gain in noise reduction of the delta-f scheme relies on the reduced variance of the marker weights, provided that the assumption $||δf||/||f||≪1$ for some definition of the norm $||·||$ is met. However, for plasma simulations exhibiting a large perturbed component δf, i.e., such that $||δf||/||f||∼1$⁠, one usually falls back to the full-f scheme, which entails using high marker numbers to achieve low enough noise levels. To still possibly retain the advantage of the delta-f scheme, one can also evolve^11–13  f₀, albeit at a slower timescale than that of the fluctuating δf. This work explores the benefits of having such a time-evolving background by constraining f₀ to be a flux-surface-dependent Maxwellian which is time-dependent via its evolving gyrocenter density and temperature profiles. Another source of statistical sampling noise is related to “weight-spreading”^12,14 resulting from the implementation of collision operators in the delta-f PIC scheme using a Langevin approach. However, this problem will not be addressed in this work as collisions are not considered.

Following the success of a previous work¹⁵ by the author using a similar approach in a simplified setup (sheared slab geometry, adiabatic electrons), the adaptive scheme is implemented in tokamak relevant axisymmetric toroidal geometry in the framework of the global gyrokinetic PIC code ORB5.¹⁶ The simulations are electrostatic and “flux-driven,” with gyrokinetic ions while electrons have a hybrid response¹⁷ and instabilities being driven by both ions and electrons in a mixed Ion-Temperature-Gradient-Trapped-Electron-Mode (ITG-TEM) regime. Background density and temperature adaptations are made for both species independently. When compared to nonadaptive cases, results from the adaptive scheme exhibit lower errors resulting from statistical sampling noise. Evolving the plasma profiles in the presence of radially localized sources up to their quasi-steady state using the adaptive scheme will be shown to require much lower marker numbers than with the nonadaptive scheme.

This paper is organized as follows. Section II introduces the physical model to be solved, namely, the Vlasov–Maxwell system in the first-order gyrokinetic approximation. A discussion of the Quasi-Neutrality-Equation (QNE) then follows, elaborating on the electrons' hybrid response. It concludes by describing the functional form of the different possible source terms used in this work. Section III elaborates on the use of a control variate in the form of either a canonical or local Maxwellian, with flux-surface-averaged (f.s.a.) density, parallel flow and temperature having an explicit time-dependence. The relaxation equations that connect the three lowest-order velocity moments of f₀ and δf are then explained. The modification of δf requires a correction term to be added to the right-hand side (r.h.s.) of the QNE. Section IV details the considered initial profiles of density and temperature, along with the form of the actual heat sources for each species used in this work. Information about grid resolution, Fourier filtering, source strengths, and adaptive scheme parameters are given. Section V A begins the result section with a discussion on time traces of various transport parameters and sampling noise diagnostics, comparing them between nonadaptive and adaptive cases. Section V B investigates the effect of the adaptive scheme on time-averaged profiles of f.s.a. density and temperature. This section also examines the need of the QNE r.h.s. correction. Section V C diagnoses the time evolution of sampled phase-space volume and marker distribution, which reveals in some simulation cases a problem of under-sampling. Section V D discusses the f.s.a. profiles of density and temperature, time-averaged over the quasi-steady state, results of which could be acquired only under the adaptive scheme. Section VI then summarizes the main merits of the adaptive scheme demonstrated in this work, along with possible future developments and applications directions.

The grid resolution for the radial s, poloidal $θ⋆$ and toroidal $φ$ is taken to be $Ns×Nθ⋆×Nφ=256×512×256$⁠, where N represents the number of intervals. Toroidal Fourier modes in the range of $0≤n≤64$ will be simulated, only keeping poloidal Fourier modes $m∈[nq(s)−Δm,nq(s)+Δm]$ with half-width $Δm=5$⁠. $(kθρs)max∼1$ and $(krρs)max∼3$⁠, respectively. To lighten the computational cost “heavy” electrons are considered, i.e., with an electron-ion mass ratio $me/mi=1/200$⁠. The time resolution is $csΔt/a=4.1×10−3$⁠, and the maximum linear growth rate is found to be $γmaxa/cs=0.931$⁠. The strength of the Krook operator $Sn$ for noise control is set to $γn=11%γmax$⁠. For the flux-driven simulations considered here, the radial heat source profile $GH(ψ)$ of Eq. (7) is approximated by fitting a Gaussian function around the peak of the time-averaged effective heat source at quasi-steady state of a previously run temperature-gradient-driven simulation with the initial profiles given by Table I. The heat source strength $γH$ is determined by equating the integrated power to that of the effective heat source in the radial range $s∈[0,sb]$⁠. The heat sink at the radial edge is replaced by the Krook buffer S_b. The calculation for the heat source is done for both ions and electrons separately. An example of the heat source profiles for ions and electrons are shown in Fig. 2.

Simulations are initialized by loading markers in phase space and initializing their weights using quasi-random distributions. All adaptive cases have adaptation rates $αn=αE=11%γmax$ and adaptation period set with N_t = 100 for both ions and electrons. Even when considering the hybrid electron model, adapted electron profiles include contributions from both passing and trapped electrons. Unless stated otherwise, all time traces have a moving time averaging window of $csΔt/a=4.085$⁠. Marker numbers N_p are displayed in millions (M). Comparison will be made between the nonadaptive (standard) case with N_p = 256M, 128M and the adaptive case with N_p = 128M, 64M. The choice of simulations with different N_p allows for the discerning of sampling noise's effect on different diagnostic measures.

We see that the shearing rate for the standard cases are consistently higher than that of the adaptive cases. As $ωE×B$ amplitudes measure the rate of turbulent eddy shearing, higher $ωE×B$ amplitudes in the standard cases suppress turbulence more efficiently, and thus, allows for steeper gradients to be maintained for a given flux level, than in the adaptive cases within $s∈[0,7,0.9]$⁠. (Notwithstanding the fact that in this work, no standard cases reached quasi-steady state, we assume that both standard and adaptive cases would share the same quasi-steady state, and thus, the same average heat flux.) This explains the consistent higher χ_i of the adaptive cases as compared to the standard cases. As the system finally reaches quasi-steady state for $cst/a≳1000$⁠, simulations under both standard and adaptive schemes show signs of asymptoting to the same χ_i and $|ωE×B|$ values. This asymptoting behavior of the standard case is not conclusive as simulations with $Np>256$ M need to be run for an even longer simulation time to reach quasi-steady state for comparison with the adaptive cases presented here. Due to the numerical cost of such simulations, these standard runs are not conducted during the work of this paper.

To see how electron diffusivity χ_e compares to χ_i, Fig. 3(b) shows the ratio $χi/χe$ for all cases considered. In the time interval $cst/a∈[0,150]$⁠, the system evolves gradually from an initially TEM-dominated regime to a mixed TEM-ITG regime, with a $χi/χe$ ratio evolving from about 0.2 to 1.2. For further times, this ratio remains at the same value. This change of regime is due essentially to the density profile evolution, as will be discussed further (see Fig. 11).

As the system approaches quasi-steady state, with the buffer being the only source/sink of particles, the particle flux is expected to reduce and reach near zero values. Note that in ORB5 the markers that leave the domain are re-injected with their weights set to zero. While this represents de facto a source term in the plasma volume, for the simulations of this paper the buffer is strong enough to have reduced marker weights before they reach the boundary. Figure 5 shows the f.s.a. ion gyrocenter flux $〈Γ→·∇ψ〉/〈|∇ψ|〉$ for cases with various N_p under standard and adaptive schemes. We first note that adaptive cases (blue and green curves) have consistent slightly higher ion gyrocenter flux compared to the standard cases (black and orange) right after the initial overshoot in the time window $cst/a∈[50,200]$⁠. This is again consistent with the lower $ωE×B$ amplitudes of the adaptive cases (see Fig. 4), thereby allowing for higher turbulence levels, and therefore, particle transport to develop. Nonetheless, the standard case with N_p = 256M (black line) seems to initially converge to the same quasi-steady state value as that of the adaptive cases.

The trends displayed in Fig. 7 imply that to reach a similar SNR as the one at the end of the adapted case with N_p = 64M, one would need at least N_p = 512M when using the nonadaptive scheme. Thus, the adaptive scheme allows for a reduction in computational cost by a factor of 8 for ensuring the same numerical quality. For a given N_p, the adaptive cases have their SNR values drop much less from their maximum values compared to the standard cases. This drop happens mostly at the initial phase $cst/a≲200$ of the simulation when profiles evolve the most. This reduction could be lessened somewhat further by increasing the adaptation rates α_n and α_E of Eqs. (12) and (14), though improvements resulting from increasing adaptation rates have been found to be marginal. This is because the adaptive scheme of this work is based on a f.s.a. control variate. Any variation of the fluctuations in the poloidal direction, for example, will not be accounted for in f₀. Nonetheless, for the cases studied, we conclude that a simulation run with N_p = 64M, or even potentially N_p = 32M, gives us significant results, which otherwise, i.e., with standard scheme, would only be obtained with at least N_p = 512M.

Finally, the importance of the correction term given by Eq. (15) to the quasi-neutrality Eq. (5) is illustrated by its effect on the zonal flow shearing rate $ωE×B$⁠. Figure 8 shows the time evolution up to $cst/a=123$ of the radial profile of $ωE×B(s)$ for four simulations run with marker number N_p = 128M. As the plots shown are still early in simulation times, according to Fig. 7, they are considered reliable, whether under the standard or adaptive schemes. Using the standard case in Fig. 8(a) as reference, we first note that strong inward avalanches are triggered from the edge, in a region where the time-averaged $ωE×B$ has positive values. The inward direction of avalanche propagation is in line with the analysis presented in Ref. 32. We also observe a stationary radial corrugation structure coinciding with the q = 1 flux surface at s = 0.538. Small transport barriers tend to develop around low-order Mode-Rational-Surfaces (MRS) including corrugated zonal flow shearing rate profiles. This effect, related to the nonadiabatic passing electron dynamics, is only partially captured by our hybrid electron model, given that it only accounts for the f.s.a. kinetic contribution from passing particles [see Eq. (5)]. Such structures have been previously analyzed¹⁸ using either a Padé or an arbitrary wavelength solver for the ion polarization density, and using fully drift kinetic electrons. The role of the correction term (15) is shown by comparing Figs. 8(b)–8(d) against the reference case in Fig. 8(a). As the control variate is close to a f.s.a. function, one expects persistent f.s.a. features to be gradually erased as low order f.s.a. velocity moments are being removed from the δf component, which appears on the r.h.s. of Eq. (5). This explains that without the correction term, Fig. 8(d), the simulation fails to resolve these persistent corrugated structures. Based on Fig. 8(c), the prescription of the adapted electron density given by Eq. (16) indeed preserves the persistent f.s.a. structures. Though not shown, adaptive runs with QNE correction term calculated by quadratures [Fig. 8(b)] or by means of Eq. (16) are indistinguishable under the simulation parameters used in this work.

We now consider evolved f.s.a. profiles, further averaged over the first time window $cst/a∈[592,613]$ (shown e.g., in Fig. 3). All profiles are derived from binning marker $f0Ω$ values and w-weights into radial bins for the contributions from the background f₀ and perturbed δf distributions, respectively. Figure 9 shows the gyrocenter temperature logarithmic gradient for all cases considered and for both ions and electrons. We first note the dip toward zero in logarithmic gradients for both species inwards from the heat sources (see Fig. 2) at s = 0.35. This dip moves further toward the magnetic axis at later times as the heat is transported from the heat source peak toward the core, thus, flattening the temperature profiles in this region (also see Fig. 10). Figure 9(a) shows a steepening of gradient in $s∈[0.9,1.0]$ (well within the buffer region) for the ions, which is barely present for the electrons [in Fig. 9(b)]. Conversely, in the adaptive cases Fig. 9(b) shows a corrugation of the electron logarithmic gradient in the vicinity of the q = 1 flux surface at $s≈0.55$⁠. This can be further seen in the inserted figure in Fig. 9(b), showing a zoom-in of this region and illustrating the ability of the evolving background to capture such fine profile features. No such corrugation is observed for the ions. This corrugated structure is related to the one of the $ωE×B$ radial profiles of Fig. 8. Finally, we see that for both species, the adaptive cases resulted in a lower logarithmic gradient in region $s∈[0.6,0.9]$⁠. This is related to the fact that adaptive cases exhibit greater local heat diffusivity values [see Fig. 3(a)] due to lower local $ωE×B$ amplitudes, as mentioned.

Ion and electron temperature profiles are shown in Fig. 10. Note that while at this time, there is no off-axis peak ion temperature profile, there is one for the electron temperature under the standard scheme, at around s = 0.45. This does not coincide with the peak of the heat source for the electrons at around s = 0.55 [see Fig. 2(b)]. The slight corrugation for the electron temperature at s = 0.55 was emphasized in its logarithmic gradient in Fig. 9(b). One further notes an increase in ion temperature at the magnetic axis under the standard scheme, with the N_p = 128M case having a larger increase than that of the N_p = 256M simulation. Only a slight increase in ion temperature at the magnetic axis under the adaptive scheme is observed. We thus attribute this problem as numerical, and is related to the under-sampling of phase space volume discussed in Sec. V C. This is already suggested by the dip in background temperature profiles (dashed) of the standard cases (black and orange). As these background profiles are not adaptive, we expect them to coincide with the initial profile (dashed gray). Between the standard cases, the case with N_p = 128M deviates more from its initial values compared to the case with N_p = 256M. This effect is not observed for electrons. This sampling problem of phase-space volume will be discussed in Sec. V C.

Figure 11 shows the gyrocenter densities of the ions and electrons time-averaged over $cst/a∈[592,613]$⁠. Focusing first on Fig. 11(a), we see that at the considered stage of the simulation, the nonadaptive background density of the standard cases (dashed black and orange) has a lower value near the magnetic axis as compared to the initial profile (dashed gray). The source of this deviation is the same as that in Fig. 10(a) for the ion temperature, i.e., under-sampling of phase space. Next, the f₀ contribution to density of the adaptive cases (dashed blue and green) is nearly identical to the one including the δf contribution, which indicates that the time-dependent background density has correctly captured the evolution of the total f.s.a. density. Furthermore, the adaptive cases seem to have converged in marker number N_p, whereas the standard cases has not, with N_p = 256M for the standard case showing a smaller deviation of the evolved density profile compared to the N_p = 128M case. Figure 11(b) shows the difference between ion gyrocenter and electron densities, where we see a larger difference between these densities for the standard cases. Between these standard cases, the case with N_p = 128M shows a larger density difference as compared to that with N_p = 256M, thus, indicating that this density difference is partly due to sampling error, as pointed out in Figs. 10 and 11(a). The same N_p trend can be observed for the adaptive cases, though the results are clearly already more converged. Note that ion gyrocenter and electron density difference is not expected to be zero, as one still needs to account for ion Finite-Larmor-Radius (FLR) effects and ion polarization density. Taking a step back, we now see in Fig. 11(a) that the apparent convergence of the standard case with N_p = 128M (solid orange) with the adaptive cases (solid blue and green) is just a coincidence.

Figure 13 shows the various radial cuts of the ion marker number per $(s,v2)$-bin corresponding to Figs. 12(b) and 12(c), respectively, i.e., at the same instant $cst/a=612$ and with N_p = 128M. We note a general accumulation of markers toward lower energies. Near the magnetic axis at $s=0.02,0.10$ (black, red) however, there is greater marker accumulation to lower energies for the adaptive case compared to the standard case, thus, concluding that the under-sampling near the magnetic axis is the result of phase-space volume corruption. The smaller extent of the under-sampling in Fig. 12(c) for the adaptive case seems to be ensured by the larger accumulation of markers shown in Fig. 13(b). Further investigation in a future work will be required to fully understand this issue.

The problem of under-sampling for the standard cases Figs. 12(a), 12(b), and 13(a), thus, explains the density drop in Fig. 11(a) and temperature rise in Fig. 10(a). The former is due to the lack of markers being accounted for, and the latter is due to the lack of the contribution of low-energy markers. As temperature is the average energy over density, the under-estimated density also contributes to the over-estimation of temperature. In terms of the gyrokinetic equations of motion (2), the difference between the standard and adaptive schemes lies in the evaluation of the E × B drift and parallel acceleration terms involving the self-consistent gyroaveraged electrostatic potential $ϕ$⁠. The adaptive scheme allows for a more accurate evaluation of the r.h.s. of the QNE (5) using a better control variate f₀, this in turn leading to a more accurate solution for $ϕ$⁠. Moreover, under the standard scheme, a better evaluation of $ϕ$ requires a higher N_p. Indeed, based on Figs. 12(a) and 12(b), the under-sampling issue near s = 0 improves with higher N_p. Even though the incompressibility of phase space should be verified for any field $ϕ$⁠, be it self-consistent or not, nonsmooth solutions of $ϕ$ to (5) as a result of poor (noisy) integration of the r.h.s. can lead to E × B drifts and parallel acceleration in the marker trajectories, which are difficult to integrate in time, thus, violating phase space incompressibility. Though not shown, no such under-sampling is observed for the electrons in both the standard and adaptive cases.

This section describes profiles for simulations that have been able to reach quasi-steady state with just one marker loading at initial time. For the considered heat sources, only adaptive cases clearly maintain adequate SNR levels to reach the simulation time $cst/a=1266$ as shown in Fig. 7.

Figure 14 shows the temperature profiles for the ions and electrons separately, contributed by the respective total distributions $f0+δf$⁠. Profiles contributed only by the background f₀ are not shown as they closely follow the total contribution, indicating that the adaptation rates are large enough and as a result $〈δf〉$ is very small. Also shown are the corresponding temperature profiles of the previous time window $cst/a=[592,613]$ (colored dashed). We see that f.s.a. temperature profile at $cst/a=1266$ has further risen mainly around s = 0.4, resulting in this region in a maximum deviation of about 20% from initial values. As already mentioned for Fig. 10(a) the increase in temperature at the magnetic axis for the ions, Fig. 14(a), is due to the under-sampling problem discussed in Sec. V C and illustrated in Figs. 12(c) and 12(d). Except for the temperature rise at the magnetic axis, the profiles with N_p = 128M and N_p = 64M of the respective species seem to be converged.

The logarithmic gradients of these temperature profiles are shown in Fig. 15. The risen peak in temperature around $s∈[0.3,0.4]$ has caused an increase and decrease in the $−〈R/LT〉$ values to the left and right of this region, respectively. Furthermore, the temperature peak has led to negative values of $−〈R/LT〉$⁠. Both ions and electrons also have increased temperature gradients in the pedestal region $s∈[0.8,1.0]$⁠. This is probably caused by edge buffer region and homogeneous Dirichlet condition at s = 1 for $ϕ$ both contribute to reduction of turbulence leading to profile steepening at the outer boundary. Logarithmic gradient values at the flat gradient region $s∈[0.6,0.9]$⁠, conversely, remains approximately constant in time after a decrease before the first time window $cst/a∈[592,613]$⁠. The corrugation of the electron temperature logarithmic gradient profile around $s≈0.55$ has also remained relatively constant in amplitude.

We now investigate the density profile evolution at quasi-steady state. Figure 16 shows the flux-surface- and time-average over $cst/a∈[1246,1266]$ (colored solid) density and its logarithmic gradient profiles for the ions and electrons contributed by $f0+δf$ under the adaptive scheme. For comparison, the density values at the previous time window $cst/a∈[592,613]$ (colored dashed) are shown in Fig. 16(a). As the ion gyrocenter flux decreases approaching quasi-steady state, so does the change in ion gyrocenter density profiles given that the simulations include no particle sources, Fig. 16(a). The adaptive cases with N_p = 128M and 64M seem to have converged density profiles everywhere except near the magnetic axis, due to above-mentioned ion under-sampling of phase-space. The difference in the profiles between the two simulations with different N_p is however very small (5%). Though not shown, the difference between ion gyrocenter and electron density for the adaptive cases (blue and green) remains approximately the same between the two time-averaging windows $cst/a∈[592,613]$ and $cst/a∈[1246,1266]$⁠, see Fig. 11(b). Taken as a whole, species densities have decreased by 50% from their initial values. This amount of deviation certainly challenges the delta-f PIC constraint of $||δf||/||f||≪1$ in the standard approach without background adaptation. We consider now both ion gyrocenter (solid) and electron (dashed) density logarithmic gradient of the adaptive cases as shown in Fig. 16(b). We first note the significant decrease in 50% from their initial values in the radial region $s∈[0.4,1.0]$⁠. For the region $s∈[0.2,0.4]$⁠, however, there is a local increase in $−〈R/Ln〉$ as the actual density values decrease. We see once again the corrugation in the electron density logarithmic gradient near the q = 1 flux surface, which is absent for the ions, which has already been observed in Figs. 9(b) and 15(b) for electron temperature logarithmic gradient profiles.