{'fraction_non_alphanumeric': 0.03775595139004495, 'fraction_numerical': 0.018411852838355253, 'mean_word_length': 4.47103825136612, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerics and computation in gyrokinetic simulations of electromagnetic turbulence with global particle-in-cell codes

Electromagnetic turbulence is addressed in tokamak and stellarator plasmas with the global gyrokinetic particle-in-cell codes ORB5 [E. Lanti et al, Comp. Phys. Comm, vol. 251, 107072 (2020)] and EUTERPE [V. Kornilov et al, Phys. Plasmas, vol. 11, 3196 (2004)]. The large-aspect-ratio tokamak, down-scaled ITER, and Wendelstein 7-X geometries are considered. The main goal is to increase the plasma beta, the machine size, the ion-to-electron mass ratio, as well as to include realistic-geometry features in such simulations. The associated numerical requirements and the computational cost for the cases on computer systems with massive GPU deployments are investigated. These are necessary steps to enable electromagnetic turbulence simulations in future reactor plasmas.


Introduction
Realistic simulations of electromagnetic turbulence are of crucial importance to understand and predict behaviour of burning plasmas before such plasmas become experimentally available. Burning plasmas are complex systems with multiple spatial and temporal scales [1]. Electromagnetic turbulence is ubiquitous in such plasmas. It is the basic component of the "scenery" involving the fast-particle dynamics, the global Magneto-Hydrodynamical (MHD) and Alfvénic activity, zonal flows, and transport. The saturation of the electromagnetic turbulence is a consequence of a complex interplay between these components. A single numerical first-principle framework, based on the global gyrokinetic electromagnetic formulation and including self-consistently all of the parts of the problem, is desirable.
There are a number of gyrokinetic codes which develop capabilities for global simulations in electromagnetic regimes, for example GEM [2], GTC [3,4], XGC [5], GENE [6], GKNET [7], GKW [8]. In this paper, we employ the well-established gyrokinetic particle-in-cell (PIC) codes ORB5 [9] and EUTERPE [10]. The ORB5 code has recently been refactored to allow simulations on heterogeneous (GPU) systems. A very good scaling of the code has been reported in Ref. [11] on various GPU systems and including electromagnetic perturbations. Our paper extends these results to massive production simulations of electromagnetic turbulence deploying GPUs on the Marconi100 computing system (CINECA). We believe that it is a step forward in the adoption of this new HPC technology by the gyrokinetic PIC simulation community. Electromagnetic simulations are known to be very challenging for the gyrokinetic particle-in-cell codes because of the numerical stability issues related to the cancellation problem [12,13]. Such simulations are also very time consuming since the fast electron dynamics has to be resolved. In ORB5 and EUTERPE, we address the numerical stability problem using the pullback mitigation technique [14,15] for the cancellation problem. Very long simulation times normally required when the electron dynamics is resolved are substantially accelerated in the case of ORB5 on the GPUs of the Marconi100 computer system. This paper reports on the first massive usage of the GPU accelerators in ORB5 production runs.
Electromagnetic turbulence has already been addressed using ORB5 in Ref. [16]. While being a groundbreaking effort, this work had a number of limitations. For example, the plasma beta β = 0.1% and the machine size L x = 2r a /ρ i = 350 were quite small. Here, r a is the minor radius of the tokamak and ρ i is the characteristics ion gyro-radius. Also, the geometry was limited to a large-aspect-ratio circular crosssection tokamak and the typical mass ratio was reduced to m i /m e = 200 (although a limited number of simulations using larger mass ratios were also performed). However, it is known [12,13] that the numerical error caused by the cancellation problem scales as δA ∼ L 2 x β m i /m e . Also, the more complex shaping of the magnetic surfaces in realistic tokamak and stellarator geometries can complicate the cancellation problem [17]. Thus, it still remains to be proven that the gyrokinetic particle-in-cell simulations of electromagnetic turbulence are possible under more challenging, realistic, conditions. This demonstration is the main purpose of our present paper. In the following, we will show that simulations of electromagnetic turbulence using more realistic parameters and in more complex geometries are possible with ORB5 [9] and EUTERPE [10] in their present form and on currently available computer systems. We have chosen ITER and W7-X geometries for the first assessment of the electromagnetic turbulence in the shaped plasmas using the global gyrokinetic particle-in-cell codes.
In the finite-beta regimes, considered in this paper, the plasma experiences a transition from the Ion Temperature Gradient driven (ITG) to the Kinetic Ballooning Mode (KBM) turbulence which becomes dominant when beta increases. An important result of this work is the demonstration of numerically stable global particle-in-cell simulations of the electromagnetic turbulence in the KBM regime. This is a novel result for global fully-gyrokinetic PIC codes. Previously, the global simulations of the ballooning and KBM-type turbulence have been performed using MHD [18] and twofluid [19] codes. Most of the gyrokinetic studies [20,21,22,23] have used the flux-tube approximation and Eulerian discretization of the nonlinear gyrokinetic system. There are only two examples of global gyrokinetic KBM turbulence simulations, reported so far: Ref. [7] using the Eulerian code GKNET and Ref. [4] where the gyrokinetic PIC code GTC has been used. The principle difference of Ref. [4] comparing to our paper is the treatment of the electrons. In Ref. [4], the electrons are treated in the fluid approximation ("fluid-kinetic hybrid electron model") whereas this papers uses the fully gyrokinetic description for all species. In the fluid-electron approximation, the cancellation problem is not an issue, in contrast to our simulations where the cancellation problem must be solved.
The novel simulations reported in this paper have been made possible by a combination of several recent developments. First, the pullback mitigation of the cancellation problem implemented in EUTERPE [14] and more recently in ORB5 [15] has stabilized numerically electromagnetic simulations and relaxed considerably the resolution requirement for such simulations. Second, the GPU enabling of ORB5 [11] made it possible to run the large number of small time steps, needed to resolve the electron dynamics, in a reasonable simulation time since the markers can now be pushed on the GPUs. The third aspect, making this paper possible, is a large amount of computing time available for ORB5 and EUTERPE simulations in the frame of PRACE 21st Call for Project Acces. In particular, all ORB5 simulations reported in this paper have been perofmed on the GPU computing systems Marconi100 installed very recently at CINECA using the computing time from the PRACE project. In this paper, we will not be able to clarify every aspect of the global electromagnetic turbulence physics which is a quite complex topic at the frontier of the ongoing research. However, we hope to contribute to this research providing an insight into the details of the simulations which may help future progress.
The structure of this paper is as follows. In Sec. 2, we present, for completeness, the equations solved by the code. In Sec. 3, simulations are discussed for the large-aspectratio circular cross-section tokamak [16], for the down-scaled ITER [24], and for the down-scaled Wendelstein 7-X (W7-X) configurations. Conclusions are made in Sec. 4.

Equations and discretization
In this section, we present for completeness equations solved by the codes. We use the mixed-variable formulation of gyrokinetic theory. This formulation is based on the splitting of the magnetic potential into the Hamiltonian and symplectic parts A = A (h) + A (s) , see Ref. [14] for details.
The equations include the gyrokinetic Vlasov equation: the equations for the gyro-center orbits: an equation for ∂A (s) /∂t (Ohm's law, see [14]): and the field equations (here we use the notation as in Ref. [14]): For the cancellation problem, the pullback mitigation [14] is applied.
(i) At the end of each time step, the splitting of the magnetic potential is redefined: (ii) As a consequence, A (new) (t i ) = 0 (iii) The new mixed-variable distribution function and the parallel velocity become [25]: The linearized version of the pullback transformation is: The velocity is not transformed when the linearized pullback is applied and the approximation b * = b * 0 is used, see Eq. (6).
(iv) The nonlinear mixed-variable system of equations (4)- (10) is explicitly solved at the next time step t i + ∆t in a usual way, but using the transformed coordinates and the transformed distribution function as the initial conditions.
Details of the pullback implementation have been described in Ref. [15] for ORB5 and in Ref. [14] for EUTERPE. The nonlinear pullback transformation has been introduced in Ref. [25] using the GYGLES code. The nonlinear equations of motion solved in GYGLES and described in Ref. [25] are considerably more complicated than the equations implemented in ORB5 [15] and EUTERPE [14] where the nonlinear terms proportional to A (h) (e. g. the terms proportional to A A (s) ) have been omitted for simplicity. It has been seen numerically (using GYGLES) that these terms do not affect the nonlinear results. The reason for that is the smallness of A (h) which is set to zero at the end of every time step. Another approximation used in ORB5 is which amounts for a small geometric correction normally unimportant. These small omitted terms may however play certain role at longer simulation times. This subject may need further assessment in future. In this present paper, we limit our formulation to keeping the dominant nonlinear terms, proportional to b×∇ A (s) ·∇ φ and corresponding to the nonlinear magnetic flutter. This nonlinearity arises in the equations of motion from the perturbed term kept in b * , Eq. (6). In the next Section, we show that these larger nonlinear terms are of crucial importance for the turbulence saturation at a higher beta.

Increasing plasma beta
We consider geometry and parameters described in Ref. [16]: a tokamak with aspect ratio 10, machine size L x = 350, circular cross-section, the safety factor q(ρ) = 1.8503 − 0.074096ρ + 0.95318ρ 2 − 5.703ρ 3 + 7.5779ρ 4 − 0.19182ρ 5 − 1.8253ρ 6 (here, ρ is the minor radius of the circular flux surfaces), and the plasma profiles: with s = ψ/ψ a , ψ the poloidal magnetic flux, ψ a the poloidal magnetic flux at the plasma edge, s 0 = 0.5, κ n = 0.3, κ T = 1.0, and ∆ T = ∆ n = 0.208. The ubiquitous nonlinear generation of the small scales in the phase space (filamentation) is controlled with the Krook operator. We have used the standard noise-correction Krook operator described in Ref. [26] for ORB5 and further detailed in Ref. [27] for EUTERPE. Following [16], γ Krook = 10 −4 ω ci is applied for the ions and for the electrons with only the zonal-flow conserving correction applied in ORB5 simulations, see Ref. [26] for further implementation details. Keeping all these parameters unchanged, we increase  the plasma beta from β = 0.1% to a more realistic β = 1.6%. In contrast to Ref. [16], we use the nonlinear pullback transformation and do not include the fast particles. In Fig. 1(a), the growth rate normalized to the ion cyclotron frequency ω ci is shown as a function of β for m i /m e = 200. One can see how the growth rate decreases with beta until a certain threshold is reached beyond which the growth rate begins to strongly increase. This is a typical ITG-to-KBM transition pattern [7,5]. Here, the growth rate is computed from the linear phase of the perturbed energy flux evolution, see below for more details. In Fig. 1(b), the growth rate is plotted for different mass ratios including the realistic one for β = 0.8%. One sees that the system is somewhat stabilized when the mass ratio increases. One sees that γ Krook /ω ci = 10 −4 is not much smaller than the growth rate shown in Fig. 1. We correct it to γ Krook = γ lin /10 in the following.
The radial heat flux density is normalized to the gyro-Bohm units Q GB = n e T e c s (ρ s /a) 2 with the sound speed c 2 s = T e /m i and the ion sound gyro-radius ρ s = c s /ω ci . To compute the fluxes and the perturbed profiles, the plasma volume is divided into N b radial bins, is computed using the expressions: The mean value of Q (j) s is plotted in Fig. 2. Note that this is the total radial heat flux density with the convective part subtracted. The growth rate shown in Fig. 1 is 1s only which is assumed to evolve as Q (j) 1s ∼ exp(2γt) during the linear phase on the evolution. The evolving plasma temperature and density are computed using the expressions: In Fig. 3, the temperature relaxation is shown for (a) the ITG regime and (b) the KBM regime. The temperature is computed using Eq. (20). One sees that the temperature relaxes much stronger in the KBM case. Accordingly, the ion heat flux evolution, plotted in Fig. 4(b), shows a sequence of the broad relaxation events in the KBM case whereas the turbulence is more localized around the unperturbed temperature gradient in the ITG case, Fig. 4(a). This localization is a consequence of the weaker temperature relaxation for the ITG turbulence observed in Fig. 4(a). This broader radial distribution of the KBM turbulence is also reflected in the shearing rate roughly estimated here as α E ≈ (1/L x ) 2 ∂ 2 φ 00 /∂s 2 with φ 00 the "zonal" component of the electrostatic potential. The shearing rate is shown in Fig. 5 for (a) the ITG regime and (b) the KBM regime. We normalize the shearing rate to the corresponding linear growth rate γ lin shown in Fig. 1(a). One can see that the distribution of α E is much broader in the KBM case although its value is smaller than in the ITG case.
The mode structure evolution of the ITG turbulence, Fig. 6, shows a typical pattern of the turbulent eddies decorrelated by the zonal flow. In contrast, the nonlinear phase      of the high-beta case is dominated by a global mode. The poloidal periodicity of this structure corresponds to the beta-induced Alfvén Eigenmode (BAE) instability considered in Ref. [28] which has the toroidal mode number n = 5 and the poloidal mode number m = 9. The fast ions have been included in Ref. [28] but it can be shown that this BAE mode is linearly unstable even in the absence of the fast particles destabilized by the background plasma gradients with the growth rate increasing with the plasma beta. At the lower beta, the global structure with the BAE periodicity appears when the fast particles (driving the BAE) are included [16]. The BAE spectral component is clearly observed in the high-beta spectrum shown in Fig. 9. Here, the absolute value of the electrostatic potential corresponding to a given toroidal mode number and averaged over all radial positions and poloidal mode numbers is shown. One sees how the spectrum, initiated at the mode numbers 20 < n < 30 evolves into the spectrum dominated by n = 0 and the toroidal mode numbers corresponding to the BAE, the principal mode number n = 5 and its nonlinear harmonics n = 10 and n = 15. Interestingly, the low-beta spectrum, shown in Fig. 8, also evolves into the state which clearly contains the BAE component although the mode numbers characteristic to the ITG turbulence, 20 < n < 30 can also be clearly seen. This is in contrast to the high-beta case Fig. 9 where this part of the spectrum is suppressed.
We will conclude this section with some observation concerning the numerical peculiarities of the high-beta electromagnetic simulations. First, we address the role of the "nonlinear pullback" in the turbulence saturation. In Fig. 10(a), it can be seen how a high-beta simulation crashes in the nonlinear phase if the linearized pullback approximation, Eq. (14), and the approximation b * = b * 0 in Eq. (6) are used. In contrast, these approximations have a marginal effect in the low-beta case (electromagnetic ITG regime) considered in Ref. [16]. It is important to understand which of these two approximations is responsible for the crash seen in Fig. 10(a). For this purpose, we run a simulation which employs the linearized transformation, Eq. (14) but uses all terms in Eq. (6), including the flutter perturbation (the term proportional to ∇ A (s) ).
In Fig. 10(b), we observe that this simulation runs without any numerical instability.  The linearized pullback transformation, Eq. (14) applied in Fig. 10(b), does not seem to make any substantial difference. This indicates importance of the nonlinear magnetic flutter for the Alfvénic turbulence saturation at higher beta.
The second numerical point to address here is related to the radial resolution necessary for the electromagnetic turbulence simulations. It was actively discussed in the early phase of the numerical electromagnetic gyrokinetics that resolution of the electron skin depth is needed for such simulations to be numerically stable [29,30]. Such a requirement would make electromagnetic simulations very difficult, in particular for large machines such as ITER, where the radial resolution N s ≈ 40 000 would be needed. Fortunately, in our high-beta simulations we see that resolving the collisionless electron skin depth is not needed for the numerical stability. All the simulations shown so far had the radial resolution N s = 256 which is significantly below the electron skin depth resolution. For example, to resolve the skin depth for β = 1.6%, N s = 1250 would be required. In Fig. 11(a), the heat flux is plotted for for β = 1.6% and different radial resolutions, including the case with the electron skin depth being resolved. One sees that the heat fluxes are very close during the overshoot and comparable later for all radial resolutions considered. In these simulations, the total number of markers did not change although it would have to be increased when using the larger grid resolutions. The problem was in the GPU memory limitations prohibiting increasing the marker resolution for the larger radial grids. The expected consequence is the increase in the "noise" level shown in Fig. 11(b). The signal corresponds here to the electrostatic potential inside the diagonal Fourier filter [31] with the half-width ∆m included in the simulations (e. g. to push the markers). The noise corresponds to the electrostatic potential outside this filter. Normally, ∆m = 5 except in Fig. 11, where ∆m = 3.
The noise generation is the third numerical aspect to mention here. In Fig. 11(b), we observe that the noise level is not negligible although one cannot say that the simulation is dominated by the noise. For the typical parameters used here, the radial grid resolution N s = 256, the number of the ion markers N i = 200 million, and the number of the electron markers N e = 500 million, the noise level is about 5% in the final phase of the simulation (it is much smaller at the beginning of the nonlinear phase). The noise level can be reduced by increasing the number of the markers or increasing the decay rate of the Krook operator (affecting the linear growth rate). We can reduce the noise level even for γ Krook = γ lin /10 down to about 4% by increasing the total number of the markers to N i = 400 million and N e = 1.2 billion. To run this case on Marconi100 GPUs, we needed 64 nodes. We shall note that this paper represents the first attempt to simulate electromagnetic turbulence in the high-beta regime using ORB5. Previous works [32,26,27] dedicated to the noise control in ORB5 simulations used the adiabatic-electron electrostatic approximation. It is clear that more effort is needed to extend these earlier results to the electromagnetic regime. These effort shall address alternative noise-control schemes (such as the "quadtree" approach [27]), conservation properties of the noise-control techniques in the electromagnetic regime, effects of the energy and particle sources, as well as dependence of the noise level on the physical parameters such as the plasma beta, the mass ratio, the machine size L x , and the magnetic equilibrium. This program goes beyond the scope of this present paper.

Increasing machine size and mass ratio
Another limitation of Ref. [16] was the rather small machine size measured in gyroradii L x = 2r a /ρ i = 350 (for ITER L x > 1000). Increasing L x is computationally more challenging than increasing β since it requires more resolution. We typically include about L x /3 poloidal modes in our simulations, covering roughly the range 0 ≤ k ⊥ ρ i ≤ 1. The number of the toroidal modes are defined by the condition of small k for a given poloidal mode number (i. e. small Landau damping) which implies m max ≈ q max n max where q max is the maximum value of the safety factor. These broad turbulence spectra require a spatial grid fine enough to resolve all the modes included (typically four grid points per Fourier mode both in the poloidal and in the toroidal directions). Finally, the radial resolution and the number of the markers for all species need to be adjusted (increased) accordingly to fulfill the condition 0 ≤ k ⊥ ρ i ≤ 1 and to have enough markers per Fourier mode and per radial grid cell. These resolution requirements may be in conflict with the available resources, such as the GPU memory, limiting both the number of the markers and the number of the grid points.
In Fig. 12, the heat flux and the signal-to-noise ratio are shown for β = 0.8% and L x = 480. One sees that the noise level is increased in the case with larger L x . This can be expected since the number of the Fourier modes to be resolved is also larger in this case: 832 Fourier modes for L x = 480 compared to 400 Fourier modes for L x = 350 considered in the previous section. We manage to reach the noise level of 7% in the late nonlinear stage in our largest simulation for L x = 480 using 64 nodes on Marconi100. The resulting electrostatic potential evolution is shown in Fig. 13. Here, one can see again the characteristic global structure with the BAE periodicity emerging from the turbulent background and similar to the structure observed in Fig. 7.
We now turn our attention to increasing the mass ratio from m i /m e = 200 used so far in this paper, except in Fig. 1(b), to higher values. An obvious difficulty of this parameter adjustment is the need to decrease the time step while keeping the simulation length in the physical time approximately the same, leading to a longer duration of the simulation. For the simulations shown in Fig. 1(b), we have observed that the scaling of the time step with m i /m e is sufficient. In Fig. 14, the heat flux and the signal-to-noise ratio are shown for β = 0.8%, L x = 350, and the mass ratio m i /m e = 1000. One sees that the noise is at the level of 1% in this case. The electrostatic potential evolution is shown for this case in Fig. 15. In this case, we can see the BAE-type structure developing similar to Figs. 7 and 13. Also in the evolution of the toroidal spectrum, shown in Fig. 16, one can see how the initial peak at 20 < n < 30 is replaced by the distinct BAE component at n = 5 and its nonlinear harmonics at n = 10 and n = 15. In Fig. 17, the temperature and the density profile relaxation, computed using Eqs. (20) and (21), are shown for β = 0.8%, L x = 350, and m i /m e = 1000. One sees that the relaxation is not very strong. The normalized shearing rate and the radial heat flux evolution are plotted in Fig. 18 showing the turbulence spreading in the regions of the large temperature gradients.

Simulations in shaped plasmas
In this subsection, we attempt to relax the limitation of the simplified tokamak geometry considered in Ref. [16]. The main motivation for this attempt is the observation that the cancellation problem may become more complicated in shaped plasmas [17]. First, we follow Ref. [24] and consider the ITER geometry with the appropriate magnetic field structure, aspect ratio, and the safety factor profile. This configuration is, however, scaled down to the parameters used in the majority of the runs of Ref. [16], meaning      [16]. One sees that the turbulence reaches a saturated state.
small L x = 350, β = 0.001, and m i /m e = 200. The temperature and density profiles employed in the simulation are shown in Fig. 19. One sees that the profiles relax as a consequence of the turbulent dynamics. In Fig. 20, the heat flux density and the structure of the electrostatic potential in the poloidal cross-section of this down-scaled ITER are shown. One sees that the turbulence is in the saturated state with the mode structure showing turbulent eddies broken by the zonal flows, qualitatively similar to the large-aspect-ratio tokamak, shown in Fig. 6(a). For a long-time numerical stability, it was necessary to increase the width of the diagonal poloidal filter using 21 poloidal modes in the filter (instead of the 11 modes used before). As a consequence, the noise level is about 20% in the nonlinear stage for the marker resolution available in the simulation, although it does not collapse into a numerical instability for quite a long time. A substantial increase in the computational effort will be needed in the future. Finally, we relax the limitation of the axisymmetric ambient magnetic field and consider the stellarator geometry of Wendelstein 7-X (W7-X). In contrast to all previous simulations, we use EUTERPE code and the SKL Irene supercomputer (TGCC) for the stellarator runs. As in the ITER case considered above, we down-scale the real W7-X to the parameters similar to those used in Ref. [16]: the number of ion gyro-radii in the simulation box L x = 200, β = 0.14%, and m i /m e = 200. The Krook operator with the relaxation rate γ Krook /ω ci = 10 −5 is used for the noise control. The diagonal filter is applied with the half-width ∆m = 17. The simulation is initialized with a random noise. The flat density profiles and the ITG-unstable temperature profiles are chosen for both the ions and the electrons. The resulting evolution of the ion temperature profile and the mode structure in the saturation phase are shown in Fig. 21. One sees that the turbulence is located in the center of the simulation domain where the temperature gradient has its maximum. The zonal flow generated by the turbulence has its maximum, as expected, in the region of the strong ITG activity. As a consequence, the weaklyelectromagnetic turbulence enters a saturation regime. The temperature gradient is nonlinearly relaxed by the turbulence. The average heat flux evolution and the signalto-noise ratio are plotted in Fig. 22. One sees that the noise becomes quite large in the late phase of the simulation although it does not dominate the signal. Increasing the marker resolution would mitigate this issue but it would be very time-consuming (120 hours on 1920 CPU cores were needed for the run shown in Fig. 22). Still, it is encouraging that the electromagnetic PIC simulation of this type in the global stellarator geometry has run for quite a long time without a collapse. In future, a careful study of the noise correction and its effect on the stellarator zonal flow and turbulence, similar to Ref. [27], will be needed in the electromagnetic regime.

Conclusions
In this paper, we have considered electromagnetic turbulence in tokamak and stellarator plasmas. The large-aspect-ratio tokamak and down-scaled ITER and W7-X geometries have been considered. In the large-aspect-ratio tokamak, we have shown that the extension of the results of Ref. [16] to larger plasma beta, larger machine size, and larger mass ratio is possible and can be affordable on existing computer systems with a careful tradeoff between noise level and simulation cost.
Simulations at a larger machine size appear to be the most challenging in terms of the computational cost since both the grid and the total number of markers have to be increased. As a consequence, a large number of nodes and GPUs is needed. Also, the usage of the GPU memory can become volatile and unstable when approaching the memory limit for any of the GPUs deployed.
For larger beta, including the perturbed magnetic flutter in Eq. (6) is mandatory to avoid numerical instabilities in the nonlinear phase. In terms of computational resources, the large-beta simulations are normally rather affordable since the grid and the marker resolution can often be the same as in the low-beta case. For larger beta, the nonlinear phase is dominated by a global mode with the periodicity of the BAE instability [28].
Increasing the mass ratio requires longer simulations since the time step has to be reduced. This also makes the overall computation cost higher but there are less issues related to the number of the GPUs and the stability of the GPU memory providing that the grid resolution and the number of the markers do not change.
In the realistically-shaped down-scaled ITER geometry [24], we have demonstrated numerically-stable saturated electromagnetic turbulence in the regime and for the parameters similar to Ref. [16]. However, the numerical noise was not small in this case for the marker resolution available in the simulation. The noise was also an issue in the stellarator (W7-X) electromagnetic turbulence in the nonlinear stage of the simulation. We will need a much larger marker resolution to suppress the noise in the realistic geometries. Presently, this requirement is in a conflict with the limitations imposed by the GPU memory. In future, substantial computational resources and an increased attention to the noise suppression techniques will be needed to simulate electromagnetic turbulence under realistic conditions with the gyrokinetic PIC codes.