{'fraction_non_alphanumeric': 0.03875213623961862, 'fraction_numerical': 0.014646358648397445, 'mean_word_length': 4.378295573651536, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The GBS code for the self-consistent simulation of plasma turbulence and kinetic neutral dynamics in the tokamak boundary

A new version of GBS (Ricci et al. Plasma Phys. Control. Fusion 54, 124047, 2012; Halpern et al. J. Comput. Phys. 315, 388-408, 2016; Paruta et al. Phys. Plasmas 25, 112301, 2018) is described. GBS is a three-dimensional, flux-driven, global, two-fluid turbulence code developed for the self-consistent simulation of plasma turbulence and kinetic neutral dynamics in the tokamak boundary. In the new version presented here, the simulation domain is extended to encompass the whole plasma volume, avoiding an artificial boundary with the core, hence retaining the core-edge-SOL interplay. A new toroidal coordinate system is introduced to increase the code flexibility, allowing for the simulation of arbitrary magnetic configurations (e.g. single-null, double-null and snowflake configurations), which can also be the result of the equilibrium reconstruction of an experimental discharge. The implementation of a new iterative solver for the Poisson and Amp\`ere equations is presented, leading to a remarkable speed-up of the code with respect to the use of direct solvers, therefore allowing for efficient electromagnetic simulations that avoid the use of the Boussinesq approximation. The self-consistent kinetic neutral model, initially developed for limited configurations, is ported to the magnetic configurations considered by the present version of GBS and carefully optimized. A new MPI parallelisation is implemented to evolve the plasma and neutral models in parallel, thus improving the code scalability. The numerical implementation of the plasma and neutral models is verified by means of the method of manufactured solutions. As an example of the simulation capabilities of the new version of GBS, a simulation of a TCV tokamak discharge is presented.


Introduction
Understanding the physical phenomena taking place in the boundary of fusion plasmas is of fundamental importance for the operation and design of future devices, such as ITER and DEMO.
In fact, the dynamics in the plasma boundary controls the overall performance of a fusion reactor, by setting the boundary conditions for the core. In addition, by regulating the plasma-wall interaction, the boundary region controls the plasma refueling, heat exhaust and impurity dynamics [1].
The plasma dynamics in the boundary region is strongly nonlinear and characterized by the presence of phenomena occurring on a wide range of length and time scales, from the ion Larmor radius to the machine size, and from the gyro-motion to the time scales of a discharge. In contrast to the tokamak core, large-amplitude turbulent fluctuations do not allow decoupling the fluctuating and background quantities. In addition, the plasma boundary, composed of the edge, where the magnetic field lines lie on closed and nested flux surfaces, and the scrape-off layer (SOL), where the magnetic field lines intersect the tokamak wall, is usually characterized by a complex magnetic field geometry that features one or more nulls of the poloidal magnetic field. The wide range of time and spatial scales, large amplitude fluctuations, and a complex magnetic geometry make the plasma boundary extremely challenging to model and simulate.
GBS is a first-principles, three-dimensional, flux-driven, global, turbulence code that evolves the Braginskii equations [17] in the drift limit (see, e.g., Ref. [18] for the derivation of the driftreduced Braginskii equations). The GBS code was initially developed to study plasma turbulence in basic plasma devices. The initial version of GBS evolved the two-dimensional plasma dynamics in the plane perpendicular to the magnetic field [19,20], mainly studying ideal-interchange turbulence in simple magnetized plasma (SMT) configurations, such as TORPEX [21]. Later, GBS was extended to include the direction parallel to the magnetic field, using a field-aligned coordinate system, and was used to perform global simulations in SMT configurations [22,23,24] and linear devices [25], such as LAPD [26]. In 2012, a new version of GBS was developed to simulate plasma turbulence in the SOL of tokamak devices in limited magnetic configurations [27]. An electrostatic model in the cold ion limit was considered. Moreover, the Boussinesq approximation was applied to simplify the numerical implementation of the divergence of the polarisation current.
A second version of GBS was developed and presented in Ref. [28]. The physical model was improved by adding the ion temperature dynamics [29] and electromagnetic effects [30]. The plasma model was coupled to a self-consistent kinetic neutral model, leading to the first plasma turbulence simulations of the SOL that self-consistently include the coupling to the neutral dynamics [31]. The Boussinesq approximation was relaxed and the code parallelisation was substantially improved by means of a three-dimensional Cartesian communicator. The implementation of the plasma model was verified by using the method of manufactured solutions (MMS), described in Ref. [32].
Finally, a non-field-aligned coordinate system was introduced in GBS to simulate complex magnetic geometries including one or more X-points and a third version of the code was reported in 2018, leading to the first GBS simulation of a diverted geometry [11]. The second-order numerical scheme was improved to a fourth-order finite difference scheme. Considering a simplified model with respect to the limited configuration, the diverted version of GBS was electrostatic, made use of the Boussinesq approximation, and did not include the neutral dynamics. Despite the use of a domain with a circular poloidal cross section that avoids the core region, limiting the choice of magnetic geometry, the version of GBS described in Ref. [11] was used to investigate plasma turbulence in both single-and double-null magnetic configurations, shading light on the properties of blobs [33,34].
The goal of the present paper is to describe in detail a new version of GBS that significantly improves the diverted version of GBS reported in Ref. [11]. First, while keeping the same fourthorder numerical scheme and non-field-aligned coordinate system, a rectangular poloidal cross section is implemented, which encompasses the whole plasma volume avoiding an artificial boundary with the core, hence retaining the core-edge-SOL turbulence interplay that is found to play a key role in determining the plasma dynamics of the tokamak boundary [35,36,37]. A new coordinate system is chosen adapted to the new geometry. GBS simulations that leverage this new domain and coordinate system have been recently used to investigate the turbulent transport regimes of the tokamak boundary [38] and to develop a theory-based scaling law of the pressure and density decay lengths in both the near and the far SOL [39]. Ultimately, the domain implemented in the present version allows for more flexibility on the choice of the magnetic configuration, which can also be loaded from an equilibrium reconstruction or a Grad-Shafranov solver. The possibility to simulate a magnetic configuration loaded from an equilibrium reconstruction has recently allowed for a direct comparison between GBS simulations and TCV experiments [40]. Moreover, leveraging the flexibility of the numerical scheme, the first turbulence simulations of snowflake magnetic configurations were performed [41]. Second, the physical model of Ref. [11] is extended by adding electromagnetic effects, and the self-consistent kinetic neutral model, initially developed in the context of the limited version of GBS [28,31], is ported to the present configuration, with an improved set of boundary conditions. Third, both the plasma and neutral implementations are carefully refactorized and optimized. In particular, the implementation of a new iterative solver for the Poisson and Ampère equations is introduced, leading to a remarkable speed-up of the code, therefore allowing for efficient electromagnetic simulations that avoid the use of the Boussinesq approximation. Fourth, the GBS parallelisation scheme is improved here by allowing the evolution in parallel of the plasma and neutral dynamics and leading to an improvement of the code scalability. In addition, the relatively simple numerical scheme used in GBS allows for an effective parallelisation of GBS through domain decomposition applied to all three coordinates and implemented with the Message Passing Interface (MPI).
In the present work, the numerical implementation of both the plasma and neutral models is also verified by means of the MMS for the first time. In addition, taking advantage of the improvements of GBS brought to the present version, we present the first GBS electromagnetic simulation of a single-null TCV equilibrium, including the self-consistent evolution of the neutral dynamics.
The present paper is organized as follows. After the Introduction, Sec. 2 describes the physical model implemented in GBS. In Sec. 3, we summarize the numerical implementation and optimization of the plasma and neutral models. The code verification is reported in Sec. 4, while parallelisation scalability tests are described in Sec. 5. Convergence properties are presented in Sec. 6. The first application of the new version of GBS described here is then presented in Sec. 7.
The conclusions follow.

Physical model
The aim of this section is to summarize the physical model implemented in GBS. In particular, we highlight the main differences with respect to the physical model that was evolved by the previous diverted version of GBS, described in Ref. [11].

Drift-reduced Braginskii equations
When the electron mean-free path is shorter than the parallel connection length, λ e L ∼ 2πqR, and small deviations from a Maxwellian distribution function are expected, the use of the Braginskii fluid model [17] to study the plasma dynamics is justified. In addition, turbulence in the tokamak boundary occurs on time scales considerably longer than the ion cyclotron time scale, defined by the ion cyclotron frequency Ω ci = eB/m i , and on spatial scales larger than the ion sound Larmor radius, ρ s = c s /Ω ci , with c s = T e /m i . This allows for the use of the drift-ordering approximation (see, e.g., Ref. [18]) to simplify the Braginskii equations.
In the new version of GBS presented here, we extend the physical model implemented in Ref. [11] by relaxing the Boussinesq approximation, adding electromagnetic effects, and coupling the Braginskii model with an improved version of the self-consistent kinetic neutral model initially developed in the limited version of GBS [28,31]. Then, the plasma model equations take the following form: and are coupled to the Poisson and Ampère equations, In Eqs. (1)- (8), is the scalar vorticity, while U e = v e + eψ/m e is the sum of the electron inertia and the electromagnetic induction contributions. Ampère law, Eq. (8), is solved for magnetic field fluctuations perpendicular to the equilibrium magnetic field. Precisely, the magnetic field perturbation, δ B, is written in terms of the fluctuating parallel component of the vector potential, ψ, i.e. δ B = −∇ × (ψb), where b = B/B is the unit vector of the unperturbed magnetic field and δ B results in a vector perpendicular to B (neglecting the parallel component of δ B is equivalent to exclude the fast compressional Alfvén wave from the dynamics [18]). In addition, we avoid to evolve the externally imposed equilibrium magnetic field, which would require to couple the plasma model to a Grad-Shafranov solver and to self-consistently simulate the plasma current. For this reason, the contribution of the equilibrium parallel current j , computed by averaging the parallel current in the toroidal direction, is subtracted to the total parallel current in Eq. (8). We therefore neglect large-scale variation of the equilibrium magnetic field and we focus on the electromagnetic effects rising from small-scale, small-amplitude magnetic perturbations [18]. The spatial operators appearing in Eqs. (1)- (8) are the E × B convective term, the curvature operator, the parallel gradient, which includes the electromagnetic flutter contribution, and the perpendicular Laplacian, where f is a general scalar function. The implementation of these operators in the GBS coordinate system is outlined in Sec. 2.3. The source terms in the density and temperature equations, s n and s T , are added to fuel and heat the plasma. The gyroviscous terms are defined as where η 0i = 0.96nT i τ i and η 0e = 0.73nT e τ e . The numerical diffusion terms, D f ∇ 2 ⊥ f , are added for numerical stability.
Plasma and neutral dynamics are coupled through the ionization, recombination and chargeexchange processes as well as elastic electron-neutral collisions. These are described by the Krook operators with collision frequencies defined as and where v e and v i are the electron and ion velocities, σ iz , σ rec , σ en and σ cx are the ionization, recombination, elastic electron-neutral, and charge-exchange collision cross sections. The effective reaction rates, vσ , are taken from the OpenADAS database [42]. In Eq. (5), E iz is the ionization energy. While Eqs. (1)- (8) are presented in physical units, they are implemented in GBS in dimensionless form. In the following, we use physical units, but we present the simulation results using GBS dimensionless units. The density, n, is normalized to the reference density n 0 . The electron and ion temperatures, T e and T i , are normalized to the reference values T e0 and T i0 , respectively. The electron and ion parallel velocities, v e and v i , are normalized to the reference sound speed c s0 = T e0 /m i . The magnetic field is normalized to its modulus B 0 at the tokamak magnetic axis. The electrostatic potential, φ , is normalized to T e0 /e and ψ is normalized to ρ s0 B 0 , with ρ s0 = c s0 /Ω ci the reference ion sound Larmor radius. Perpendicular lengths are normalized to ρ s0 and parallel lengths are normalized to the tokamak major radius, R 0 . Time is normalized to R 0 /c s0 . By normalizing Eqs. (1)-(8), the following dimensionless parameters that regulate the system dynamics are identified: the normalized ion sound Larmor radius, ρ * = ρ s0 /R 0 , the ion to electron reference temperature ratio, τ = T i0 /T e0 , the normalized electron and ion parallel thermal conductivities, and the reference electron plasma β , β e0 = 2µ 0 n 0 T e0 /B 2 0 , and the normalized Spitzer resistivity, ν = e 2 n 0 R 0 /(m i c s0 σ ) = ν 0 T −3/2 e , with σ = 1.96 n 0 e 2 τ e m e n = 5.88 4 √ 2π and where λ is the Coulomb logarithm.

Kinetic model for neutral atoms
The effect of neutral particles on plasma turbulence has been approached by modelling the neutral dynamics with fluid and kinetic models. For example, two-dimensional fluid models that include the neutral gas continuity equations, accounting for neutral gas ionization and charge exchange processes, have been used to investigate the role of the neutral gas in the SOL [43,44]. More recently, three-dimensional turbulent simulations of the tokamak boundary, coupling fluid plasma and fluid neutral models, have been carried out with the BOUT++ code [45]. However, since the neutral mean-free path in the SOL can be large, a kinetic approach is often preferred for the evolution of neutral dynamics. The first attempt to couple a kinetic neutral model to a plasma turbulence model was presented in Ref. [46], where the EIRENE kinetic code [47] based on a Monte Carlo algorithm was coupled to the two-dimensional turbulence code TOKAM2D. The first three-dimensional turbulent simulation of the SOL with the self-consistent evolution of a kinetic neutral model was carried out by using GBS in limited geometry [31]. Recently, the three-dimensional turbulence code TOKAM3X has been coupled to the EIRENE code to simulate plasma turbulence and kinetic neutrals in a diverted geometry [48].
We summarize here the kinetic model for neutral atoms implemented in GBS. The model is described in detail in Ref. [31] and extended here to the simulation of diverted configurations to include the contribution of the E × B and diamagnetic fluxes to the ion flux to the wall. We consider a single mono-atomic neutral species represented by a distribution function f n with its dynamics being described by the following kinetic equation, where f i and n i are the ion distribution function and the ion density, respectively. The neutralneutral collisions, which have a lower reaction rate than the charge-exchange and ionization processes in medium size tokamaks or attached conditions, are neglected. The elastic electron-neutral collisions, retained in Eq. (3), are neglected in Eq. (24) because of the electron-to-neutral mass ratio.
The boundary conditions for f n at the wall are derived under the assumption that the impacting neutrals and ions are either reflected or absorbed. If absorbed, the neutral particle is immediately released with a velocity that depends on the wall properties and is independent of the impacting particle velocity. The distribution function of the neutrals flowing from wall to the plasma volume (i.e. neutrals with velocity such that v ·n > 0, withn the unit vector normal to the wall) is therefore given by where α refl is the reflection fraction, assumed the same for neutrals and ions, x b indicates the boundary position, Γ out = Γ out,n + Γ out,i is the sum of the neutral and ion fluxes to the wall and projected in the direction perpendicular to it, v p = v pn is the velocity perpendicular to the boundary, with v p = v ·n, and χ in is the inflowing velocity distribution function given by the Knudsen cosine law, with θ = arccos(v p /v) and T b the wall temperature [31].
In the limit where the turbulent time scale is much longer than the typical time of flight of neutrals, τ turb τ n , the neutral adiabatic approximation can be applied. This corresponds to impose ∂ t f n = 0 in Eq. (24). Moreover, we assume that the neutral mean free path is shorter than typical parallel lengths of the plasma structures. Under these assumptions, the formal solution of Eq. (24) can be obtained by using the method of characteristics [31], where we express a position x in terms of x ⊥ , the coordinate on the plane perpendicular to B, and x , the coordinate parallel to B, r ⊥ is the coordinate along the neutral characteristic defined by r ⊥b denotes the distance along the characteristic from the position x and the wall, v ⊥ is the component of the velocity perpendicular to B, and ν eff = ν iz + ν cx is the effective collision frequency for neutral loss. In the following, we drop the parametric dependencies on t and x to simplify the notation. The volumetric source term in Eq. (27) results from charge-exchange and recombination processes and is given by [31] S where is the ion velocity distribution. The source of neutrals at the wall, δ (r ⊥ − r ⊥b ) f n (x ⊥b , x , v,t), is given by the boundary conditions in Eq. (25). The ion recycling term present in the boundary conditions and the recombination term appearing in S(x ⊥ , v) (see Eqs. (25) and (28)) do not depend on the neutral distribution function and can be computed directly from the plasma quantities. On the other hand, the charge-exchange term in S(x ⊥ , v) as well as the reflected and re-emitted neutrals in the boundary conditions (see Eqs. (25) and (28)) depend on n n (x ⊥ ) = f n dv, the neutral density. This suggests to integrate Eq. (27). In fact, Ref. [31] shows that, by integrating Eq. (27), a linear integral equation for n n (x ⊥ ) is obtained, where dA is the infinitesimal area in the poloidal plane D, da b is the infinitesimal length along the boundary ∂ D, and Γ out,n (x ⊥b ) is the neutral flux towards the wall, The contribution to the neutral density due to the ion recycling at the wall and recombination events in Eq. (29), n n[out,i] and n n[rec] , as well as the corresponding contribution to the neutral flux in Eq. (30), Γ out,n[out,i] and Γ out,n [rec] , are defined as Generalizing Ref. [31], the ion flux projected in the direction perpendicular to the wall, Γ out,i , is evaluated by considering the ion parallel flux and the perpendicular fluxes due to the E × B and the diamagnetic drifts The kernel functions, appearing in Eqs. (29) and (30), involve integrals in the velocity space, where we consider the direct path between two points as well as the paths that include one reflection at the wall (we neglect paths with multiple reflections as λ n /L < 1, with λ n the averaged neutral mean free path and L the typical machine size in poloidal plane) and we define, for path = {dir, refl}, , being K 0 (x) the modified Bessel function of the second kind. The vector x ⊥ indicates a position along the path from the source to the We note that a neutral flux can be externally imposed by means of a gas puff that introduces a localized source of neutrals. Similarly, a pumping region on the wall can also be considered. This can be simply implemented by multiplying the kernel functions, K b→p and K b→b , by a recycling coefficient smaller than one.

Differential operators
The differential operators in Eqs. (9)- (12) are written in the (R, ϕ, Z) cylindrical non-fieldaligned coordinate system, where R is the distance from the axis of symmetry of the torus, Z is the vertical coordinate, and ϕ is the toroidal angle. The poloidal cross section has a rectangular shape, particularly suitable for the simulation of the TCV tokamak, and the domain encompasses the whole plasma volume. A representation of the GBS domain is shown in Fig. 1.
The toroidally symmetric equilibrium magnetic field, used to compute the GBS differential operators, is written in terms of the poloidal magnetic flux Ψ, as The poloidal magnetic flux is a function of R and Z, and can be provided as an analytical function, an equilibrium reconstruction, or as the result of a Grad-Shafranov solver. As an example, the magnetic separatrix of a TCV equilibrium reconstruction is shown in Fig. 1.
The spatial differential operators in Eqs. (1)- (8) are written in (R, ϕ, Z) coordinates. They are expanded in the large aspect ratio limit, ε ∼ r/R 0 1, assuming the poloidal component of the magnetic field smaller than the toroidal one, δ ∼ B p /B ϕ 1, with a safety factor at the midplane q ∼ ε/δ of order unity. Only the leading order terms in ε and δ are retained.
In the following, we derive the expression of the operators as they are implemented in GBS. As a first step, we note that, at zeroth-order in ε and δ , the modulus of the magnetic field is constant, where the magnetic field is written as with e R , e Z , and e ϕ the basis vectors. Hence, the Poisson brackets, Eq. (9), can be written as which, in dimensionless units, leads to Neglecting local current and therefore assuming ∇ × B = 0, the curvature operator, Eq. (10), is expanded in ε and δ as The spatial derivatives ∂ Z B 2 and ∂ R B 2 can be determined by using Eq. (45), where Finally, by keeping only the leading order terms in ε and δ , the curvature operator, normalized to 1/(R 0 ρ s0 ), becomes The parallel gradient, Eq. (11), is normalized to 1/R 0 and, in dimensionless units, is given by where The perpendicular laplacian, Eq. (12), can be developed as The second term on the right-hand side of Eq. (54) is one order ε smaller than the first one, which can be written in cylindrical coordinates as and expanded in ε and δ , Thus, the leading order terms of Eq. (56) can be expressed as In addition to the differential operators given by Eqs. (9)- (12), there are three operators appearing in Eqs. (1)-(6) that can be derived from Eqs. (9)-(12) by applying the same ordering in ε and δ : the parallel laplacian, ∇ 2 f , the curvature of the parallel gradient, C(∇ f ), and the parallel gradient of the curvature, ∇ [C( f )]. In dimensionless units, the parallel laplacian is given by the curvature of the parallel gradient by and the parallel gradient of the curvature by

Boundary and initial conditions
The large electric field established on the ρ s scale in the magnetic pre-sheath violates the hypothesis behind the drift-approximation. Consequently, the boundary conditions of GBS are applied at the magnetic pre-sheath entrance. The set of generalized Bohm-Chodura sheath boundary conditions implemented in GBS was originally derived in Ref. [49] in the cold ion limit and extended in Ref. [29] to include warm ions. In this version of GBS, we neglect the correction terms arising from the variation of the density and electrostatic potential in the direction tangent to the wall, and we apply the magnetic pre-sheath boundary conditions at the walls containing the strike points, i.e. the walls where the divertor legs terminate. For instance, in the case of the TCV magnetic equilibrium shown in Fig. 1, magnetic pre-sheath boundary conditions are applied at the bottom and inner walls. The boundary conditions assume the following form: where Λ 0 = log m i /(2πm e ) 3 for hydrogen plasmas, and s denotes the derivative in the direction perpendicular to the wall (s = Z for the top and bottom walls, and s = R for the inner and outer walls). The top (bottom) sign refers to the magnetic field pointing towards (away from) the target plate. In addition to the boundary conditions in Eqs. (61)-(66), we simply consider ψ = 0 at the magnetic pre-sheath entrance. In order to avoid the discontinuity of the parallel velocities at the location where the magnetic field is tangent to the wall, a smoothing function from +c s to −c s is applied to guarantee that v i and v e vary without strong discontinuities. In addition, we impose ∂ s |v i | = 0 and ∂ 2 ss |v i | = 0 in Eqs. (61)- (66) in these smoothing regions where the normal derivative to the wall of |v i | can become negative, leading to unphysical boundary conditions for n, φ and Ω. We highlight that the smoothing affects less than 2% of the wall surface and is at certain distance from the strike points, therefore it has a negligible effect on the overall dynamics.
The boundary conditions in Eqs. (61)- (65) are also applied at the walls that do not contain strike points, while the Dirichlet boundary condition φ = ΛT e /e is used for φ , with Λ = Λ 0 − log 1 + T i /T e . In fact, Poisson equation is ill-defined if a Neumann boundary condition for the electrostatic potential is applied at the four walls of the domain.
As long as the quasi-steady state that is established after a transient is statistically independent of the initial conditions, the initial profiles do not present physical interest. As initial condition for our simulations, we impose φ = ΛT e /e, with T e the initial electron temperature, which is uniform over the entire GBS domain. Also Ω, n, and T i have initial uniform profiles. The electron and ion parallel velocities are properly designed functions to satisfy the boundary conditions v e,i | wall = ±c s 1 + T i /T e with ∂ s v e,i | wall = 0, being s the coordinate perpendicular to the wall. The constraint ∂ s v i | wall = 0 guarantees vanishing Neumann boundary conditions for n, φ , ω, in agreement with the uniform initial profiles of these quantities. Random noise is added to all fields to trigger the turbulent dynamics.

Implementation and optimization
While the time evolution is provided by a fourth-order Runge-Kutta algorithm and the spatial discretization by a fourth-order centered finite differences scheme, similarly to the GBS diverted version of Ref. [11], GBS now relies on an optimized iterative solver for the Poisson and Ampère equations. In addition, the kinetic neutral model, initially developed in the limited version of GBS [28,31] and ported here to the diverted geometry, is significantly optimized.
We show in Fig. 2 the workflow of GBS. Following the initialization of the simulation, plasma and neutral modules are run simultaneously. Each Runge-Kutta substep, used to advance the GBS plasma equations, performs three tasks (see the grey area of Fig. 2). First, the boundary conditions are applied to every plasma quantity and the necessary spatial operators are applied to the plasma fields. Then, the right-hand side of the drift-reduced Braginskii equations is evaluated and Eqs. (1)-(6) are advanced, updating the values of n, T e , T i , Ω, U e and v i . Finally, the Poisson and Ampère equations, Eqs. (7) and (8), are numerically solved to update φ and ψ. After every Runge-Kutta time step (see the blue area of Fig. 2), the plasma module checks if the neutral calculation is completed and updates the value of the neutral-related terms in Eqs. (1)-(6) with the result of the new neutral calculation, if available. In parallel to the plasma evolution, the neutral density, temperature and velocity are computed by the neutral module (red area in Fig. 2). At the beginning of each neutral calculation, GBS computes the reaction rates, Eqs. (15)- (19), by using the current plasma quantities. Then, it evaluates the kernel functions, Eqs. (36)- (39), and finally it computes the neutral density, temperature and velocity by solving Eq. (29). After every neutral calculation, the plasma and neutral modules synchronize, and the values of the plasma density, electron and ion temperatures used to compute the reaction rates is updated. A minimum neutral calculation frequency can be imposed in order to guarantee a good convergence of the simulation results, as detailed in Sec. 6.
The remainder of the chapter is structured as follows. First, we describe the spatial discretization of the plasma equations. Then, we focus on the implementation and optimization of the Poisson and Ampère laws. Finally, we report on the implementation and optimization of the neutral model.

Spatial discretization of the plasma equations
All plasma quantities are evaluated on a uniform Cartesian grid, which discretizes the R, ϕ and Z coordinates, with a size N R × N Z × N ϕ , where N R , N Z and N ϕ are the number of grid points in the radial, vertical and toroidal direction, respectively. The neutral quantities are also discretized on a uniform grid in the (R, ϕ, Z) cylindrical coordinate system, N R × N Z × N ϕ , that is in general coarser than the plasma grid in the poloidal plane (N R < N R and N Z < N Z ), while the neutral and plasma grid resolutions along the toroidal direction are the same (N ϕ = N ϕ ). The grid spacing is denoted as ∆R, ∆Z and ∆ϕ. The grid for n, T e , T i , Ω, and φ , denoted as φ -grid, is staggered in the vertical and toroidal directions with respect to the grid for v e , v i and ψ, which is denoted as v-grid. More precisely, the φ -grid is shifted by ∆Z/2 and ∆ϕ/2, along the vertical and toroidal directions, respectively, with respect to the v-grid. The use of staggered grids prevents the formation of checkerboard patterns that can appear when treating an advection problem with a finite centered difference scheme, as shown in Ref. [50]. In the following, we denote as f i, j,k a quantity at the position (R i , ϕ j , Z k ) in the v-grid and as f i, j+ 1 2 ,k+ 1 2 a quantity at the position The differential operators in Eqs. (10)- (12) are computed as a linear combination of first and second derivatives along R, Z and ϕ, with subsequent interpolation between staggered grids, if needed. For instance, the first derivative along Z is discretized at the fourth-order by means of a 5-point stencil as if both the field and its derivative are evaluated on the v-grid. On the other hand, we have if the field is evaluated on the v-grid and its derivative on the φ -grid. Analogous expressions to Eqs. (67) and (68) hold for fields evaluated on the φ -grid, if their derivatives are evaluated on the φ -grid or v-grid, respectively. A quantity evaluated on the v-grid is interpolated to the φ -grid by using a fourth-order interpolation, and analogous expression is used to interpolate from the φ -grid to the v-grid. The second derivative along Z is given by if both the field and its derivative are evaluated on the v-grid. Analogous expression holds if both the field and its derivative are evaluated on the φ -grid. The Poisson brackets, Eq. (9), are discretized by means of a fourth-order Arakawa scheme [51]. We note that the discretization at fourth-order of the curvature-related contributions in the gyroviscous terms, Eqs. (13) and (14), requires a 7-point stencil because of the presence of second derivatives on the φ -grid of a quantity evaluated on the v-grid and vice versa. In order to use a 5-point stencil for all GBS operators, these derivatives are implemented at the second-order. For example, the second derivative on the φ -grid of a field evaluated on the v-grid is given by

Implementation and optimization of the Poisson and Ampère equations
In Ref. [11], the Boussinesq approximation is applied, thus considerably simplifying the implementation of the Poisson equation, and the electrostatic potential is computed by a direct inversion of the perpendicular laplacian operator, Eq. (12), through LU factorization using the external MUMPS library [52]. Despite being computationally demanding, the LU factorization is carried out once for all at the beginning of the simulation and therefore does not significantly impact the cost of the simulations.
In the version of GBS described here, the Boussinesq approximation is relaxed, requiring the implementation at fourth-order of the with and similarly for ∂ Z (n∂ Z φ ) i, j,k . Since n depends on time, the δ coefficients vary in time, therefore the matrix that discretizes the Laplacian operator has to be assembled and factorized every time step, leading to a dramatic increase of the computational effort, which becomes prohibitively large already for a medium size grid if a LU factorization is used to solve Poisson equation. In fact, an initial profiling of GBS, carried out on a medium size poloidal grid (N R × N Z = 300 × 600), corresponding to the one used to simulate a medium size tokamak such as TCV, shows that, if the MUMPS direct solver is used, more than 50% of the time used to advance the plasma equations is spent in the factorization of the matrix that discretizes Poisson equation. We note that a similar matrix needs to be assembled and inverted every time step to solve the Ampère equation (see Eq. (8)). However, the factorization of the Poisson matrix is computationally more expensive than the factorization of the matrix that discretizes Ampère equation.
In order to avoid the LU factorization and optimize GBS, an iterative solver is used. We rely on the Data Management for Structured Grids (DMDA) of the PETSc library [53] that provides a flexible framework to implement a large number of solvers and preconditioners. By carrying out a scan on different solvers and preconditioners, analysing more than 2000 combinations, we find that the best performance is achieved with the algebraic multigrid preconditioner boomerAMG provided by the external package HYPRE [54] and the deflated GMRES solver [55]. In order to fine tune the preconditioner parameters, we focus on the main options of boomerAMG, i.e. the strong threshold, the maximum number of levels, the coarsen type, and the interpolation type (see Ref. [54] for details). The best performance is achieved by using a strong threshold of 0.25, a maximum of levels of 30, the Falgout coarsen type, and the ext+i interpolation type.
In order to evaluate the performance gain arising from the iterative solver implemented with the PETSc library with this set of optimized parameters, we run a set of GBS simulations considering the N R × N Z = 150 × 300, N R × N Z = 300 × 600 and N R × N Z = 600 × 1200 grids. We refer to these grids as the half-TCV, TCV and double-TCV grids, since N R × N Z = 300 × 600 is the typical poloidal grid used for a simulation of a TCV discharge with toroidal magnetic field of 0.9 T. The time to solution per time step is shown in Fig. 3 for the different grid sizes. Across all sizes considered here, we observe a speed-up of, approximately, a factor of 40 when the iterative solver is used, with respect to a direct solver, reducing considerably the cost of our simulations and making the simulation at the TCV scales possible, otherwise prohibitively expensive. As an aside, we also note that the solution of the Poisson and Ampère equations using the iterative solver and the evaluation of the right-hand side of Eqs. (1)-(6) have approximately the same computational cost in the half-TCV and TCV simulations. On the other hand, more than half of the computational time is spent to solve the Poisson equation in the double-TCV simulation, even using the iterative solver (see Fig. 3). Simulations at scales larger than TCV require therefore further optimization.

Implementation and optimization of neutral model in diverted geometry
We summarize here the implementation and the optimization of the neutral model. The discretization of the kernel functions, Eqs. . (74) The linear system of Eq. (74) can then be solved by inverting the kernel matrix, which, in the limit of large aspect ratio considered here, has [N R N z + 2(N R + N Z )] 2 elements and is inverted for each poloidal plane at every neutral calculation. Since our simulation domain encompasses the whole plasma volume, contrary to the previous neutral version of GBS, the K matrix is a dense matrix with, in principle, all its elements strictly positive. However, since the value of the matrix elements decays exponentially with the distance between the two connected points, becoming negligible when their distance is several neutral mean free paths, we introduce a threshold value below which the matrix element is considered as vanishing. Since the evolution of the neutral dynamics is computationally more expensive than the evolution of the plasma dynamics, we recalculate the neutral quantities over a time interval longer than the time step used to evolve the plasma equation, but still smaller than the turbulence timescale. The solution of the kinetic equation is a two steps process. First, the elements of the K matrix are computed, requiring the evaluation of the kernel functions, Eqs. (36)- (39). In the present version of GBS, the evaluation of K is optimized by improving the code vectorisation and avoiding unnecessary computation. Then, the K matrix is inverted to compute the neutral density by using Eq. (74). The matrix inversion is improved by avoiding an expensive LU factorization.
We first focus on the evaluation of the K matrix. This requires to compute [N R N Z + 2(N R + N Z )] 2 elements, all of them involving the evaluation of one of the four kernel functions defined in Eqs. (36)-(39), a computationally expensive evaluation. In fact, approximately 65% of a neutral step is spent in the computation of the K matrix elements (the remaining time is spent in the solution of the associated linear system). In the version of GBS presented in Ref. [31], each element of the K matrix was evaluated by using the same routine that included logic conditions used to identify points belonging to the plasma volume or the wall. Here, the logic conditions are moved outside the routine that evaluates the elements of K, which is split into four specialized routines, each of them computing one of the kernel functions in Eqs. (36)- (39). Therefore, points belonging to the plasma and the wall are separated at the beginning of the simulation and a whole subset of K is evaluated instead of computing each single element. In addition, the modified Bessel function of second kind involved in the K b→p and K b→p kernels is pre-computed at the beginning of the simulation and later re-used at the cost of a memory access. Finally, the computation of the integral r ⊥ 0 ν eff (x ⊥ )dr ⊥ appearing in the four kernel functions, Eqs. (36)-(39), is improved. This requires first the interpolation of the plasma quantities necessary to compute ν eff along the neutral trajectory and then the numerical evaluation of the integral. The interpolation is performed by using a second-order method that is carefully implemented to enable the compiler auto-vectorisation. In contrast to Ref. [31], here we use the same routine to integrate along the neutral trajectory and over the velocity space. Moreover, integrals can be evaluated by either the (left) rectangular or midpoint rules.
We focus now on the inversion of the K matrix in Eq. (74). We note that the K matrix evolves in time, similarly to the matrix associated with the Poisson and Ampère equations. While a direct solver based on the MUMPS library was used in Ref. [31], requiring an expensive LU factorization, here we implement an iterative solver. Also in this case, we choose to use the framework provided by the PETSc library, opting for the use of the GMRES solver without preconditioner.
In order to evaluate the improvement of performance arising from the new implementation, we analyze the time to carry out one neutral step with various neutral grid sizes. We consider a coarse neutral grid of size N R × N Z = 25 × 50, a medium neutral grid of size N R × N Z = 50 × 100, and a fine neutral grid of size N R × N Z = 100 × 200, which is typical of a simulation of a TCV discharge. In Fig. 4, the time to solution per neutral time step of the implementation that evaluates single elements of the K matrix and uses the direct solver is compared to the one of the implementation that evaluates a whole subset of the K matrix and uses the iterative solver. We note that the speedup arising from the new implementation increases with the size of the system, ranging from a factor of four to a factor of ten.

Verification of GBS implementation
We report here on the verification of the plasma and neutral implementation. We verify the two modules separately, as presented below.

Verification of the plasma model
Similarly to the past versions, we verify the plasma model implementation in the present version of GBS by using the method of manufactured solutions. Here, we briefly describe the MMS. A detailed description can be found in Ref. [32] where the MMS was applied for the first time to verify a plasma turbulence code.
In order to verify the numerical implementation of a model M, the analytical solution s of the model M, which satisfies M(s) = 0, is compared to the numerical solution s h of the discretized model M h , where h is the discretization parameter. The code is successfully verified if the discretization error e h = ||s − s h || → h p as h → 0, with p the order of the adopted discretization scheme. Since the analytical solution of a model equations is in general unknown, the main idea behind the MMS is to manufacture an arbitrary analytical function, u, and evaluate S = M(u). We remark that, in general, S = 0, since u is not the solution of M. On the other hand, the arbitrary function u is solution of the model N, defined as N(u) = M(u) − S = 0. Since S can be analytically evaluated, N h and M h are affected by the same discretization error. Therefore, verifying N h is equivalent to verify M h , with the discretization error given by e h = ||u − u h ||.
In the following, we report for the first time the verification of the implementation in GBS of the plasma model that avoids the Boussinesq approximation and includes electromagnetic effects. We consider a simulation domain with a rectangular poloidal cross section of size L R = 37.5 ρ s0 and L Z = 50 ρ s0 , in the radial and vertical directions, respectively. The verification is carried out with ρ −1 * = 100, ν = 1, β e0 = 10 −4 , τ = 1 and m i /m e = 1. The dimensionless parameters are chosen so that the terms on the right-hand side of Eqs. (1)- (6) are all of the same order of magnitude.
The magnetic field considered for the verification is analytically obtained by solving the Biot-Savart law in the infinite aspect ratio limit for a current density with a Gaussian distribution inside the simulation domain, which guarantees continuity and derivability of the magnetic field, and an additional current filament outside the simulation domain, which produces the X-point. This leads to where EI(x) is the exponential integral function, with I 0 = 40 ρ 2 s0 B 0 , σ 0 = 6.25 ρ s0 , R 1 = 100 ρ s0 , Z 1 = 0, and Z 2 = −40 ρ s0 . The manufactured solutions for the evolved scalar quantities u = n, T e , T e , v e , v i , Ω, φ , ψ are chosen as where A u , B u , C u , D u , E u , F u , α u , β u , and γ u are arbitrary constants whose value is chosen to excite all terms in the right-hand side of Eqs. (1)- (8), ensuring that none of them provides a dominating contribution to the numerical error. The source terms can be computed by substituting Eq. (77) into Eqs. (1)- (8). This process is carried out by using the symbolic calculation offered by the Mathematica software package [56] and the analytical expressions of the source terms are directly converted to the Fortran language of the GBS code. In order to decouple the tolerance error associated with the iterative solver from the discretization error, we first consider the verification of GBS where Poisson and Ampère equations are solved with a direct method. As a second step, we discuss the verification of the iterative solver. We do not include the curvature-related contributions appearing in the gyroviscous terms, Eqs. (13) and (14), since they are implemented at second-order. These terms are verified independently.
Figs. 5 (a) and (b) show the L 2 and L ∞ norms of the discretization errors at t = 0.01. We consider various grid refinements as a function of the discretization parameter h = ∆R/∆R 0 = ∆Z/∆Z 0 = ∆ϕ/∆ϕ 0 = ∆t/∆t 0 , where ∆R 0 , ∆Z 0 , and ∆ϕ 0 is the spacing of the N R × N Z × N ϕ = 256 × 256 × 256 grid, with time step ∆t 0 = 6.25 × 10 −6 . The coarsest grid considered for the verification is N R × N Z × N ϕ = 8 × 8 × 8 with ∆t = 2 × 10 −4 . The order of accuracy, with rh indicating the coarsening of the temporal and spatial grids by a factor of r, is shown in Figs. 5 (c) and (d). Since a fourth-order numerical scheme is used to discretize both space and time, we expect the discretization error to decrease as h 4 when h → 0, in good agreement with the result of Fig. 5 that shows the convergence to p = 4 for both the L 2 and L ∞ norms as h decreases.
We discuss now the verification of the iterative solver. The error affecting the solution of an iterative solver is the combination of the discretization error and the tolerance error, the last depending on the tolerance threshold. We consider the solution of Poisson and Ampère equations and we compare the results of the direct and iterative solvers based on the MUMPS and PETSc libraries. Fig. 6 shows the global and local errors of φ , for Poisson equation, and ψ, for Ampère equation, as a function of the relative tolerance of the iterative solver. We note that the errors decrease with the relative tolerance until a value of 10 −7 . Below this value, the discretization error dominates, independently of the solver tolerance.

Verification of the neutral model
We verify here, for the first time, the implementation of the neutral model. We proceed in two steps. First, we verify the correct implementation of the routine used to perform integrals over a neutral trajectory and over the velocity space, which appear in the kernel functions, Eqs. (36)- (39). Second, we focus on the construction and inversion of the matrix K, Eq. (75), used to evaluate the neutral density.
To carry out the first verification step, we perform unit tests that allow us to verify individually the routines used to numerically integrate and to interpolate along a neutral trajectory. The interpolating routine is tested on an analytical function f , considering various resolutions. Fig. 7 shows the typical discretization error in a point of f domain as a function of the grid resolution. The order of convergence is correctly retrieved as h decreases. The integration routine is tested by comparing the numerical and analytical result of the integral of an analytical function f over an arbitrary interval. The discretization error as well as the order of convergence is shown in Fig. 8. The results of the unit tests point out the correct implementation of the routines used to evaluate integrals along a neutral trajectory.
We focus now on the verification of the correct construction and inversion of the matrix K. For this purpose, we compare the numerical solution of the neutral density, n n , to the analytical manufactured solution, n n,M = A n + B n sin(α n R) cos(β n Z) , where A n , B n , α n , and β n are arbitrary constants. In addition, we consider the manufactured neutral flux as vanishing, i.e. Γ out,n,M = 0. The analytical source term required by the MMS is directly computed by substituting n n,M and Γ out,n,M in Eq. (29), where the contribution from the ion flux is not included to decouple the neutral and plasma modules. We note that Eq. (80) requires the analytical evaluation of the kernel function K p→p (x ⊥ , x ⊥ ), Eq. (36). For simplicity, we choose The choice of these functions allows for the straightforward analytical evaluation of K p→p (x ⊥ , x ⊥ ), resulting in K p→p (x ⊥ , x ⊥ ) = const. As a consequence, all the elements of the matrix associated to K p→p arising from the discretization are equal. The integral over the poloidal plane in Eq. (29), evaluated through the matrix inversion (see Eq. (74)), is discretized using the rectangular integration rule. Therefore, the numerical solution of the neutral density is expected to converge to the manufactured solution with the order of convergence p = 1 as h decreases. The local discretization   error of the neutral density and the corresponding order of convergence are shown in Fig. 9 as a function of h. The expected order of convergence of the neutral density to the analytical solution is retrieved, hence verifying the correct implementation in GBS of Eq. (74).

Parallelisation scalability tests
The use of finite differences on a uniform (R, ϕ, Z) grid allows for an effective implementation of GBS on parallel high-performance computers. Domain decomposition is applied to all three spatial dimensions and implemented with the Message Passing Interface (MPI) for both the plasma and neutral equations. Communication between different processors is carried out by means of ghost cell passing, using standard MPI functions. In addition, the global MPI communicator is split into two independent ones in order to evolve the plasma and neutral equations in parallel. After every neutrals calculation, the neutral and the plasma modules synchronize, i.e. the neutral terms appearing in the plasma equations and the plasma quantities in the neutral module are updated (see Fig. 2). The splitting of the MPI communicators has some important advantages with respect to the serial approach: (i) the neutral density is evolved continuously, updating the neutral sources in the plasma equations at the desired frequency by choosing the number of MPI tasks for the plasma and neutral modules, (ii) the time to solution is reduced, and (iii) the code scalability is improved allowing us to run on a larger number of computing nodes, which is essential on large systems due to the large amount of memory required by the neutral module. We show here that this approach leads to a very efficient strong and weak scalability.
Being independent, the parallelization properties of the plasma and neutral modules are analyzed separately. Focusing, first, on the plasma module, we consider a typical grid size of a TCV simulation, with n pR , n pZ and n pϕ the number of MPI tasks in the R, Z and ϕ direction, respectively. To avoid the demanding communications on poloidal planes required by the Poisson and Ampère equations, the strong scalability test is carried out by only increasing the MPI resources allocated to the ϕ dimension, i.e. by solving whole RZ-planes on a computing node. As shown in Fig. 10 (a), GBS scales with an almost perfect efficiency up to n pϕ = 32, while efficiency decreases to 0.75 at n pϕ = 64 (two toroidal planes per computing node), pointing out the high scalability of GBS especially on the toroidal direction, mainly thanks to the fact that the Poisson and Ampère equations are solved independently on each poloidal plane. The weak scalability of the plasma module is carried out on grids of size N R × N Z × N ϕ = 300 × 600 × 2n pϕ , where n pϕ is the number of MPI tasks in the toroidal direction. The results are shown in Fig. 10 (b), confirming the good scalability properties of the plasma module in GBS. The strong scalability of the neutral module is carried out on a neutral grid typical of medium size simulations such as TCV. Similarly to the plasma scalability, we consider one computing node in the RZ-plane, while increasing n pϕ , starting from the minimum number of nodes allowed by memory requirements. In fact, the solution of the neutral system, Eq. (74), requires the allocation of matrices of sizes [N R × N Z + 2(N R + N Z )] 2 in each poloidal plane, thus being highly memory consuming 1 . As shown in Fig. 11 (a), the speed-up is almost ideal up to n pϕ = 128. The weak scalability is performed on neutral grids of size N R × N Z × N ϕ = 300 × 600 × 2n pϕ . Fig. 11 (b) shows an almost perfect efficiency as the system size increases. and is tested starting from 32 nodes, which is the minimum number of nodes required by memory constraints. The speed-up is almost ideal up to n pϕ = 128, corresponding to a total of 4608 cores (128 nodes). The weak scalability is carried out on neutral grids of size N R × N Z × N ϕ = 300 × 600 × 2n pϕ , where n pϕ is the number of MPI tasks along ϕ.

Convergence properties
We test here the convergence of GBS results with respect to the plasma and neutral grid refinement. We consider the magnetic field given by the analytical flux function in Eq. (76) and we compare the results of three simulations with increasing plasma grid resolution and no coupling with neutrals, and of three simulations with increasing neutral grid resolution and same plasma grid resolution. The results are analyzed when the simulations reach a quasi-steady state, resulting from the balance between the plasma and heat sources and losses at the wall. We also study the convergence with respect to the frequency of the neutral computation.
The simulations considered for studying the convergence with respect to the plasma grid are characterized by a coarse (∆R = ∆Z = 3.75 ρ s0 and R 0 ∆ϕ = 65.4 ρ s0 ), a medium (∆R = ∆Z = 2.50 ρ s0 and R 0 ∆ϕ = 49.1 ρ s0 ) and a fine (∆R = ∆Z = 1.67 ρ s0 and R 0 ∆ϕ = 32.7 ρ s0 ) resolution grid. In order to investigate the convergence, we consider the time and toroidally averaged radial profiles at the outer midplane of n, T e , T i and φ . As shown in Fig. 12, the radial profiles clearly converge to the ones obtained from the simulation with the finest resolution. Assuming that ρ s0 is the characteristic spatial length along the direction perpendicular to the magnetic field for phenomena occurring in our simulations, we conclude that ∆R = ∆Z 2.50 ρ s0 is the minimum spatial resolution of the plasma grid on the poloidal plane that guarantees convergence of simulation results. In order to guarantee stability of the simulations, the resolution along ϕ is chosen according to the ratio of the poloidal to the toroidal component of the magnetic field, The convergence with respect to the neutral grid refinement is tested by comparing the results of three simulations with the same plasma resolution, ∆R = ∆Z = 1.67 ρ s0 and R 0 ∆ϕ = 32. . The resolution along ϕ is given by the plasma grid. The neutral density is evaluated every ∆t = 0.08 R 0 /c s0 . The time traces of the spatially averaged neutral density, neutral temperature and neutral parallel velocity, v n = v n · b, show a clear convergence to the results of the finest neutral grid (see Fig. 13). We also analyze the time and the toroidal averages of the neutral density and ionization source, S iz , at the bottom wall, where the two strike points are located. Fig. 14 shows that there is no noticeable difference in the neutral density and ionization source profiles between simulations with different neutral grid resolutions, pointing out that the resolution ∆R = ∆Z 0.075 λ n is sufficient to guarantee the convergence of the simulation results. Finally, we study the convergence of GBS results with respect to the frequency of the neutral r/a calculation. We consider three simulations with a plasma grid resolution of ∆R = ∆Z = 1.67 ρ s0 and R 0 ∆ϕ = 32.7 ρ s0 , and a neutral grid resolution of ∆R = ∆Z = 0.075 λ n . The neutral density is evaluated every ∆t = 0.04 R 0 /c s0 , ∆t = 0.08 R 0 /c s0 and ∆t = 0.16 R 0 /c s0 , respectively. Fig. 15 shows that the time and toroidally averaged radial density profile at the outer midplane and neutral density at the bottom wall do not display any significant difference among the three simulations considered here. We conclude that evaluating the neutral density every ∆t 0.1 R 0 /c s0 guarantees the correct convergence of the results. We note that, since in our simulations turbulence is mainly driven by ballooning modes with a maximum growth rate of γ B = √ 2c s / L p R 0 [57], with L p = |p e /∂ R p e | the equilibrium pressure gradient length at the separatrix on the outer midplane, the minimum neutral calculation frequency that guarantees the convergence of the simulation results is γ B ∆t 0.7. Since the neutral and the plasma models can be evolved simultaneously, the number of MPI tasks for the plasma and neutral modules can be chosen to guarantee this neutral calculation frequency.

First simulation results in an experimentally-relevant scenario
As an example of application of the GBS version presented here, we report on the results of the first GBS electromagnetic simulation that includes the coupling with neutral dynamics of a lower single-null discharge performed in the TCV tokamak (tokamak major radius R 0 = 0.9 m and minor radius measured from the tokamak magnetic axis to the separatrix at the outer midplane a = 0.20 m). The magnetic equilibrium considered in the simulation is given by the equilibrium reconstruction of the TCV discharge #65402 at 1 s (plasma elongation and triangularity at the separatrix κ = 1.71 and δ = 0.35, respectively, toroidal magnetic field at the tokamak magnetic axis B 0 = 0.9 T, plasma current I p =146 kA). The discharge considered here is in forward field (ion-∇B drift direction pointing toward the active X-point). The upstream density and electron temperature at the separatrix, taken as reference density and temperature in the simulation, are n 0 = 0.6 × 10 19 m −3 and T e0 = 35 eV, respectively. With these reference values, the ion sound Larmor radius is ρ s0 1 mm, the sound speed is c s0 4.1 × 10 4 m/s, and the reference time is t 0 = R 0 /c s0 0.02 ms. In order to reduce the computational cost of the simulation, we consider here a domain corresponding to half size of TCV, i.e. L R = 300 ρ s0 and L Z = 600 ρ s0 . The dimensionless simulation parameters are ρ −1 * = 450, τ = 1, η 0e = 3 × 10 −4 , η 0i = 1, χ e = 20, χ i = 1, D f = 7 for f = {n, T e , T i , Ω,U e , v i }, m i /m e = 3000, β e0 = 2 × 10 −6 , and ν 0 = 0.05. The amplitude of the electron temperature source is chosen so that T e /T e0 1 at the separatrix. This leads to a power source with an intensity, in physical units, of approximately 150 kW, which is close to the experimental estimated value of the power crossing the separatrix, P sep = 120 kW. No external ion temperature source is used in the simulation and the ion heating is provided by the equipartition term in Eq. (6). The value of neutral reflection coefficient is α refl = 0.2. Although we consider a discrete gas puff in the private flux region, most of the neutral particles are generated from the ion recycling at the wall. A discrete pump at the bottom wall in the region R > 500ρ s0 is also considered.
The plasma spatial grid is N R × N Z × N ϕ = 150 × 300 × 64 and the neutral spatial grid is N R × N Z × N ϕ = 50 × 100 × 64. The time step is dt = 10 −5 R 0 /c s0 and the neutral kinetic equation is solved every ∆t = 0.04 R 0 /c s0 . The simulation reaches the turbulent quasi-steady state after approximately 80 R 0 /c s0 , starting from flat initial profiles of density and temperature.
We briefly discuss here the simulation results, showing that they are in agreement with typical experimental observations. Typical snapshots on a poloidal plane of n, T e , T i , φ , j and ψ are presented in Fig. 16. As experimentally observed (see, e.g., Refs. [58,59]), the tokamak core is characterized by low-amplitude density and temperature fluctuations, which increase from the core to the edge. A wave-like turbulent dynamics, experimentally observed [60,61], is clearly visible in the edge (see Fig. 16). When crossing the separatrix, the turbulent eddies experience a strong E × B shear and detach from the main plasma, forming filaments that propagate radially in the SOL, in agreement with experimental observations described, e.g., in Refs. [60,61]. In addition, as it is also experimentally observed [62,63,64,65], turbulence in the far SOL is characterized by intermittent events due to coherent plasma filaments [66].
The strong electric field gradient observed across the separatrix also reflects experimental observations [67,68]. Because of the ambipolarity of the plasma flow at the sheath, the electrostatic potential in the SOL is positive and proportional to the electron temperature [69]. Therefore, the electrostatic potential increases from the far SOL to the separatrix. On the other hand, the electrostatic potential inside the separatrix is negative with an electric field proportional to the ion pressure gradient. As a consequence, a strong electric field gradient and poloidal E × B shear form in the region across the separatrix. The E × B shear can play an important role in suppressing turbulence at the tokamak edge, as discussed in detail in Ref. [38], where GBS simulations are used to investigate the turbulent transport regimes in this region. Here we show that the E × B shear persists also when the neutral dynamics is included.
Typical snapshots of n n , S iz and v n are shown in Fig. 17. As expected, the neutral density is larger in the proximity of the target plates, where most of the recycling takes place, although a non-negligible fraction of neutrals can be found in the core, mainly in the region close to the Xpoint. In fact, a significant fraction of neutral particles reach the separatrix from the target on the high-field side before being ionized (see Fig. 17 (b)). In our simulation, core ionization accounts for approximately 25% of the total ionization source, which is in agreement with the fact that the TCV discharge considered here is in the low-density and low-recycling regime. We note that v n is positive in the outer and negative in the inner divertor regions (see Fig. 17 (c)). Since the direction of the magnetic field points towards (away from) the outer (inner) target, the parallel neutral flux is directed towards the wall at both targets.

Summary and conclusions
The present paper describes a new version of GBS, a three-dimensional, flux-driven, global, two-fluid turbulence code, developed for the self-consistent simulation of plasma turbulence and kinetic neutral dynamics in the tokamak boundary. In the version of GBS presented here, the simulation domain is extended to encompass the whole plasma volume, avoiding an artificial boundary with the tokamak core and retaining the core-edge-SOL turbulence interplay. Both the plasma and neutral implementations are carefully optimized leading to a significant speed-up of the code. In particular, a new iterative solver based on the PETSc library is implemented and optimized for the solution of the Poisson and Ampère equations, allowing us to efficiently carry out electromagnetic simulations at TCV size, while avoiding the use of the Boussinesq approximation. Also the neutral module is refactored and optimized by implementing an iterative solver based on the PETSc library.
The implementation of the plasma and neutral models is then carefully verified by means of the MMS, including, for the first time, the verification of the electromagnetic terms and of the kinetic neutral model. The verification of the neutral model is completed by a set of unit tests to verify the routines used to compute the integrals over the neutral trajectories.
Our tests show the efficient scalability of GBS on parallel high-performance computers, mainly thanks to the numerical scheme used to discretize the differential operators. The uniform Cartesian grid allows for a massive parallelisation through MPI domain decomposition. The splitting of the plasma and neutral communicators improves the parallelisation of GBS.
The convergence properties with respect to the plasma and neutral grid refinement are tested in typical turbulence simulations, showing that convergence of the simulation results is achieved with a plasma and neutral grid spacing of approximately 2.5 ρ s0 and 0.075 λ n , respectively. The convergence with respect to the neutral calculation frequency is also studied, showing that evaluating the neutral density every ∆t = 0.1R 0 /c s0 is sufficient to guarantee the convergence of the simulation results.
Finally, the results of the first GBS electromagnetic simulation of a TCV lower single-null discharge, which include the self-consistent evolution of neutral dynamics, are presented. A preliminary analysis of the simulation results show that GBS results agree well with experimental observations.
In the future, we plan to address three aspects to extend the flexibility and applicability of GBS. First, GBS will be adapted to account for the complex geometry of the tokamak wall, which may include baffles, such as those installed in TCV [70], or a complex divertor region. Second, by leveraging the relatively simple numerical scheme used in the version of GBS presented here, GBS will be ported to GPU, allowing for efficient simulations of large scale magnetic fusion devices, such as ITER and DEMO. Third, GBS plasma model equations will be extended by adding further moments of the electron and ion distribution functions [71], allowing us to retain kinetic effects that are expected to play an important role in the hot boundary of large tokamaks. arbitrary collisionality. Journal of Plasma Physics, 83(6), 2017.