2.1. The physical model

The multi-species neutral model is independent of the magnetic geometry and it does not require any modification. Although included in this code version, the coupling to neutrals is neglected here for simplicity. The plasma model is similar to the one described in Giacomin et al. (Reference Giacomin, Ricci, Coroado, Fourestey, Galassi, Lanti, Mancini, Richart, Stenger and Varini2022b ), with a few modifications that are described in this section. The plasma model equations implemented in the GBS version used herein are (assuming a single ion species with $Z_i=1$ )

(2.1) \begin{align} \frac {\partial n}{\partial t} &= -\frac {1}{B}[\phi ,n]+\frac {2}{eB} [\mathcal{C}(p_e)-en\mathcal{C}(\phi )] -\boldsymbol{\nabla} _{\parallel }(n v_{\parallel e}) + D_n{\nabla} _{\perp }^2 n\nonumber \\ &\quad +s_n+\nu _{\text{iz}}n_{\text{n}}-\nu _{\text{rec}}n,\\[-10pt]\nonumber \end{align}

(2.2) \begin{align} \frac {\partial \varOmega }{\partial t} & = -\frac {1}{B}\boldsymbol{\nabla }\boldsymbol{\cdot }\left(\frac {1}{\varOmega _{ci}}[\phi ,\boldsymbol{\omega }]\right) - \boldsymbol{\nabla }\boldsymbol{\cdot }\left( \frac {v_{\parallel i}}{\varOmega _{ci}}\boldsymbol{\nabla} _\parallel \boldsymbol{\omega }\right)+ \frac {1}{e}\boldsymbol{\nabla} _{\parallel }j_{\parallel }\nonumber \\ & \quad + \frac {2}{eB} \mathcal{C}(p_e + p_i) + D_{\varOmega }{\nabla} _\perp ^2 \varOmega -\frac {n_{\text{n}}}{n}\nu _{\text{cx}}\varOmega , \end{align}

(2.3) \begin{align} \frac {\partial U_{\parallel e}}{\partial t} &= -\frac {1}{B}[\phi ,v_{\parallel e}] + \frac {e}{m_e}\left( \frac {j_\parallel }{\sigma _\parallel }+\boldsymbol{\nabla} _\parallel \phi -\frac {1}{en}\boldsymbol{\nabla} _\parallel p_e-\frac {0.71}{e}\boldsymbol{\nabla} _\parallel T_e -\frac {4\eta _{0e}}{3en}{\nabla} _\parallel ^2 v_{\parallel e}\right) \nonumber \\[5pt] & \quad - v_{\parallel e}\boldsymbol{\nabla} _\parallel v_{\parallel e} + D_{v_{\parallel e}}{\nabla} _\perp ^2 v_{\parallel e}+\frac {n_{\text{n}}}{n}(\nu _{\text{en}}+2\nu _{\text{iz}})(v_{\parallel n}-v_{\parallel e}),\\[-7pt]\nonumber \end{align}

(2.4) \begin{align} \frac {\partial v_{\parallel i}}{\partial t} &= -\frac {1}{B}[\phi ,v_{\parallel i}] - v_{\parallel i}\boldsymbol{\nabla} _\parallel v_{\parallel i} - \frac {1}{m_i n}\boldsymbol{\nabla} _\parallel (p_e+ p_i) -\frac {4\eta _{0i}}{3m_in}{\nabla} _\parallel ^2 v_{\parallel i} + D_{v_{\parallel i}}{\nabla} _\perp ^2 v_{\parallel i}\nonumber \\[5pt] & \quad +\frac {n_{\text{n}}}{n}(\nu _{\text{iz}}+\nu _{\text{cx}})(v_{\parallel n}-v_{\parallel i}),\\[-7pt]\nonumber \end{align}

(2.5) \begin{align} \frac {\partial T_e}{\partial t} &= -\frac {1}{B}[\phi ,T_e] - v_{\parallel e}\boldsymbol{\nabla} _\parallel T_e + \frac {2}{3}T_e\left[0.71\frac {\boldsymbol{\nabla} _\parallel j_\parallel }{en} - \boldsymbol{\nabla} _\parallel v_{\parallel e}\right] + \boldsymbol{\nabla} _\parallel (\chi _{\parallel e}\boldsymbol{\nabla} _\parallel T_e) + D_{T_e}{\nabla} _\perp ^2 T_e \nonumber \\[5pt] &\quad + \frac {4}{3}\frac {T_e}{eB}\left[\frac {7}{2}\mathcal{C}(T_e)+\frac {T_e}{n}\mathcal{C}(n)-e\mathcal{C}(\phi )\right] + s_{T_e}-\frac {n_{\text{n}}}{n}\nu _{\text{en}}m_e\frac {2}{3}v_{\parallel e}(v_{\parallel n}-v_{\parallel e})\nonumber \\[5pt] &\quad - \frac {4}{3}\frac {m_e}{m_i}\frac {1}{\tau _e} (T_e- T_i) + \frac {n_{\text{n}}}{n}\nu _{\text{iz}}\left[-\frac {2}{3}E_{\text{iz}}-T_e+m_e v_{\parallel e}\Bigl (v_{\parallel e}-\frac {4}{3}v_{\parallel n}\Bigr )\right] ,\\[-7pt]\nonumber \end{align}

(2.6) \begin{align} \frac {\partial T_i}{\partial t} &= -\frac {1}{B}[\phi ,T_i] - v_{\parallel i}\boldsymbol{\nabla} _\parallel T_i + \frac {4}{3}\frac {T_i}{eB}\left[\mathcal{C}\left(T_e +\frac {T_e}{n}\mathcal{C}(n)-e\mathcal{C}(\phi )\right.\right] - \frac {10}{3}\frac {T_i}{eB}\mathcal{C}(T_i) \nonumber \\[5pt] &\quad + \frac {2}{3}T_i\left[(v_{\parallel i}-v_{\parallel e})\frac {\boldsymbol{\nabla} _\parallel n}{n} -\boldsymbol{\nabla} _\parallel v_{\parallel e}\right] + \boldsymbol{\nabla} _\parallel (\chi _{\parallel i}\boldsymbol{\nabla} _\parallel T_i) + D_{T_i}{\nabla} _\perp ^2 T_i + s_{T_i} \nonumber \\[5pt] &\quad + \frac {4}{3} \frac {m_e}{m_i}\frac {1}{\tau _e} (T_e- T_i) + \frac {n_{\text{n}}}{n}(\nu _{\text{iz}}+\nu _{\text{cx}})\left[T_n-T_i+\frac {1}{3}(v_{\parallel n}-v_{\parallel i})^2\right], \end{align}

which are coupled to the Poisson and Ampère equations

(2.7) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\Bigl ( \frac {n}{B\varOmega _{ci}} \boldsymbol{\nabla} _\perp \phi \Bigr ) & = \varOmega - \boldsymbol{\nabla }\boldsymbol{\cdot }\left(\frac {1}{m_i\varOmega _{ci}^2}\boldsymbol{\nabla} _\perp p_i\right),\\[-8pt]\nonumber \end{align}

(2.8) \begin{align} \left( {\nabla} _\perp ^2 - \frac {e^2\mu _0}{m_e}n\right)v_{\parallel e} & = {\nabla} _\perp ^2 U_{\parallel e} - \frac {e^2\mu _0}{m_e}n v_{\parallel i}. \end{align}

In (2.1)–(2.8), $n$ is the plasma density ( $n=n_e=n_i$ for a single ion species with $Z_i=1$ ), $T_e$ and $T_i$ are the electron and ion temperatures, $p_e=nT_e$ and $p_i=nT_i$ are the electron and ion pressures, $U_{\parallel e} = v_{\parallel e} + e \psi /m_e$ is the sum of the electron inertia and the electromagnetic induction contributions, $v_{\parallel i}$ is the ion parallel velocity, $\varOmega = \boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{\omega }/\varOmega _{ci}) = \boldsymbol{\nabla }\boldsymbol{\cdot }\bigl [(n \boldsymbol{\nabla} _\perp \phi + \boldsymbol{\nabla} _\perp p_i/e)/(\varOmega _{ci}B)\bigr ]$ is the scalar vorticity, $\varOmega _{ci}=eB/m_i$ is the ion cyclotron frequency, $e$ is the elementary charge, $m_i$ and $m_e$ $\nu _{\textrm {iz}}$ , $\nu _{\textrm {rec}}$ , $\nu _{\textrm {cx}}$ and $\nu _{\textrm {en}}$ are the ionisation, recombination, charge-exchange and electron–neutral collisions frequencies, respectively (see Giacomin et al. Reference Giacomin, Ricci, Coroado, Fourestey, Galassi, Lanti, Mancini, Richart, Stenger and Varini2022b for details on the plasma–neutral coupling), $E_{\textrm {iz}}$ is the ionisation energy, $\sigma _\parallel = 1.96 n e^2 \tau _e/m_e$ is the parallel electric conductivity and $\chi _e$ and $\chi _i$ are the electron and ion parallel thermal conductivities, as defined in the Braginskii model (Braginskii Reference Braginskii1965).

The spatial operators in (2.1)–(2.8) are the $\boldsymbol{E}\times \boldsymbol{B}$ convective term (Poisson brackets)

(2.9) \begin{align} [\phi ,f]=\boldsymbol{b}\ \boldsymbol{\cdot }\ \bigl (\boldsymbol{\nabla }\phi \times \boldsymbol{\nabla }f\bigr ), \end{align}

the curvature operator

(2.10) \begin{align} \mathcal{C}(f)=\frac {B}{2}\left(\boldsymbol{\nabla }\times \frac {\boldsymbol{b}}{B}\right)\boldsymbol{\cdot }\boldsymbol{\nabla }f, \end{align}

the parallel gradient, which includes the electromagnetic flutter contribution

(2.11) \begin{align} \boldsymbol{\nabla} _\parallel f=\boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{\nabla }f + \frac {1}{B}[\psi ,f], \end{align}

and the perpendicular Laplacian

(2.12) \begin{align} {\nabla} _\perp ^2 f=\boldsymbol{\nabla }\boldsymbol{\cdot }\bigl [(\boldsymbol{b}\times \boldsymbol{\nabla }f)\times \boldsymbol{b}\bigl ], \end{align}

where $f$ is a general scalar function, $\mathbf{\textit{E}}=-\boldsymbol{\nabla}\phi$ is the electric field, and $\mathbf{\textit{b}} = \mathbf{\textit{B}}/B$ is the unit vector along the magnetic field direction.

Equations (2.1)–(2.8) are implemented in GBS in their dimensionless form, where $n$ is normalised to the reference density $n_0$ , $T_e$ and $T_i$ are normalised to the reference electron temperature $T_{e0}$ , $v_{\parallel e}$ and $v_{\parallel i}$ are normalised to the reference sound speed $c_{s0}=\sqrt {T_{e0}/m_i}$ , the magnetic field is normalised to the reference magnetic field strength $B_0$ , the electrostatic potential is normalised to $T_{e0}/e$ , $\psi$ is normalised to $\rho _{s0}B_0$ , with $\rho _{s0}=c_{s0}/\varOmega _{ci0}$ , perpendicular and parallel lengths are normalised to $\rho _{s0}$ and $R_0$ , respectively, with $R_0$ the machine major radius, and the time is normalised to $R_0/c_{s0}$ . The reference electron plasma $\beta$ , $\beta _{e0}=2\mu _0n_0T_{e0}/B_0^2$ , and the normalised Spitzer resistivity, $\nu =e^2n_0R_0/(m_ic_{s0}\sigma _\parallel ) = \nu _0T_e^{-3/2}$ , stem from the normalisation process, with

(2.13) \begin{align} \nu _0 = \frac {4\sqrt {2\pi }}{5.88}\frac {e^4}{(4\pi \epsilon _0)^2}\frac {\sqrt {m_e}R_0n_0\lambda }{m_ic_{s0}T_{e0}^{3/2}} \end{align}

the reference resistivity, which depends on the reference density and reference electron temperature ( $\lambda$ is the Coulomb logarithm).

The main difference between the physical model considered in the present work and the one described in Giacomin et al. (Reference Giacomin, Ricci, Coroado, Fourestey, Galassi, Lanti, Mancini, Richart, Stenger and Varini2022b ) is related to the spatial dependence of the modulus of the magnetic field, $B$ . In fact, under the assumption of weak poloidal magnetic, the magnetic field modulus is approximated to $B^2 = B_\varphi ^2 + B_p^2 = B_\varphi ^2(1+B_p^2/B_\varphi ^2)\simeq B_\varphi ^2$ , namely $B\simeq |B_\varphi |$ . In the large-aspect-ratio limit and by considering $B_0$ as the magnetic field strength at the magnetic axis, the modulus of the magnetic field is approximated by $|B_\varphi | = B_0R_0/R\simeq B_0$ . This leads to $B\simeq B_0$ and $\varOmega _{ci} \simeq \varOmega _{ci0}$ . Consequently, the magnetic field dependence in the scalar vorticity is neglected in Giacomin et al. (Reference Giacomin, Ricci, Coroado, Fourestey, Galassi, Lanti, Mancini, Richart, Stenger and Varini2022b ), i.e. $\varOmega \propto \boldsymbol{\nabla }\boldsymbol{\cdot }(n \boldsymbol{\nabla} _\perp \phi + \boldsymbol{\nabla} _\perp p_i/e)$ . This approximation simplifies the implementation of the vorticity equation (see (2.2)) and Poisson’s law (see (2.7)). In addition, the factor $1/B$ in front of the Poisson brackets in (2.1)–(2.6) is neglected in Giacomin et al. (Reference Giacomin, Ricci, Coroado, Fourestey, Galassi, Lanti, Mancini, Richart, Stenger and Varini2022b ) when these equations are written in dimensionless form ( $B/B_0=1$ ). These approximations are justified in conventional tokamaks, where $B$ varies weakly within the plasma volume. On the other hand, they do not hold in general, and the spatial dependence of $B$ is therefore retained in the present work.

2.2. The extended GBS operators

The GBS domain is discretised on a cylindrical uniformly spaced $(R, \varphi , Z)$ grid. The spatial derivatives are implemented with a fourth-order centred finite-difference scheme and the time evolution is provided by a fourth-order Runge–Kutta algorithm. An iterative solver based on data management for structured grids (DMDA) of the PETSc library (Balay et al. Reference Balay2019) is used for the implementation of the Poisson and Ampère equations. While the spatial discretisation and temporal advancement schemes are identical to those employed in the main GBS version, the differential operators in (2.9)–(2.12) are here derived in their most general form, explicitly avoiding the asymptotic expansions in $\epsilon \equiv r/R_0 \ll 1$ and $\delta \equiv B_p/B \ll 1$ .

The magnetic field is written in cylindrical coordinates, $\boldsymbol{B} = B_R {\boldsymbol{e}_R} + B_Z \boldsymbol{e}_Z + B_\varphi \boldsymbol{e}_\varphi$ , where $\boldsymbol{e}_R$ , $\boldsymbol{e}_Z$ and $\boldsymbol{e}_\varphi$ denote the radial, vertical and toroidal directions, respectively. Whereas the differential operators in Giacomin et al. (Reference Giacomin, Ricci, Coroado, Fourestey, Galassi, Lanti, Mancini, Richart, Stenger and Varini2022b ) are expressed in terms of the poloidal flux function $\varPsi$ and its spatial derivatives, in the present work they are instead formulated directly in terms of the magnetic field components $B_R$ , $B_Z$ , $B_\varphi$ and their corresponding spatial derivatives. This choice enables turbulence simulations also in magnetic configurations without well-converted magnetic flux surfaces, such as those with magnetic chaos.

In cylindrical coordinates, the Poisson brackets operator, (2.9), becomes

(2.14) \begin{align} [\phi , f] &= \frac {1}{R}\frac {B_R}{B}(\partial _\varphi \phi {\partial _Z} f - {\partial _Z}\phi \partial _\varphi f)\nonumber\\[7pt] & \quad +\frac {1}{R}\frac {B_Z}{B}(\partial _R\phi \partial _\varphi f- \partial _\varphi \phi \partial _R f)\nonumber\\[7pt] & \quad +\frac {B_\varphi }{B}({\partial _Z}\phi \partial _R f- \partial _R\phi {\partial _Z} f), \end{align}

which is normalised to $1/\rho _{s0}^2$ . Only the last term is retained in Giacomin et al. (Reference Giacomin, Ricci, Coroado, Fourestey, Galassi, Lanti, Mancini, Richart, Stenger and Varini2022b ) when the operator is expanded in $\epsilon$ and $\delta$ . As in Giacomin et al. (Reference Giacomin, Ricci, Coroado, Fourestey, Galassi, Lanti, Mancini, Richart, Stenger and Varini2022b ), the Poisson brackets operator is discretised using a fourth-order Arakawa scheme (Peterson & Hammett Reference Peterson and Hammett2013).

The curvature operator, (2.10), is written as

(2.15) \begin{align} \mathcal{C}(f) = \frac {1}{2 B^3 R}&\left(\partial _\varphi f \left(\partial _Z B_R \big(B^2-2 B_R ^2\big)-\partial _R B_Z \big(B^2-2 B_Z ^2\big) \right.\right.\nonumber\\[3pt] &\left. -2 B_R \bigl (B_Z (\partial _Z B_Z -\partial _R B_R )+B_\varphi \partial _Z B_\varphi \bigr )+2 B_Z B_\varphi \partial _R B_\varphi \right)\nonumber \\[3pt] &+\partial _Z f \bigl (-\partial _\varphi B_R (B^2-2 B_R ^2)+B_\varphi (B^2+2 B_R (\partial _\varphi B_\varphi -R\partial _R B_R )\nonumber\\[3pt] &-2 R B_Z \partial _R B_Z) + RB^2 \partial _R B_\varphi +2 B_R B_Z \partial _\varphi B_Z -2 R B_\varphi ^2 \partial _R B_\varphi \bigr )\nonumber \\[3pt] &+\partial _R f \bigl (\partial _\varphi B_Z \big(B^2-2 B_Z ^2\big)-R\partial _Z B_\varphi (B^2-2 B_\varphi ^2)\nonumber \\[3pt] &\left. +B_R (2 R B_\varphi \partial _Z B_R -2 B_Z \partial _\varphi B_R )-2 B_Z B_\varphi (\partial _\varphi B_\varphi -R\partial _Z B_Z )\bigr )\right), \end{align}

and it is normalised to $1/(R_0\rho _{s0})$ . If only the leading-order terms in $\epsilon$ and $\delta$ are retained and the local current is neglected, as in Giacomin et al. (Reference Giacomin, Ricci, Coroado, Fourestey, Galassi, Lanti, Mancini, Richart, Stenger and Varini2022b ), the normalised curvature operator is drastically simplified to $\mathcal{C}(f) = (B_\varphi /B_0)\partial _Z f$ , which corresponds to the curvature associated with the toroidal magnetic field.

The parallel gradient operator, (2.11), becomes

(2.16) \begin{align} \boldsymbol{\nabla} _\parallel f = \frac {B_R}{B}\partial _R f + \frac {B_Z}{B}\partial _Z f + \frac {B_\varphi }{BR}\partial _\varphi f, \end{align}

and it is normalised to $1/R_0$ . The parallel Laplacian operator is obtained by applying twice the $\boldsymbol{\nabla} _\parallel$ operator, leading to

(2.17) \begin{align} {\nabla} _\parallel ^2 f = \,& a_R \partial _R f + a_Z \partial _Z f + a_\varphi \partial _\varphi f + a_{RR} \partial _R^2f +a_{ZZ}\partial _Z^2f+a_{\varphi \varphi }\partial _\varphi ^2 f \nonumber\\&+ a_{RZ}\partial _{RZ}f + a_{R\varphi }\partial _{R\varphi }f+ a_{Z\varphi }\partial _{Z\varphi }f, \end{align}

where the coefficients $a_R$ , $a_Z$ , $a_\varphi$ , $a_{RR}$ , $a_{ZZ}$ , $a_{\varphi \varphi }$ , $a_{RZ}$ , $a_{R\varphi }$ , $a_{Z\varphi }$ are given in Appendix A.

The perpendicular Laplacian, (2.12), is written as

(2.18) \begin{align} {\nabla} _\perp ^2 f & = c_R \partial _R f + c_Z \partial _Z f + c_\varphi \partial _\varphi f + c_{RR} \partial _R^2f+c_{ZZ}\partial _Z^2f+c_{\varphi \varphi }\partial _\varphi ^2 f \nonumber\\& \quad + c_{RZ}\partial _{RZ}f + c_{R\varphi }\partial _{R\varphi }f+ c_{Z\varphi }\partial _{Z\varphi }f, \end{align}

where the coefficients $c_R$ , $c_Z$ , $c_\varphi$ , $c_{RR}$ , $c_{ZZ}$ , $c_{\varphi \varphi }$ , $c_{RZ}$ , $c_{R\varphi }$ , $c_{Z\varphi }$ are given in Appendix A. The perpendicular Laplacian is normalised to $1/\rho _{s0}^2$ . The Poisson equation, (2.7), involves a modified perpendicular Laplacian operator of the form $\boldsymbol{\nabla }\boldsymbol{\cdot }(f \boldsymbol{\nabla} _\perp \phi )$ , with $f=n/B^2$ on the left-hand side and $f=1/B^2$ on the right-hand side. This operator is written as in (2.18) with the coefficients given in Appendix A.

At leading order in $\epsilon$ and $\delta$ , the perpendicular Laplacian in (2.18) undergoes a substantial simplification, such that $c_{RR}$ and $c_{ZZ}$ are the only coefficients that remain non-vanishing. Consequently, the Laplacian operator appearing in the Poisson and Ampère equations is reduced to a two-dimensional operator acting on the $(R, Z)$ plane. However, in the most general case (including RFP plasmas), the coefficients in (2.18) that multiply the derivatives in the toroidal direction do not vanish, thus requiring the use of three-dimensional solvers in (2.7)–(2.8). This is the major code modification with respect to Giacomin et al. (Reference Giacomin, Ricci, Coroado, Fourestey, Galassi, Lanti, Mancini, Richart, Stenger and Varini2022b ). The three-dimensional iterative solvers are implemented using the DMDA of the PETSc library. As in the two-dimensional case, the best performance is achieved with the algebraic multigrid preconditioner boomerAMG provided by the HYPRE library (Falgout & Yang Reference Falgout and Yang2002) and the flexible generalised minimal residual (flexible GMRES) solver (Saad Reference Saad1993). The preconditioner parameters are adapted to the three-dimensional solver, changing the coarsen type from Falgout to HMIS and increasing the strong threshold from 0.25 to 0.7.

2.3. Verification and scalability

The modifications applied to the differential operators and, more critically, the implementation of three-dimensional solvers, necessitate a renewed verification of the code. This verification is performed using the method of manufactured solutions (MMS) (Riva et al. Reference Riva, Ricci, Halpern, Jolliet, Loizu and Mosetto2014), as described below.

The MMS is a numerical technique that consists of verifying the correct implementation of a model $M$ . The verification is successful if the numerical solution $s_h$ of the discretised model $M_h$ converges to the analytical solution $s$ of the model $M$ as the discretisation error $h$ tends to zero, i.e. $e_h = \|s-s_h\| \to h^p$ as $h \to 0$ with $p$ the order of the discretisation scheme ( $p=4$ in GBS since both space and time are discretised with a fourth-order method). Since the analytical solution of the model equations is in general unknown, the fundamental idea underlying the MMS is to prescribe an arbitrary analytical function, denoted by $u$ , and to evaluate the corresponding source term $S = M(u)$ , where $M$ represents the continuous model operator. The manufactured function $u$ is then, by construction, an exact solution of the modified model $N$ , defined as $N(u) = M(u) - S = 0$ . Since $S$ can be analytically evaluated, the discrete operators $N_h$ and $M_h$ are affected by the same discretisation error. Consequently, the verification of the correctness of the discrete formulation $N_h$ is equivalent to the verification of the discrete model operator $M_h$ .

The verification is carried out by choosing the same numerical set-up as in Giacomin et al. (Reference Giacomin, Ricci, Coroado, Fourestey, Galassi, Lanti, Mancini, Richart, Stenger and Varini2022b ), i.e. $\rho _*^{-1}= R_0/\rho _{s0} = 100$ , $\nu _0=1.0$ , $m_e/m_i=1.0$ , $\beta _{e0} = 10^{-4}$ , and domain sizes in the radial and vertical directions of $L_R=37.5\,\rho _{s0}$ and $L_Z=50\,\rho _{s0}$ . The magnetic field is defined as

(2.19) \begin{align} B_R &= \sin \biggl (\frac {2\pi R}{R_{\mathrm{max}}-R_{\mathrm{min}}}\biggr )\sin \left(\frac {2\pi Z}{Z_{\mathrm{max}}-Z_{\mathrm{min}}}\right)\cos (\varphi ) ,\\[-7pt]\nonumber\end{align}

(2.20) \begin{align} B_Z &= \cos \biggl (\frac {2\pi R}{R_{\mathrm{max}}-R_{\mathrm{min}}}\biggr )\cos \left(\frac {2\pi Z}{Z_{\mathrm{max}}-Z_{\mathrm{min}}}\right)\cos (\varphi ) ,\\[-7pt]\nonumber\end{align}

(2.21) \begin{align} B_\varphi &= 2 + \cos \biggl (\frac {2\pi R}{R_{\mathrm{max}}-R_{\mathrm{min}}}\Bigr )\cos \left(\frac {2\pi Z}{Z_{\mathrm{max}}-Z_{\mathrm{min}}}\right)\sin (\varphi ) ,\end{align}

in such a way that all coefficients of the differential operators are non-zero. We remark that the verification is independent of the choice of the equilibrium magnetic field, as long as it leads to non-trivial differential operators. Similarly to Giacomin et al. (Reference Giacomin, Ricci, Coroado, Fourestey, Galassi, Lanti, Mancini, Richart, Stenger and Varini2022b ), the prescribed solutions $u=n, T_e, T_i, \phi , \varOmega , v_{\parallel e}, v_{\parallel i}, \psi$ are chosen as

(2.22) \begin{align} u (R, \varphi , Z) = A_u\bigl [B_u + \sin (C_uZ+\alpha _u)\sin (D_u\varphi +\beta _u)\sin (E_u t + F_u R + \gamma _u) \bigr ] ,\end{align}

where $A_u$ , $B_u$ , $C_u$ , $D_u$ , $E_u$ , $F_u$ , $\alpha _u$ , $\beta _u$ and $\gamma _u$ are arbitrary constants. The source term is computed by substituting the prescribed solution into (2.1)–(2.7), with the magnetic field given by (2.19)–(2.21).

The verification is performed by considering five grids of increasing spacing, with the finest grid $(N_R^*=128, N_\varphi ^*=128, N_Z^* = 128)$ and ${\rm d}t^*=1.25\times 10^{-5}\,R_0/c_{s0}$ . The discretisation parameter $h\in \{1, 2, 4, 8\}$ is defined as $h=N_R^*/N_R=N_Z^*/N_Z=N_\varphi ^*/N_\varphi = {\rm d}t/{\rm d}t^*$ . Figure 1 shows the $\mathcal{L}_2$ and the $\mathcal{L}_\infty$ norms of the discretisation error as a function of $h$ , along with the order of convergence $p=\log (e_{rh}/e_h)/\log {r}$ , where $r=2$ is the grid refinement. The discretisation error scales as $h^4$ , hence verifying the implementation of the modified differential operators and three-dimensional solvers.

Figure 1.

Results of the verification test. The $\mathcal{L}_\infty$ (a) and $\mathcal{L}_2$ (b) norms of the discretisation error as functions of $h$ . Order of convergence $p$ for the $\mathcal{L}_\infty$ (c) and the $\mathcal{L}_2$ (d) norms of the discretisation errors. The order of convergence is evaluated from two consecutive grid resolutions.

Figure 1. Long description

Panel A: A line graph shows the norms of the discretization error as functions of grid resolution (h). The y-axis represents the norms of the discretization error on a logarithmic scale, and the x-axis represents the grid resolution (h) on a linear scale. Different lines represent various variables such as phi, n, omega, v_parallel_i, T_e, T_i, v_parallel_e, and psi. Panel B: Another line graph shows the norms of the discretization error as functions of grid resolution (h). The y-axis represents the norms of the discretization error on a logarithmic scale, and the x-axis represents the grid resolution (h) on a linear scale. Different lines represent various variables such as phi, n, omega, v_parallel_i, T_e, T_i, v_parallel_e, and psi. Panel C: A line graph shows the order of convergence for the L_infinity norms of the discretization errors. The y-axis represents the order of convergence, and the x-axis represents the grid resolution (h) on a linear scale. Different lines represent various variables such as phi, n, omega, v_parallel_i, T_e, T_i, v_parallel_e, and psi. Panel D: A line graph shows the order of convergence for the L_2 norms of the discretization errors. The y-axis represents the order of convergence, and the x-axis represents the grid resolution (h) on a linear scale. Different lines represent various variables such as phi, n, omega, v_parallel_i, T_e, T_i, v_parallel_e, and psi.

The parallelisation in GBS is implemented through the domain decomposition by using the message passing interface (MPI) library, which, thanks to the use of a uniform cylindrical grid, offers a very efficient scalability up to a large number of MPI tasks. While the code parallelisation is not affected by the modifications on the differential operators, the implementation of new iterative solvers may impact the code performance. It is therefore important to test how the wall time of the solvers scales with the number of cores (strong scaling). The scalability test is carried on Discoverer HPC (2 $\times$ AMD EPYC 7H12 64-Core Processor) considering a grid size of $N_R\times N_\varphi \times N_Z = 160 \times 128 \times 160$ , which is a typical grid for a half-size RFX-mod RFP simulation (§ 3). Figure 2 shows the wall-clock time of one step, computed as time average over 100 steps, as a function of the number of nodes. The time taken by the three-dimensional Poisson and Ampère solvers is also indicated in figure 2. With the chosen grid, GBS scales almost ideally up to 32 nodes (4092 cores). Importantly, the newly implemented three-dimensional solvers scale well with the number of MPI tasks. For comparison, the time-to-solution related to Poisson and Ampère two-dimensional solvers is also shown in figure 2(a). The simulation with two-dimensional solvers is approximately a factor 1.5 faster than the same simulation with three-dimensional solvers, although this factor is sensitive to the computing architecture and MPI domain decomposition, and it can rise above 2.

Figure 2.

(a) Strong scalability test performed on Discoverer HPC (2 $\times$ AMD EPYC 7H12 64-Core Processor) considering a grid size of $N_R\times N_\varphi \times N_Z = 160 \times 128 \times 160$ and increasing the number of computing nodes from 1 to 32 (1 node = 128 cores). The wall-clock execution times of the Poisson (orange) and Amperè (green) solvers are shown along with the total wall-clock time required for a single time step (blue). The ideal scaling is denoted by the red markers, while the two-dimensional solver times are indicated by the black markers. (b) Time to solution for a single time step for different MPI domain decompositions in GBS simulations employing three-dimensional and two-dimensional solvers. The numbers of MPI tasks distributed along the $R$ , $Z$ and $\varphi$ directions are denoted by $n_{pR}$ , $n_{pZ}$ and $n_{p\varphi }$ , respectively.

Because of the intense parallel communication required by Poisson and Ampère solvers, the GBS version with two-dimensional solvers scales more efficiently in the toroidal direction. Indeed, the parallelisation strategy usually adopted is to distribute the MPI tasks within a computing node (or a fraction of it) on the poloidal plane, and increase the number of MPI tasks along $\varphi$ consequently. The implementation of three-dimensional solvers makes the choice of the MPI decomposition irrelevant, as the three directions are equivalent. Indeed, figure 2(b) shows that the time to solution weakly changes with the chosen MPI decomposition. A larger variation is instead observed if two-dimensional solvers are employed, with the most efficient case corresponding to the one with the largest number of MPI tasks along $\varphi$ .

Although the newly implemented three-dimensional iterative solvers do not degrade the excellent scalability properties of GBS, their heavier computational cost makes GBS simulations more expensive. Therefore, it is advisable to use this modified version of GBS only when the approximations $\delta \ll 1$ and/or $\epsilon \ll 1$ are violated or when it is not possible to define a flux function $\varPsi$ .

3.1. Simulation set-up

The magnetic field considered in this work is reconstructed with the NCT code (Zanca & Terranova Reference Zanca and Terranova2004) from the RFP plasma discharge no. 29187 at $t=209$ ms of the RFX-mod device (minor radius $a=0.46$ m and major radius $R_0=2.0$ m). This is the same magnetic field considered in Giacomin et al. (Reference Giacomin, Momo, Predebon, Vianello and Zuin2026) and it corresponds to a QSH state with plasma current $I_p\simeq 1.5$ MA and shallow reversal. The magnetic field reconstruction includes the axisymmetric equilibrium and all the MHD modes with poloidal mode numbers $m\in \{-1, 0, 1, 2\}$ and toroidal mode numbers up to $n=24$ . The dominant mode is $(m, n) = (1, 7)$ , which yields an approximately helical magnetic field configuration characterised by a toroidal periodicity of $n = 7$ (see Giacomin et al. Reference Giacomin, Momo, Predebon, Vianello and Zuin2026 for additional details on the magnetic field). The Poincaré map of the magnetic field on the $(R,\varphi)$ -midplane, depicted in figure 3, reveals pronounced magnetic stochasticity, which originates from the overlap of magnetic islands associated with adjacent resonant surfaces. The magnetic island near the wall ( $r/a\simeq 1$ ) is caused by the $m=0,\, n=7$ unstable tearing mode.

The number of grid points in the radial, toroidal and vertical directions is $N_R=160$ , $N_\varphi = 128$ and $N_Z=160$ , respectively. The size of the simulation domain is $L_R= 300\,\rho _{s0}$ in the radial direction, $L_Z=300\,\rho _{s0}$ in the vertical direction and $L_\varphi = 2\pi /7$ in the toroidal direction, with $\rho _*^{-1}= R_{s0}/\rho _{s0}=750$ . In particular, the quasi-periodicity associated with the dominant MHD mode with toroidal mode number $n=7$ is exploited in this work to restrict the computational toroidal domain to $0 \leqslant \varphi \lt 2\pi /7$ . By considering typical values of electron temperature and magnetic field in the boundary of RFP discharges, $T_{e0} \simeq 60$ eV and $B_0 \simeq 0.6$ T, the reference Larmor radius is $\rho _{s0} \simeq 1.5$ mm. This yields a poloidal cross-section of size $L_R=L_Z\simeq 0.45$ m. Consequently, the size of the poloidal plane in this set of GBS simulations corresponds to half of the physical size of the RFX-mod device. The major radius in the simulations is also half of the RFX-mod one. Adopting a reduced system size enables comprehensive parameter-space exploration that would otherwise be computationally prohibitive at full scale.

The set of GBS simulations considered here includes a scan on the normalised resistivity, $\nu _0\in \{0.05,\, 0.5,\, 5.0,\, 50\}$ at $\beta _{e0}=10^{-3}$ , and a scan on the reference electron $\beta$ , $\beta _{e0}\in \{0.001,\, 0.005,\, 0.01,\, 0.05\}$ at $\nu _0=0.05$ . We note that varying $\nu _0$ is effectively equivalent to changing the collisionality through the reference density $n_0$ (see (2.13)). The reference quantities are evaluated at the plasma boundary and the experimental values of $\nu _0$ and $\beta _{e0}$ in the considered RFX-mod plasma discharge are approximately $\nu _0 = 0.05$ and $\beta _{e0} = 0.001$ . We highlight that values of $\nu _0$ greater than 1 are unrealistic for a RFX-mod RFP plasma, even in a poorly confined single helicity regime. Similarly, $\beta _{e0}$ values exceeding 0.005 are also not realistic. In this work, GBS simulations employing such non-physical values of $\nu _0$ and $\beta _{e0}$ are carried out solely to test the theoretical framework presented in § 4.

These parameter scans are performed at fixed temperature and density sources, which are uniform functions along $\varphi$ defined as

(3.1) \begin{align} s_n &= s_{n0} \exp \left[-\frac {(r-r_{n0})^2}{\Delta r_n^2}\right],\\[-7pt]\nonumber \end{align}

(3.2) \begin{align} s_{T_e} &= \frac {s_{T0}}{2}\left[1-\tanh \left(\frac {r-r_{T0}}{\Delta r_T}\right)\right], \end{align}

with $s_{n0} = 0.01\,n_0c_{s0}/R_0$ , $r_{n0}=0.5\,a$ , $\Delta r_n = 0.1\,a$ , $s_{T0}=0.05\,T_{e0}c_{s0}/R_0$ , $r_{T0}=0.2\,a$ and $\Delta r_T=0.05\,a$ . These parameter values are chosen such that the total particle fuelling rate is approximately $10^{21}\,\textrm {particles}\,{\rm s}^{-1}$ and the applied heating power is $1\,\mathrm{MW}$ (the heating power corresponds to the power flowing to the boundary region). Since the latter is a function of density and temperature profiles, its value may vary across the simulation scan. No ion temperature source is considered. The functional form of the source terms in (3.1)–(3.2) is similar to that considered in the tokamak simulations of Giacomin & Ricci (Reference Giacomin and Ricci2022). Following the initial transient phase, the simulations reach a quasi-steady state where sources and losses balance each other. The time averaging and the statistical analysis are carried out over the final $\Delta t = 10\ R_0/c_{s0}$ interval in all the simulations.

3.2. Simulation results

Representative two-dimensional snapshots of the density and electrostatic potential on the $(R,\varphi)$ -midplane, obtained from four different GBS simulations characterised by $(\nu _0, \beta _{e0}) \in \{(0.05, 0.001), (5.0, 0.001), (50, 0.001), (0.05, 0.05)\}$ , are presented in figure 3. In the vicinity of the reversal surface, the $(R,\varphi)$ -midplane of RFP plasmas is almost perpendicular to the magnetic field and thus plays a role analogous to that of the poloidal plane in tokamak plasmas. Turbulence filaments develop in the outermost plasma region in proximity of the reversal surface. Turbulence structures observed on the high-field side ( $r\lt 0$ ) and on low-field side ( $r\gt 0$ ) are similar. This is expected since the projection of the curvature vector onto the minor radius direction in RFP plasmas is always positive, i.e. there is no good curvature region as in tokamaks. Although the magnetic flux surfaces are disrupted, the electrostatic potential and density still conform to the structure of the underlying helical magnetic field linked to the dominant $m=1,\, n=7$ resonant mode (see the reminiscent magnetic island around $r/a\simeq 0.4$ in figure 3). The GBS simulations indicate that the size and propagation of turbulent eddies around the $m=0, n=7$ magnetic island are essentially independent of the toroidal angle. In other words, they look the same whether they are located near the O-point or the X-point of the island. On the other hand, the electron temperature and density radial profiles are rather flat across the island. As a result, turbulence is only weakly driven inside the magnetic island, and the majority of turbulent structures are generated outside the island and merely pass through it. Nonetheless, this result is in qualitative agreement with experimental observations from the thermal helium beam diagnostic (Agostini et al. Reference Agostini, Scarin, Spizzo, Vianello and Carraro2014), where no clear difference is observed between filaments around the O-point and the X-point of the boundary magnetic island $(m, n)=(0, 7)$ .

At $\nu _0=0.05$ , $\beta _{e0}=0.001$ (reference case), turbulent eddies are relatively small and limited to the region $r/a\gtrsim 0.8$ (and the symmetric region $r/a\lesssim -0.8$ ). The magnitude of turbulent fluctuation levels in the plasma core is negligible. The electrostatic potential follows the ambipolarity response to stochastic transport, which generates an outward radial electric field proportional to the radial gradient of the electron temperature and density, as discussed in Giacomin et al. (Reference Giacomin, Momo, Predebon, Vianello and Zuin2026). The toroidal size of the turbulent eddies is larger at $\nu _0=5.0$ , and turbulence spreads over a slightly wider radial region. Despite the two orders of magnitude increase in resistivity, the two-dimensional density profile remains essentially unchanged relative to the reference case. The turbulent transport enhancement, limited to the region $r/a\gt 0.8$ , is insufficient to cause a significant change on density. On the other hand, the ambipolarity response of the radial electric field discussed in Giacomin et al. (Reference Giacomin, Momo, Predebon, Vianello and Zuin2026) is reduced by the increased resistivity through the term $j_\parallel /\sigma _\parallel$ in (2.3), leading to a weaker radial electric field. At $\nu _0=50.0$ , turbulent eddies extend well into the core region, yielding a large increase of turbulent transport, which degrades the overall plasma confinement. This state of very large turbulent transport that appears at large collisionality shares similarities with the degraded confinement regime found in tokamak simulations (Giacomin & Ricci Reference Giacomin and Ricci2020). However, this regime is attained here only at substantially larger values of $\nu _0$ than those typically required in tokamak simulations (see § 4 for a detailed discussion). Large-scale turbulent fluctuations appear also at large values of $\beta _{e0}$ . In this case, turbulence affects the entire plasma region. This regime is similar to the high- $\beta$ tokamak boundary turbulence regime identified in Giacomin & Ricci (Reference Giacomin and Ricci2022), which corresponds to the crossing of the ideal ballooning instability boundary, but substantially larger values of $\beta _{e0}$ are required in RFP boundary turbulence simulations to attain this turbulence regime. We note that the electrostatic potential is significantly modified at large $\nu _0$ and large $\beta _{e0}$ values. This behaviour arises because turbulent transport exceeds stochastic transport in the confined plasma region at sufficiently large values of $\nu _0$ and $\beta _{e0}$ , thereby reducing the ambipolarity response of the electric field.

The dependence of the particle and energy confinement times on $\nu _0$ and $\beta _{e0}$ , shown in figure 4, supports the previous observations. The particle and energy confinement times remain nearly constant as $\nu _0$ increases from 0.05 to 5, while they drop by almost a factor of two at $\nu _0=50$ . No dependence on $\beta _{e0}$ is observed until reaching $\beta _{e0} = 0.05$ , where the confinement times are reduced by approximately 30 %. The trend of the particle and energy confinement times is similar. Notably, boundary turbulence does not affect the overall plasma confinement at intermediate and moderate values of $\nu _0$ , marking a clear distinction with respect to tokamak boundary turbulence simulations (Giacomin & Ricci Reference Giacomin and Ricci2020, Reference Giacomin and Ricci2022). This simulation outcome is further discussed in § 4.

Figure 4.

Particle and energy confinement times computed from the $\nu _0$ (a) and $\beta _{e0}$ (b) simulation scans. The values are normalised to the ones computed from the reference simulation ( $\nu _0 = 0.05$ and $\beta _{e0} = 0.001$ ).

Figure 5.

Power spectrum of the density fluctuations as a function of $k_\varphi \rho _{s0}$ at various values of $\nu _0$ (a) and $\beta _{e0}$ (b), averaged over radius in the region $0.7\lt r/a\lt 1$ and over time in the quasi-steady-state phase. The power spectra are normalised to their maximum value. The zonal component ( $k_\varphi =0$ ) is excluded. The experimental power spectrum, taken from Rea et al. (Reference Rea2015), is represented by green star-shaped markers.

Figure 6.

Probability distribution function of the phase shift $\alpha _{\tilde {n}, \tilde {\phi }}$ between density and electrostatic potential fluctuations in the region $0.7\lt r/a\lt 1$ , weighted by the density fluctuation power spectrum and averaged over time in the quasi-steady-state phase, at various values of $\nu _0$ (a) and $\beta _{e0}$ (b). The experimental value at $k_\varphi \rho _s\simeq 0.1$ from Rea et al. (Reference Rea2015) is depicted as vertical green dashed line.

Figure 5 shows the $k_\varphi$ power spectrum of the density fluctuations, averaged over the radial interval $0.7 \lt r/a \lt 1$ and over time in the quasi-steady-state phase, from the $\nu _0$ and $\beta _{e0}$ simulation scans. The maximum of the power spectrum shifts at smaller $k_\varphi$ values as $\nu _0$ increases, although the power spectrum maximum is rather broad at $\nu _0=0.05$ and $\nu _0=0.5$ . The dependence of the dominant binormal wave vector on $\nu _0$ is consistent with resistive ballooning modes (RBMs), in agreement with Giacomin et al. (Reference Giacomin, Momo, Predebon, Vianello and Zuin2026), which identifies the RBM as the dominant instability in the reference case. The $\beta _{e0}$ scan reveals a weak dependence of the density power spectrum on $\beta _{e0}$ , except for the case at $\beta _{e0}=0.05$ , where the density power spectrum maximum shifts at the smallest finite $k_\varphi \rho _{s0}$ . Turbulence in the $\beta _{e0}=0.05$ simulation is most likely affected by the reduced toroidal domain size ( $L_\varphi = 2\pi /7$ ), which limits the minimum finite $k_\varphi \rho _{s0}$ value. Large-amplitude global modes develop across the full domain at large $\beta _{e0}$ , with toroidal mode numbers that are comparable to the dominant low- $n$ MHD tearing modes responsible for the helical magnetic field. An accurate description of turbulence at $\beta _{e0}\gtrsim 0.05$ would likely require (i) evolving the equilibrium magnetic field, while only the perpendicular magnetic field fluctuations are typically evolved in two-fluid turbulence codes, and (ii) extending the simulation domain to the full torus. Nonetheless, the high- $\beta _{e0}$ simulation is consistent with the analysis at high- $\beta _{e0}$ reported in Giacomin & Ricci (Reference Giacomin and Ricci2022) for tokamak boundary turbulence simulations.

3.3. Experimental comparison

The experimental power spectrum, derived from fluctuation measurements acquired by a diagnostic probe inserted into the plasma boundary region in similar conditions (Rea et al. Reference Rea2015; Spolaore et al. Reference Spolaore2017), is superimposed in figure 5. The experimental power spectrum peaks around $k_\varphi \rho _{s0}\simeq 0.05$ and agrees fairly well with the maximum of the numerical spectrum at $\nu _0=0.5$ and $\beta _{e0}= 0.001$ . The widths of the density fluctuation power spectra at $\nu _0 \simeq 0.05$ and $\nu _0 = 0.5$ are comparable to those obtained experimentally. However, the decay of the experimental power spectrum for $k_\varphi \rho _{s0} \gt 0.2$ is weaker than the corresponding decay observed in the numerical power spectra. This discrepancy may indicate the existence of an intermediate- $k_y$ microinstability that is not captured by the GBS simulations.

The cross-field $\boldsymbol{E}\times \boldsymbol{B}$ particle turbulent fluxes (and analogously for the turbulent heat flux) depends on the amplitude of the density and electrostatic potential fluctuations and on their phase shift, $\varGamma _{E,r} \propto k_\varphi |\tilde {n}| |\tilde {\phi }| \sin (\alpha _{\tilde {n},\tilde {\phi }})$ . At given fluctuation amplitudes, turbulent transport is maximum when $\alpha _{\tilde {n},\tilde {\phi }} = \pi /2$ , which is the case for purely ballooning mode-driven turbulence (Mosetto et al. Reference Mosetto, Halpern, Jolliet, Loizu and Ricci2013). Since the amplitudes of $\tilde {n}$ and $\tilde {\phi }$ are similar across the $\nu _0$ and $\beta _{e0}$ simulation scans, the phase shift $\alpha _{\tilde {n}, \tilde {\phi }}$ provides a direct measure of the impact of turbulent transport. The distribution probability of $\alpha _{\tilde {n}, \tilde {\phi }}$ in the region $0.7 \lt r/a \lt 1$ , averaged over time, is shown in figure 6 at various values of $\nu _0$ and $\beta _{e0}$ . Starting with the $\nu _0$ scan, the most probable value of $\alpha _{\tilde {n},\tilde {\phi }}$ in the simulations with $\nu _0=0.05$ and $\nu _0=0.5$ is approximately $\pi /4$ , which is in good agreement with the most probable value of $\alpha _{\tilde {n},\tilde {\phi }}$ measured experimentally (Rea et al. Reference Rea2015). The $\alpha _{\tilde {n},\tilde {\phi }}$ probability distribution function broadens as $\nu _0$ increases at $\nu _0\gt 0.5$ . At $\nu _0=50$ , the $\alpha _{\tilde {n},\tilde {\phi }}$ probability distribution peaks around $\alpha _{\tilde {n},\tilde {\phi }}\simeq \pi /2$ . In this case, turbulence is purely driven by RBMs. The parameter $\beta _{e0}$ weakly affects the probability distribution function of $\alpha _{\tilde {n}, \tilde {\phi }}$ , except at $\beta _{e0}=0.05$ , where it becomes broader and peaks at larger $\alpha _{\tilde {n}, \tilde {\phi }}$ values.

Boundary turbulence emerging from this set of RFP turbulence simulations exhibits properties analogous to those typically observed in tokamak boundary turbulence simulations. This is consistent with previous experimental studies that reported similar statistical characteristics of turbulence in the boundary region of RFP and tokamak plasmas (Vianello et al. Reference Vianello, Spolaore, Agostini, Cavazzana, De Masi, Martines, Momo, Scarin, Spagnolo and Zuin2016). On the other hand, turbulent transport in RFP simulations is localised in the proximity of the reversal surface and it impacts very weakly the overall plasma confinement. Figure 7 displays the radial electron pressure profile at the outboard midplane obtained at $\varphi =0$ from the $\nu _0$ and $\beta _{e0}$ simulation scans, averaged over time in the quasi-steady-state phase. The pressure profile is rather flat at $r/a\gtrsim 0.9$ at all values of $\nu _0$ and $\beta _{e0}$ , whereas the overall confinement is mostly set by the pressure gradient at $r/a \lt 0.8$ . Since turbulence develops only over a narrow plasma region $r/a\gtrsim 0.8$ , cross-field transport in the outer plasma core is mainly driven by other mechanisms, including the stochastic transport from magnetic chaos, which is independent of $\nu _0$ , thus explaining the weak dependence of the radial electron pressure profile on $\nu _0$ at small and intermediate $\nu _0$ values. At $\nu _0=50$ , turbulence strongly develops over a wide radial region. The enhanced turbulent transport causes a collapse of the pressure profile and the energy and particle confinement times drop substantially (see figure 4). This phenomenology closely resembles the behaviour observed in the vicinity of the tokamak density limit, as characterised in Giacomin et al. (Reference Giacomin, Pau, Ricci, Sauter and Eich2022a ).

Figure 7.

Electron pressure radial profile in the outer midplane at $\varphi =0$ from the $\nu _0$ (a) and $\beta _{e0}$ (b) simulation scans, averaged over time in the quasi-steady-state phase. The dashed lines represent an exponential fit of the electron pressure profile in the region $0.7\lt r/a\lt 1$ .

Figure 7. Long description

Panel A: A line graph shows the electron pressure radial profile in the outer midplane for different values of nu_0. The x-axis represents the normalized radius (r/a) ranging from 0.5 to 1.0, and the y-axis represents the normalized electron pressure (pe/p0) ranging from 0 to 10. The graph includes data points and dashed lines representing exponential fits for nu_0 values of 0.05, 0.5, 5.0, and 50.0. Panel B: Another line graph shows the electron pressure radial profile for different values of beta_e0. The x-axis represents the normalized radius (r/a) ranging from 0.5 to 1.0, and the y-axis represents the normalized electron pressure (pe/p0) ranging from 0 to 10. The graph includes data points and dashed lines representing exponential fits for beta_e0 values of 0.001, 0.005, 0.01, and 0.05.

The $\beta _{e0}$ simulation scan reveals a weak dependence of the edge pressure profile on $\beta _{e0}$ at $r/a\gt 0.7$ , even at large values of $\beta _{e0}$ . In contrast, a significant flattening of the pressure profile occurs at $r/a \lt 0.7$ in the simulation with $\beta _{e0}=0.05$ . This flattening is associated with the onset of ideal ballooning modes in the plasma core, where $\beta _e\propto p_e$ reaches values of the order of 20 %. At these large values of $\beta _{e0}$ , the magnetic flutter becomes important, enhancing the stochastic transport caused by the underlying magnetic chaos. A test simulation at even higher $\beta _{e0}$ shows destabilisation of ideal ballooning modes over the whole plasma volume, including the boundary region. Subsequently, the entire pressure profile is flattened. However, describing the plasma dynamics at these large $\beta _{e0}$ values would required evolving the full equilibrium magnetic field, as done in MHD simulations, and not just the magnetic field perturbation, as typically done in two-fluid turbulence codes.

In summary, while the RFP turbulence simulations presented here show properties similar to those typically observed in tokamak boundary turbulence simulations, they clearly point out a secondary role played by boundary turbulence in plasma confinement, marking a clear distinction with respect to tokamak boundary turbulence simulations carried out in the same geometry (Giacomin et al. Reference Giacomin, Vianello, Cavazzana, Molisani, Spolaore and Zuin2025). This aspect is further investigated in the following section.

4.1. Turbulence phase space in the RFP plasma boundary

The key quantity at the base of the theoretical framework developed in Giacomin & Ricci (Reference Giacomin and Ricci2022) is the equilibrium pressure gradient length $L_p$ across the separatrix. An analytical scaling of $L_p$ has been derived in Giacomin et al. (Reference Giacomin, Stagni, Ricci, Boedo, Horacek, Reimerdes and Tsui2021) by balancing the cross-field turbulent heat flux, the heat source and the parallel losses to the target plates, assuming (i) turbulence to be driven by RBMs and (ii) parallel transport in the scrape-off layer (SOL) to be mostly convective, leading toFootnote ¹

(4.1) \begin{align} L_p \simeq \frac {5^{8/17}\pi ^{12/17}}{2^{13/17}}\big(\rho _*^3\nu _0^6 q^{12}a^{12}(1+\kappa ^2)^6n^{10}S_p^{-4}\big)^{1/17}, \end{align}

where $L_p$ is normalised to $\rho _{s0}$ , $q$ is the edge safety factor, $a$ is the plasma minor radius normalised to $\rho _{s0}$ , $\kappa$ is the edge elongation, $n$ is the normalised separatrix density and $S_p= \rho _* \int s_p \mathrm{d}R\mathrm{d}Z$ is the normalised power source. Unless specified otherwise, dimensionless quantities are used in the following.

Equation (4.1) assumes that the binormal size of RBMs is limited by the parallel current, i.e. the curvature and the parallel gradient terms in (2.2) balance each other (Giacomin & Ricci Reference Giacomin and Ricci2020), $2\mathcal{C}(\tilde {p}_e) \sim \boldsymbol{\nabla} _\parallel \tilde {j}_\parallel$ . The pressure fluctuation amplitude, $\tilde {p}_e$ , is obtained from the linearised electron pressure equation, $\tilde {p}_e\simeq i k_y \tilde {\phi }p_e/(\rho _*L_p)$ , while the electrostatic potential fluctuation amplitude, $\tilde {\phi }$ , is estimated from the linearised Ohm’s law, $i k_\parallel \tilde {\phi }\sim \nu \tilde {j}_\parallel$ , where $k_y$ is the binormal wave vector (normalised to $1/\rho _{s0}$ ), $k_\parallel$ is the parallel wave vector (normalised to $1/R_0$ ) and $\gamma =\sqrt {2T_e/(\rho _*L_p)}$ is the maximum growth rate of the interchange instability. Therefore, $k_y$ is simply approximated by $k_y \sim k_\parallel /\sqrt {n\nu \gamma }$ . The parallel connection length of an axisymmetric toroidal plasma is defined as $L_\parallel = 2\pi R\sqrt {q^2+\epsilon ^2}$ , with $\epsilon =r/R$ . Within the boundary region of a tokamak, $L_\parallel$ reduces to $L_\parallel ^{\mathrm{tok}} \simeq 2\pi R_0 q$ , which yields $k_\parallel \simeq 1/q$ and, consequently, leads to (4.1). In contrast, in the boundary region of an axisymmetric RFP plasma, where $q \simeq 0$ , $L_\parallel$ can be approximated as $L_\parallel ^{\mathrm{RFP}} \simeq 2\pi a$ .

The GBS simulations presented in § 3 show that turbulence is mainly driven by RBMs. Therefore, we might expect the quasi-linear theory derived in Giacomin & Ricci (Reference Giacomin and Ricci2020) to hold in the boundary of RFP plasmas as long as (i) $k_\parallel$ is correctly evaluated and (ii) turbulence provides the dominant cross-field transport channel. By making the $k_\parallel$ dependence explicit in (4.1), we obtain

(4.2) \begin{align} L_p \simeq 2^{5/17} 5^{8/17}\pi ^{12/17}\big(\rho _*^3\nu _0^6 k_\parallel ^{-12}a^{12} n^{10}S_p^{-4}\big)^{1/17}, \end{align}

which is identical to (4.1) except for the dependence on the safety factor and on the elongation, which is taken as $\kappa =1$ (RFP plasmas are typically circular). The difference on the upfront numerical factor is partially due to $\kappa =1$ in (4.2) and partially comes from considering turbulent transport uniform on the poloidal angle, while, in Giacomin & Ricci (Reference Giacomin and Ricci2020), the turbulent flux is assumed to vanish on the high-field side, i.e. in the good curvature region. By neglecting the scattering of magnetic field lines due to magnetic chaos, $k_\parallel$ in the proximity of the reversal surface can be approximated by $k_\parallel \simeq 1/(\rho _*a)$ . This can be considered as an upper limit for $k_\parallel$ , since magnetic chaos and the helical structures breaking axisymmetry are likely to increase the parallel connection length (see § 4.2 for a discussion on the parallel connection length). We note that the primary difference between (4.2) and the $L_p$ scaling derived in Giacomin & Ricci (Reference Giacomin and Ricci2020) concerns $k_\parallel$ , whereas the effects of toroidal shaping associated with $n=7$ , $m \in \mathbb{N}$ tearing modes, as well as the resulting magnetic chaos, are neglected.

Three turbulence regimes have been identified in Giacomin & Ricci (Reference Giacomin and Ricci2022) depending on the value of plasma collisionality and $\beta _{e0}$ : a regime at low collisionality, where turbulence is driven by drift waves, a regime at intermediate and large collisionality, where turbulence is driven by RBMs, and a regime at high $\beta _{e0}$ , where turbulence is driven by ideal ballooning modes. In the RBM regime, turbulent transport increases with the edge collisionality, leading to a collapse of the pressure profile at large density, characterised by $L_p\sim a$ . This condition has been identified as the crossing of a density limit (Giacomin et al. Reference Giacomin and Ricci2022a ). Strongly enhanced turbulent transport and the subsequent flattening of the edge pressure profile are observed in the large $\nu _0$ simulation in § 3, thus suggesting a possible similar phenomenology for RFP plasmas. Hence, this limit is revisited here using the equilibrium pressure gradient length scaling in (4.2). The condition $L_p\sim a$ leads to

(4.3) \begin{align} \frac {\nu _0^{3/2}}{S_p} \sim \frac {1}{2^{5/4}25\pi ^3}\frac {a^{5/4}k_\parallel ^3}{\rho _*^{3/4}n^{5/2}}. \end{align}

Except for the upfront numerical factor, (4.3) is equivalent to (45) of Giacomin & Ricci (Reference Giacomin and Ricci2022) if $k_\parallel =1/q$ and $\kappa =1$ . We emphasise that $\nu _0^{3/2}/S_p$ at the transition depends strongly on $k_\parallel$ . The corresponding density limit condition for RFP plasmas is obtained by substituting $k_\parallel = 1/(\rho _*a)$ in (4.3), leading to

(4.4) \begin{align} \frac {\nu _0^{3/2}}{S_p} \sim \frac {1}{2^{5/4}25\pi ^3}\frac {1}{\rho _*^{15/4}n^{5/2}a^{7/4}}. \end{align}

Figure 8(a) illustrates the dependence of the ratio between the numerically evaluated pressure gradient scale length, obtained from GBS simulations by fitting the radial profiles of the electron pressure in the interval $0.7 \lt r/a \lt 1$ (see figure 7), and the analytical prediction given by (4.2), as a function of the parameter $\nu _0^{3/2}/S_p$ . At low $\nu _0^{3/2}/S_p$ , the numerical $L_p$ exceeds significantly the value predicted by (4.2), while the ratio $L_{p,\mathrm{GBS}}/L_p$ approaches unity at large $\nu _0^{3/2}/S_p$ . This is consistent with RBM-driven turbulent transport becoming increasingly important in the region $0.7 \lt r/a \lt 1$ as $\nu _0$ increases. In other words, the $L_p$ scaling in (4.2) strictly applies only at $\nu _0^{3/2}/S_p\gt 0.1$ , when RBM turbulence is the dominant cross-field transport channel. It is worth noting that, because of the large $k_\parallel$ value, (4.2) predicts a very steep pressure gradient at low $\nu _0$ , i.e. the pressure profile would be much steeper if exclusively limited by RBM-driven turbulent transport. Analogously, figure 8(b) shows the energy confinement time reduction associated with the increase of the RBM-driven turbulent transport. The energy confinement time drops by a factor of two before crossing the limit predicted by (4.4), which is derived from (4.3) by assuming $k_\parallel = 1/(\rho _*a)$ (corresponding to $k_\parallel R_0 = R_0/a=4$ for a RFX-mod plasma). This is expected as (4.4) provides only an upper limit. The value of $k_\parallel =2\pi /L_\parallel$ might be smaller than that estimated from magnetic field geometry considerations (see the discussion in § 4.2). For instance, the limit occurs at much smaller $\nu _0^{3/2}/S_p$ values with $k_\parallel =1$ , and the simulation with $\nu _0=50$ would be above this limit.

The phase space of boundary turbulence in figure 8 retrieves the connection between boundary turbulence and energy confinement previously found in tokamak turbulence simulations (Giacomin & Ricci Reference Giacomin and Ricci2022). On the other hand, the density limit in figure 8 is crossed only at unrealistically large values of $\nu _0$ , even when $k_\parallel = 1$ . In fact, the density limit previously derived for tokamak plasmas is approximately two orders of magnitude smaller than that predicted by (4.4)

(4.5) \begin{align} \frac {\big(\nu _0^{3/2}/S_p\big)_{\mathrm{RFP}}}{\big(\nu _0^{3/2}/S_p\big)_{\mathrm{tokamak}}} \sim \left(\frac {q}{\rho _*a}\right) \gtrsim 100. \end{align}

This result clearly questions the applicability of this theoretical framework at the edge region of RFP plasmas, where transport mechanisms other than RBM-driven turbulence appear to control the edge pressure gradient at realistic values of $\nu _0$ .

Figure 8.

Phase space of RFP boundary turbulence. Each dot in the figures corresponds to a GBS simulation, while the $\beta _{e0}$ and $\nu _0$ scans are identified with a label. The colour scale represents the ratio of the numerical to analytical edge pressure gradient length (a) and the energy confinement time (b). The vertical and horizontal dashed lines indicate the limits predicted by (4.3) and (4.7), respectively, for two different values of $k_\parallel$ .

The phase space of tokamak edge turbulence of Giacomin & Ricci (Reference Giacomin and Ricci2022) provides also a limit for the maximum achievable value of $\beta _{e0}$ before a substantial reduction of the energy confinement time is observed in tokamak turbulence simulations. This limit was derived from the parameter $\alpha _{\mathrm{MHD}} = -R_0 q^2 \mathrm{d}\beta /\mathrm{d}r$ , which is ill defined in the boundary of RFP plasmas where $q$ vanishes. Alternatively, this limit can be derived directly from the condition corresponding to the onset of ideal ballooning modes, which is obtained from the generalised Ohm’s law (see (2.3)), by balancing the electromagnetic induction and the resistive current terms

(4.6) \begin{align} \gamma \psi \sim \nu \tilde {j}_\parallel , \end{align}

where $\gamma$ is the maximum growth rate of the ballooning instability, $\gamma = \sqrt {2T_e/(\rho _* L_p)}$ . By using (2.8), the magnetic fluctuation $\psi$ can be approximated by $\psi \sim \beta _{e0} j_\parallel /(2 k_\perp ^2)$ , where $k_\perp$ is the perpendicular wave vector. In the boundary region of RFP plasmas, $k_\perp$ can be replaced by $k_\varphi$ . Furthermore, before the onset of ideal ballooning modes, turbulence is driven by RBMs and $k_\varphi ^2$ can be approximated by $k_\varphi ^2 \simeq k_\parallel ^2/(n\nu \gamma )$ . Within these assumptions, (4.6) leads to $\beta _{e0}\sim k_\parallel ^2\rho _* L_p/p_e$ . The $\beta$ limit is finally obtained by replacing $L_p$ with (4.2)

(4.7) \begin{align} \frac {\beta _{e0}S_p^{18/17}}{\nu _0^{10/17}} \sim 2^{3/17}5^{2/17}\pi ^{55/51}\rho _*^{22/17}n^{11/17}k_\parallel ^{14/17}a^{20/17}, \end{align}

where the key parameters are isolated in the left-hand side of (4.7). Except for the numerical pre-factor, (4.7) is identical to (50) of Giacomin & Ricci (Reference Giacomin and Ricci2022) if $k_\parallel = 1/q$ and $\kappa = 1$ .

As discussed in § 3, turbulent transport in the region $0.7\lt r/a\lt 1$ is weakly affected by the parameter $\beta _{e0}$ , and the GBS simulations belonging to the $\beta _{e0}$ scan return similar values of $L_{p, \mathrm{GBS}}$ (see figure 7). Consequently, figure 8(a) shows that $L_{p, \mathrm{GBS}}/L_p$ is nearly independent from $\beta _{e0}$ and it largely exceeds unity. Moreover, figure 8(b) shows that the energy confinement time drops as the $\beta$ limit, given by (4.7), is approached, providing support to the applicability of (4.7). We note that, the $\beta$ limit predicted by (4.7) is less sensitive on the $k_\parallel$ assumption than the density limit predicted by (4.3).

4.2. Limits of the two-fluid model approach in RFP

The statistical properties of turbulence fluctuations in the plasma boundary of the GBS simulations presented in § 3 compare reasonably well with experiments. Although a more comprehensive validation would require full-scale GBS simulations (at least in the poloidal plane), this comparison provides supporting evidence for the reliability of the turbulence properties predicted by GBS simulations. Such RBM-driven turbulence is also observed in RFX-mod GBS tokamak simulations (Giacomin et al. Reference Giacomin, Vianello, Cavazzana, Molisani, Spolaore and Zuin2025), thus proving a strong connection between RFP and tokamak boundary turbulence, which supports the experimental finding of Vianello et al. (Reference Vianello, Spolaore, Agostini, Cavazzana, De Masi, Martines, Momo, Scarin, Spagnolo and Zuin2016). Furthermore, the turbulence phase space developed in § 4.1 for RFP GBS turbulence simulations shares the same theoretical framework as the one developed for tokamak edge turbulence (Giacomin & Ricci Reference Giacomin and Ricci2022). These results underscore the significant insights that two-fluid, global, flux-driven turbulence simulations can offer into disentangling the complex turbulent dynamics operating at the edge region of RFP plasmas. However, there are a few criticalities. First, the equilibrium profiles and the energy confinement depend weakly on $\nu _0$ and $\beta _{e0}$ at realistic values of $\nu _0$ and $\beta _{e0}$ , while these are key parameters controlling tokamak boundary turbulence. Second, turbulence develops only in a very narrow region at the plasma boundary unless unrealistically large values of $\nu _0$ and $\beta _{e0}$ are considered. Finally, the crossing of the density limit predicted by (4.4) occurs at unrealistically large values of density.

The fluid closure used in Braginskii models requires the electron mean free path, $\lambda _e = v_{th,e}\tau _e$ ( $v_{th,e} =\sqrt {T_e/m_e}$ is the thermal electron velocity and $\tau _e = 0.51 m_e \sigma _\parallel /(ne^2)\propto \nu ^{-1}$ is the electron collision time), to be smaller than the typical length scale in the parallel direction, known as the parallel connection length, $L_\parallel$ . This condition is typically satisfied in present-day tokamak devices, although kinetic effects, which are not captured by a collisional closure, may still play an important role (Zholobenko et al. Reference Zholobenko, Body, Manz, Stegmeir, Zhu, Griener, Conway, Coster and Jenko2021; Oliveira et al. Reference Oliveira2022; Pitzal et al. Reference Pitzal, Stegmeir, Zholobenko, Zhang and Jenko2023; Frei et al. Reference Frei, Ulbl, Pitzal, Zholobenko and Jenko2025).

Figure 9.

(a) Two-dimensional snapshots of density fluctuations on the $(\theta , \varphi )$ surface at two different radial locations from the reference simulation. The black curves represent the magnetic field lines at the surface. (b) Comparison between the ensemble-averaged parallel correlation length and the electron mean free path obtained from simulations with $\nu _0=0.05$ , $\beta _{e0} = 10^{-3}$ and with $\nu _0=0.5$ , $\beta _{e0} = 10^{-3}$ . The shaded region denotes the standard deviation of $L_\parallel$ , while the black dotted line marks $2\pi r$ .

The definition of a parallel connection length in RFP plasmas is non-trivial. In the vicinity of the reversal surface, simple magnetic field geometry considerations would indicate an $L_\parallel$ value of $2\pi a$ . A turbulent eddy localised on the reversal surface forms a closed structure after completing a single poloidal turn, therefore constraining its parallel extension to $2\pi a$ . However, when moving radially inwards from the reversal surface, additional factors must be taken into account, including the helical deformation of the magnetic field and the presence of magnetic field chaos. Another key factor that enhances ballooning modes in RFP plasmas is the lack of regions with favourable magnetic curvature. Consequently, the commonly invoked argument that ballooning modes extend along magnetic field lines until they encounter a region of good curvature, where they are stabilised, does not apply to RFP configurations. In this geometry, the parallel correlation length of RBMs can exceed one full poloidal transit without the modes being quenched on the high-field side. This feature effectively increases the role of ballooning instabilities in the plasma, thus providing some support to the GBS results. Nonetheless, it is reasonable to assume that $L_\parallel$ at the edge of RFP plasmas remains shorter than $L_\parallel$ in tokamak plasmas. Hence, at similar values of electron temperature and density (similar $\lambda _e$ ), the collisionality in RFP plasmas is expected to be smaller than the one in tokamaks.

The parallel correlation length of turbulent fluctuations can be directly obtained from the simulation results. Figure 9(a) presents two snapshots of $\tilde {n}$ on the $(\theta , \varphi )$ surface, taken at two distinct radial positions of the reference simulation. On both flux surfaces, the density perturbations display a characteristic parallel scale that extends beyond one poloidal turn. The parallel correlation length, which can be effectively identified with the characteristic parallel scale $L_\parallel$ , is here defined as the parallel spatial separation between two consecutive zeros of the density fluctuation field. The ensemble average of $L_\parallel$ is displayed in figure 9(b) at different radial positions. The geometrical estimate $L_\parallel \simeq 2\pi r$ is also displayed in figure 9(b) for reference. Although $L_\parallel$ exhibits substantial scatter, the majority of the fluctuation structures possess a parallel correlation length that exceeds the value $2\pi r$ by a factor of approximately 2–3. This yields a value of $k_\parallel$ such that $1 \lt k_\parallel R_0 \leqslant 4$ , thereby justifying the bounds adopted in § 4.1.

The surface-averaged value of $\lambda _e$ obtained from the reference simulation, together with that computed from the simulation with $\nu _0 = 0.5$ and $\beta _{e0} = 10^{-3}$ , is also presented in figure 9(b). The electron mean free path in the reference simulation exceeds the parallel connection length across the entire radial region of figure 9(b), although the significant scattering affecting $L_\parallel$ yields comparable values in the region $0.8\lesssim r/a \lesssim 0.85$ . Conversely, the simulation with higher $\nu _0$ , i.e. higher $n_0$ , yields $\lambda _e\lt L_\parallel$ . In the latter case, $\lambda _e$ is smaller than the geometrical $L_\parallel$ . This is expected as $\lambda _e$ is inversely proportional to $\nu _0$ . Hence, the condition $L_\parallel /\lambda _e\gt 1$ is fulfilled for values of $\nu _0$ larger than the reference one. Moreover, at $\nu _0\gtrsim 5$ , $\lambda _e$ becomes significantly smaller than $L_\parallel$ over a wide radial region, thus allowing RBM turbulence to easily develop in the core plasma region and therefore controlling the overall plasma confinement. This explains why the energy confinement time is found to depend on boundary turbulence only at (unrealistically) large values of $\nu _0$ . It is worth remarking that the $k_\parallel$ stabilisation mechanism remains strong even at large $\nu _0$ , in agreement with (4.2). Notably, the stabilisation mechanism associated with finite $k_\parallel$ is sufficiently strong that the pressure gradient is primarily limited by transport channels other than RBM-induced turbulent transport, except in the regime where $\nu _0$ is increased to unrealistically large values. This occurs despite the fact that RBM turbulence does develop and produces power-density fluctuation spectra that exhibit reasonably good quantitative agreement with experimental measurements (see figure 5).

In conclusion, the short parallel connection length in RFP plasmas clearly restricts the validity of two-fluid turbulence simulations. This result calls for future investigations with more advanced modes, such as the one recently implemented in GRILLIX (Pitzal et al. Reference Pitzal, Stegmeir, Zholobenko, Zhang and Jenko2023).

Figure 10.

Growth rate (a) and real frequency (b) as functions of $k_y\rho _s$ from GENE local linear simulations at three different values of the electron collision frequency. The reference value of the electron collision frequency, which is directly evaluated from the electron temperature and density profiles of the GBS reference simulation, is $\nu _{e, \mathrm{ref}} \simeq 0.6\ c_s/a$ . Only unstable modes are shown.

Figure 10. Long description

The image contains two line graphs. Panel A shows the growth rate as a function of k_y rho_s. The vertical axis is labeled (a/c_s) gamma and the horizontal axis is labeled k_y rho_s. Four different data series are plotted, each representing a different electron collision frequency. The blue line represents nu_e = (1/2) nu_e, ref, the orange line represents nu_e = nu_e, ref, the green line represents nu_e = (3/2) nu_e, ref, and the red line represents nu_e = 2 nu_e, ref. The growth rate peaks and then decreases for each series. Panel B shows the real frequency as a function of k_y rho_s. The vertical axis is labeled (a/c_s) omega and the horizontal axis is labeled k_y rho_s. The same four data series are plotted as in Panel A. The real frequency increases for each series as k_y rho_s increases.