2 The parallel boundary condition

The parallel boundary condition generally used by nonlinear flux tube simulations is called the ‘twist-and-shift’ condition (Beer et al. Reference Beer, Cowley and Hammett1995). As with the perpendicular boundaries, it requires the two ends of the flux tube to have statistically identical turbulence. Finding statistically identical turbulence is non-trivial because toroidicity and plasma shaping mean that the statistical properties of turbulence change with poloidal angle. Thus, to be statistically identical the two ends of the flux tube must be at the same poloidal angle $\unicode[STIX]{x1D712}$, although, due to axisymmetry, they are still free to be at different toroidal angles $\unicode[STIX]{x1D701}$. This constraint requires the length of the domain to be an integer number of poloidal circuits, which we represent by $N_{\text{pol}}$.

While $N_{\text{pol}}$ is free to be any positive integer, it has become overwhelmingly standard in the fusion community to use $N_{\text{pol}}=1$ (e.g. figure 1). While this is likely due to the desire to minimize the computational cost of simulations, there also have been concerns about the validity of flux tubes with $N_{\text{pol}}>1$. For example, such simulations are often not ‘globally consistent’ (Scott Reference Scott1998). This means that, when $N_{\text{pol}}>1$, the narrow flux tube can permit Fourier modes that do not exist on the full, doubly periodic flux surface. However, in the true $\unicode[STIX]{x1D70C}_{\ast }\rightarrow 0$ limit for which gyrokinetics is valid, the charged particles have no way of knowing that they are on a doubly periodic flux surface. Particles can never communicate information between the different poloidal turns by moving in the toroidal direction because the distance separating the different poloidal turns is proportional to the minor radius $a$, while the turbulent correlation length is proportional to the ion gyroradius $\unicode[STIX]{x1D70C}_{i}$ (see figure 2). Thus, the distance separating the poloidal turns is an asymptotically large number of turbulent correlation lengths, so the particles will perceive the perpendicular direction within the flux surface as infinite in the $\unicode[STIX]{x1D70C}_{\ast }\equiv \unicode[STIX]{x1D70C}_{i}/a\rightarrow 0$ limit. This means that it is not relevant which modes are allowed on a doubly periodic flux. While we should still be concerned about respecting any periodicity in the parallel direction (i.e. do not use a domain with $N_{\text{pol}}=2$ to model a $q=1$ surface), violating global perpendicular periodicity will only introduce an error that is small in $\unicode[STIX]{x1D70C}_{\ast }\ll 1$. This is intuitive – global consistency should be unimportant in local simulations. Another concern is that flux tubes with $N_{\text{pol}}>1$ include multiple locations at the same poloidal angle, which must all have the same value of the zonal flows (i.e. fluctuations with a binormal wavenumber $k_{y}=0$). Fortunately, it appears that the physics of the gyrokinetic model ensures that this is fairly well satisfied in the simulations we have run (e.g. figure 3). However, we should point out that we have no proof that this must always be the case, so it remains an open issue.

Figure 2.

An example flux surface showing two poloidal turns of a magnetic field line (red arrows) and a global perturbation (grey and white stripes) with toroidal and poloidal mode numbers of 9 and 21 respectively. Also shown is a $N_{\text{pol}}=2$ flux tube (thin blue and purple lines) with the parallel boundary condition applied between the two ends (thick black horizontal lines). As shown, $\unicode[STIX]{x1D712}$ is taken to be a straight-field line poloidal angle.

Figure 3.

The zonal amplitude of the turbulent electrostatic potential $\unicode[STIX]{x1D719}$ as a function of the time $t$, as calculated by a gyrokinetic simulation with the parameters of tables 1 and 2 using kinetic electrons, a binormal domain width $L_{y}=125\unicode[STIX]{x1D70C}_{i}$, and $N_{\text{pol}}=3$. The colour (i.e. black or grey) indicates two different values of the radial wavenumber (i.e. $k_{x}\unicode[STIX]{x1D70C}_{i}\in \{0.05,0.1\}$ respectively) and the symbol (i.e. line, circle or cross) indicates the three different toroidal locations at the outboard midplane (i.e. $\unicode[STIX]{x1D712}\in \{-2\unicode[STIX]{x03C0},0,2\unicode[STIX]{x03C0}\}$). Here, the reference velocity $v_{r}\equiv \sqrt{T_{r}/m_{r}}$, temperature $T_{r}$ and mass $m_{r}$ are taken to be the ion values and $e$ is the proton charge.

Table 1.

The nominal GENE (Jenko et al. Reference Jenko, Dorland, Kotschenreuther and Rogers2000) coordinate grids used for the scans in $L_{y}$ and $N_{\text{pol}}$, where $L_{y0}\equiv 125\unicode[STIX]{x1D70C}_{r}$ is a common value used by the community. Note that all grids are equally spaced, the reference length $a$ is the tokamak minor radius, the reference magnetic field $B_{r}$ is the toroidal field in the flux tube at $R_{0}$, $R_{0}$ is the average major radius of the flux tube, and the CFL condition is explained in Courant, Friedrichs & Lewy (Reference Courant, Friedrichs and Lewy1967). The other reference quantities indicated by the subscript $r$ are taken to be the ion quantities (except for the electron-scale simulation in figure 6) and $v_{s}=\sqrt{2T_{s}/m_{s}}$ is the thermal velocity of species $s$.

Table 2.

The nominal Cyclone Base Case (Dimits et al. Reference Dimits, Bateman, Beer, Cohen, Dorland, Hammett, Kim, Kinsey, Kotschenreuther and Kritz2000) parameters used for the scans in $L_{y}$ and $N_{\text{pol}}$. Simulations with kinetic electrons use $\unicode[STIX]{x1D716}_{v\Vert }=\unicode[STIX]{x1D716}_{c}=0.5$ (Pueschel, Dannert & Jenko Reference Pueschel, Dannert and Jenko2010) and a plasma beta of $\unicode[STIX]{x1D6FD}=0.001$ for numerical reasons.

Unlike with the perpendicular boundary conditions, substituting statistical periodicity with exact periodicity across the parallel boundary is made more complicated by the effect of magnetic shear. In local simulations, the safety factor profile is approximated to be linear across the flux tube, i.e.

(2.1)$$\begin{eqnarray}q(x)=q_{0}+\frac{\text{d}q}{\text{d}x}(x-x_{0})=q_{0}\left(1+{\hat{s}}\frac{x-x_{0}}{x_{0}}\right),\end{eqnarray}$$

where $q$ is the safety factor, $x\equiv r$ is the minor radial coordinate within the domain, ${\hat{s}}\equiv (x_{0}/q_{0})\,\text{d}q/\text{d}x$ is the magnetic shear and the subscript 0 indicates the quantity is evaluated at the centre of the radial domain. This means that the two ends of the flux tube are deformed into parallelograms that have opposite tilts. Hence, the shapes do not overlap nicely, which must be resolved by the parallel boundary condition (see figure 4). In essence, the parallel boundary condition determines how the magnetic field lines are connected between these two parallelograms. The twist-and-shift condition typically implemented by gyrokinetic codes starts by using axisymmetry to justify toroidally shifting the two parallelograms until they are centred on one another (Beer et al. Reference Beer, Cowley and Hammett1995). Then, again using axisymmetry, we can make periodic copies of the flux tube, shifted in the toroidal direction, and completely cover the opposing parallelogram (see figure 4). Thus, any turbulent structure that extends past one end of the domain will be copied into the other end while maintaining the twisting caused by magnetic shearFootnote ³.

One consequence of this boundary condition is that it forces at least one flux surface in the domain to have field lines that close on themselves after going through the parallel boundary just once or, equivalently, after just $N_{\text{pol}}$ poloidal turns. Given that we have chosen to centre the parallelograms with respect to one another, the flux surface in the centre of the radial domain will always be such a surface (see figure 5). However, figures 4 and 5 show that there can be more than one – sometimes there are many. We will indicate the total number of surfaces with field lines that close on themselves after just one pass through the parallel boundary with $N_{\text{si}}$. The subscript $si$ stands for ‘self-interaction’ because these surfaces will exhibit the strongest parallel self-interaction. Now if $N_{\text{pol}}=1$, the parallel boundary condition will be applied after one poloidal turn, so these field lines will close after just one poloidal turn. Accordingly, we will call them ‘pseudo-integer’ surfaces – they close on themselves after one poloidal turn like a normal integer surface, but they are usually an artefact of the boundary condition and do not actually correspond to integer flux surfaces in the real device. Additionally, we see that there will be ‘pseudo-rational’ surfaces at other radial locations, where the field lines exactly close on themselves after two or more poloidal turns. For example, the ‘pseudo-half-integer’ surfaces close on themselves after two poloidal turns and will always occur midway between two neighbouring pseudo-integer surfaces. Alternatively, if $N_{\text{pol}}$ is increased from 1 to 2, all of the pseudo-integer surfaces will be replaced by pseudo-half-integer surfaces. Thus, the lowest-order rational surfaces in an $N_{\text{pol}}=2$ simulation will be pseudo-half-integer and there will be $N_{\text{si}}$ of them.

Figure 4.

The shape of the parallel boundary (thin solid black parallelogram) on the inboard midplane of the example flux tube shown in figure 1. The parallel boundary condition is supplied by shifting the other end of the flux tube (thick solid black parallelogram) toroidally until it is centred on the first and then making periodic copies offset in the toroidal direction (two thick dashed black parallelograms). The flux surfaces are also shown (thin solid grey lines).

Figure 5.

A Poincaré puncture plot at the outboard midplane for the flux tube from figures 1 and 4, which has $N_{\text{pol}}=1$. Following a magnetic field line that starts at different radial positions (smallest filled circles) will take us to different toroidal locations after one poloidal circuit of the device (smaller empty circles) and two poloidal circuits (large empty circles). Two radial locations have field lines that close on themselves after passing through the parallel boundary just once, so $N_{\text{si}}=2$.

These pseudo-rational surfaces are also involved in an important constraint that arises from the combination of the parallel and radial boundary conditions – the aspect ratio of the domain becomes discretized (i.e. the radial domain width $L_{x}$ divided by the binormal domain width $L_{y}$). This can be found from a derivation in real space (see appendix A), a derivation in Fourier space (Beer et al. Reference Beer, Cowley and Hammett1995) or by the following explanation. First, note that radial periodicity implies an equivalence between pairs of field lines, one on each of the two radial boundaries. In other words, any field line at $x=x_{0}-L_{x}/2$ has a matching field line at $x=x_{0}+L_{x}/2$ with which it is equivalent. These field lines are the same and must remain matched as we move in the parallel direction. Otherwise, a particle moving purely along a field line would find itself transported across field lines. However, from figures 4 and 5, we see that the twist-and-shift condition, which is applied at the parallel ends of the domain, connects different field lines at different toroidal locations. This can cause problems. If we are not careful, the twist-and-shift boundary condition will connect the two field lines from a matching pair to two other field lines that are not matching. This would cause spurious cross-field transport. We must ensure that, given any matching pair of field lines, twist-and-shift always connects them to two other field lines that also form a matching pair. From figures 4 and 5, we see that this will only be achieved if

(2.2)$$\begin{eqnarray}\unicode[STIX]{x0394}\unicode[STIX]{x1D701}=N_{\text{si}}L_{\unicode[STIX]{x1D701}},\end{eqnarray}$$

where $\unicode[STIX]{x0394}\unicode[STIX]{x1D701}$ is the total shift in toroidal angle caused by magnetic shear along the length of the flux tube, $N_{\text{si}}$ is some integer, and $L_{\unicode[STIX]{x1D701}}$ is the width of the flux tube in toroidal angle such that $L_{\unicode[STIX]{x1D701}}=2\unicode[STIX]{x03C0}$ corresponds to the full toroidal domain. For the example shown in figure 4, we can see that $N_{\text{si}}=2$ because the radial boundaries of the thick solid black parallelogram are offset in the toroidal direction by $1\times L_{y}$ and this represents half the total shift (i.e. the shift in the range $\unicode[STIX]{x1D712}\in [-\unicode[STIX]{x03C0}N_{\text{pol}},0]$, but not $\unicode[STIX]{x1D712}\in [0,\unicode[STIX]{x03C0}N_{\text{pol}}]$). Additionally, figure 5 shows that, for this example, the twist-and-shift condition connects every field line on the radial boundary to itself. This trivially ensures that matching pairs remain matching as particles move across the parallel boundary. If $N_{\text{si}}=1$, the field lines on the radial boundary would no longer close on themselves, but instead to a different matching pair offset in the $y$ direction by $L_{y}/2$.

To see how equation (2.2) discretizes the domain aspect ratio, we invert the definition of the binormal coordinate

(2.3)$$\begin{eqnarray}y(x,\unicode[STIX]{x1D701},\unicode[STIX]{x1D712})\equiv \frac{x_{0}}{q_{0}}(q(x)\unicode[STIX]{x1D712}-\unicode[STIX]{x1D701})\end{eqnarray}$$

(where $\unicode[STIX]{x1D712}$ is a straight-field line poloidal angle) and use (2.1) to calculate the toroidal shift across the sheared domain to be

(2.4)$$\begin{eqnarray}\displaystyle \unicode[STIX]{x0394}\unicode[STIX]{x1D701} & \equiv & \displaystyle \unicode[STIX]{x1D701}\left(x=x_{0}+\frac{L_{x}}{2},y=y_{0},\unicode[STIX]{x1D712}=2\unicode[STIX]{x03C0}N_{\text{pol}}\right)-\unicode[STIX]{x1D701}\left(x=x_{0}-\frac{L_{x}}{2},y=y_{0},\unicode[STIX]{x1D712}=2\unicode[STIX]{x03C0}N_{\text{pol}}\right)\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(2.5)$$\begin{eqnarray}\displaystyle & = & \displaystyle 2\unicode[STIX]{x03C0}N_{\text{pol}}\frac{\text{d}q}{\text{d}x}L_{x}=2\unicode[STIX]{x03C0}N_{\text{pol}}\frac{q_{0}}{x_{0}}{\hat{s}}L_{x},\end{eqnarray}$$

where $y_{0}$ is any arbitrary binormal location. Since

(2.6)$$\begin{eqnarray}L_{y}=\frac{x_{0}}{q_{0}}L_{\unicode[STIX]{x1D701}}\end{eqnarray}$$

follows from evaluating equation (2.3) at constant $x$ and $\unicode[STIX]{x1D712}$, we can use (2.2) and (2.5) to see that the possible values for the radial domain width are

(2.7)$$\begin{eqnarray}L_{x}=\frac{N_{\text{si}}}{2\unicode[STIX]{x03C0}N_{\text{pol}}|{\hat{s}}|}L_{y}.\end{eqnarray}$$

Thus, we find that $L_{x}$ must be a integer multiple of the distance between lowest-order pseudo-rational surfaces. Importantly, we see that if the magnetic shear is very small, the radial width will be forced to be large because we cannot decrease $N_{\text{si}}$ below 1. This causes simulations with very low, but finite magnetic shear to become expensive. Equation (2.7) also shows that $N_{\text{si}}$ can be used to specify the radial width of the domain. In fact, $N_{\text{si}}$ is a common input parameter to gyrokinetic codes and goes by various names such as nexc in GENE, jtwist in GS2, IKXSPACE in GKW and m_j in GKV.

Lastly, it is important to note that the pseudo-rational surfaces created by the parallel boundary condition can be physical (Waltz et al. Reference Waltz, Austin, Burrell and Candy2006; Weikl et al. Reference Weikl, Peeters, Rath, Seiferling, Buchholz, Grosshauser and Strintzi2018; C.J. et al. Reference C.J., Brunner, McMillan, Ball, Dominski and Merlo2020). For example, a $q=3$ surface really closes on itself after one poloidal turn and should be modelled using $N_{\text{pol}}=1$, while a $q=5/2$ surface should use $N_{\text{pol}}=2$. Moreover, the physics occurring at these surfaces may be crucial for modelling certain behaviour such as microtearing modes (Hazeltine, Dobrott & Wang Reference Hazeltine, Dobrott and Wang1975; Guttenfelder et al. Reference Guttenfelder, Candy, Kaye, Nevins, Wang, Bell, Hammett, LeBlanc, Mikkelsen and Yuh2011). However, to accurately model the overall impact of self-interaction from rational surfaces, the spacing between pseudo-rational surfaces must also be correct. For example, the domain shown in figure 5 has $N_{\text{si}}=2$ pseudo-integer surfaces, so, to be physical, the radial width of the simulation domain must correspond to twice the distance between the actual integer surfaces in the deviceFootnote ⁴. To put it another way, the simulation must have the same number of ion gyroradii separating the integer surfaces as the experiment does. Unfortunately, accomplishing this turns out to require a simulation domain that corresponds to the full flux surface (as is discussed further at the end of § 4). To see this, take the change in safety factor across the simulation domain $\unicode[STIX]{x0394}q=L_{x}\,\text{d}q/\text{d}x$ and set it equal to the number of integer surfaces in the simulation $\unicode[STIX]{x0394}q=N_{\text{si}}$. Then, substituting (2.6) and (2.7) with $N_{\text{pol}}=1$ shows that $L_{\unicode[STIX]{x1D701}}=2\unicode[STIX]{x03C0}$ (i.e. the toroidal domain size is equal to the entire toroidal domain). Thus, we see that accurately modelling self-interaction on rational surfaces requires full flux surface simulations with radial domain widths on the scale of the tokamak minor radius $a$. Flux tubes of this size offer little computational savings compared to global simulations, but if we try to get by with a smaller domain of $L_{\unicode[STIX]{x1D701}}=2\unicode[STIX]{x03C0}/N$, it will have a radial density of pseudo-integer surfaces that is $N$ times larger than reality.

3 Neglecting pseudo-rational surfaces

While the pseudo-rational surfaces created by the parallel boundary condition are unphysical unless the domain size is very large, this is not always a problem. Typically gyrokinetic simulations are used to model plasma far from low-order rational surfaces. To do this accurately we do not need to ensure that the simulation has the correct spacing of rational surfaces, we just need to ensure that their presence is negligible. Broadly speaking, this can be accomplished by verifying that the results of the simulation do not depend on the spacing of pseudo-rational surfaces. From (2.7), we see that the spacing is given by

(3.1)$$\begin{eqnarray}\frac{L_{x}}{N_{\text{si}}}=\frac{L_{y}}{2\unicode[STIX]{x03C0}N_{\text{pol}}|{\hat{s}}|},\end{eqnarray}$$

so we can vary either $L_{y}$ or $N_{\text{pol}}$ and test for convergence. Additionally, we can specifically check the low-order pseudo-rational surfaces to see if they are exhibiting any unusual behaviour. The low-order pseudo-rational surfaces are only distinct from the other surfaces in that they close on themselves after a small number of poloidal turns. If the parallel correlation length of the turbulence is sufficiently long, an eddy can interact with itself along the field line, thereby altering its statistical properties. Hence, to ensure that their presence is not affecting things, we can verify that the flux tube has statistically homogeneous turbulence in the perpendicular plane.

For example, the observation of localized, steady $E\times B$ flow shear layers (Candy Reference Candy2005; Waltz et al. Reference Waltz, Austin, Burrell and Candy2006; Dominski et al. Reference Dominski, Brunner, Goerler, Jenko, Told and Villard2015; Weikl et al. Reference Weikl, Peeters, Rath, Seiferling, Buchholz, Grosshauser and Strintzi2018; C.J. et al. Reference C.J., Brunner, McMillan, Ball, Dominski and Merlo2020) is a clear indicator that pseudo-rational surfaces may be adversely affecting the accuracy of many present-day simulations. These structures are a general consequence of self-interaction and are a fairly universal in conventional tokamak simulations with kinetic electrons. Although it is difficult to see them in tokamak simulations with adiabatic electrons, figure 6 shows that they can be found by using a slab geometry (instead of toroidal geometry). While they have been identified as an electron-scale phenomenon in linear simulations (Dominski et al. Reference Dominski, Brunner, Goerler, Jenko, Told and Villard2015), in nonlinear simulations they manifest as fairly narrow ion-scale structures. Figure 6 demonstrates that their width does not scale with the electron mass in slab ion-scale simulations and the same behaviour has been seen in previous toroidal ion-scale simulations (Dominski et al. Reference Dominski, McMillan, Brunner, Merlo, Tran and Villard2017). Furthermore, figure 6 shows that similar structures can be found at electron scales in slab electron temperature gradient simulations.

Figure 6.

The time-averaged zonal $E\times B$ shearing rate (i.e. $\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D719}/\unicode[STIX]{x2202}x^{2}$ averaged over time and $y$) as a function of minor radial location within the flux tube for slab simulations (though toroidal simulations with kinetic electrons look similar (Dominski et al. Reference Dominski, Brunner, Goerler, Jenko, Told and Villard2015; C.J. et al. Reference C.J., Brunner, McMillan, Ball, Dominski and Merlo2020)). Data are shown for ion-scale simulations using kinetic electrons with $m_{i}/m_{e}=3672$ (black thick solid), kinetic electrons with $m_{i}/m_{e}=918$ (red thick solid) and adiabatic electrons (grey thick dashed). An electron-scale simulation is shown with adiabatic ions (blue thin solid). The light grey regions indicate the portion significantly affected by pseudo-integer surfaces, while pseudo-half-integer surfaces are located half-way in between. The parameters used are given in tables 1 and 2, except $a/L_{Ts}=10$, $a/L_{ns}=2.0$ and several of the geometry parameters are not needed. Note the reference magnetic field $B_{r}$ and reference length $a$ become arbitrary in slab geometry.

These $E\times B$ flow shear layers are particularly interesting because they are constant in time – i.e. the plasma is spontaneously moving momentum around within the flux tube in order to modify the steady background flow profile. Like large-scale steady $E\times B$ flows, the tokamak will set up a parallel flow in order to convert the $E\times B$ flow layers into a purely toroidal flow (Abel et al. Reference Abel, Plunk, Wang, Barnes, Cowley, Dorland and Schekochihin2013). Such flow represents intrinsic rotation, which is generally prohibited by an underlying symmetry of the gyrokinetic model (Peeters et al. Reference Peeters and Angioni2005; Parra, Barnes & Peeters Reference Parra, Barnes and Peeters2011; Sugama et al. Reference Sugama, Watanabe, Nunami and Nishimura2011). Thus, there must be a symmetry-breaking mechanism present. One plausible mechanism is variation in turbulence characteristics (Parra & Barnes Reference Parra and Barnes2015). This could be introduced by the existence of pseudo-rational surfaces (e.g. the turbulence on integer surfaces is different than on the neighbouring irrational surfaces). Normally, this symmetry-breaking effect is small in $\unicode[STIX]{x1D70C}_{\ast }\ll 1$ because the gradual change in background gradients is the usual cause for the variation in turbulence characteristics. However, the variation caused by the pseudo-integer surfaces occurs on the scale of the gyroradius, so it is not small in $\unicode[STIX]{x1D70C}_{\ast }$. There are only two physical effects that cause variation in turbulence characteristics to break the symmetry (see Section 9 of Parra & Barnes (Reference Parra and Barnes2015)): finite gyroradius effects and the radial magnetic drifts. Since figure 6 shows the shear layers persist in slab geometry, we believe that the finite gyroradius effect is the dominant symmetry-breaking effect. This is intuitive as the symmetry breaking occurs on a small scale and appears to be dipolar (i.e. a region of positive momentum flux is always next to a region of negative), so the radial drift orbits should average over it more effectively because they are significantly larger than the gyroradius.

Importantly, these flow shear layers are just one possible symptom of self-interaction at pseudo-rational surfaces. More generally, we can investigate spurious correlations by looking at the two-point parallel correlation function

(3.2)$$\begin{eqnarray}C_{\Vert }(x,y,\unicode[STIX]{x1D712}_{1},\unicode[STIX]{x1D712}_{2})\equiv \frac{\langle \unicode[STIX]{x1D719}_{NZ}(x,y,\unicode[STIX]{x1D712}_{1},t)\unicode[STIX]{x1D719}_{NZ}(x,y,\unicode[STIX]{x1D712}_{2},t)\rangle _{t}}{\sqrt{\langle \unicode[STIX]{x1D719}_{NZ}^{2}(x,y,\unicode[STIX]{x1D712}_{1},t)\rangle _{t}\langle \unicode[STIX]{x1D719}_{NZ}^{2}(x,y,\unicode[STIX]{x1D712}_{2},t)\rangle _{t}}},\end{eqnarray}$$

where the subscript $NZ$ signifies the non-zonal portion of the quantity and $\left\langle \ldots \right\rangle _{u}$ indicates an average over any coordinate $u$. The quantity $C_{\Vert }$ indicates the degree of correlation between two points $\unicode[STIX]{x1D712}_{1}$ and $\unicode[STIX]{x1D712}_{2}$ on the same field line (i.e. at constant $x$ and $y$). Figure 7 shows the $y$-averaged correlation between adjacent outboard and inboard midplanes $\left\langle C_{\Vert }(x,y,\unicode[STIX]{x1D712}_{1}=0,\unicode[STIX]{x1D712}_{2}=\unicode[STIX]{x03C0})\right\rangle _{y}$ for various toroidal simulations. Note that all simulations with $N_{\text{pol}}=1$ have sharp spikes in the parallel correlation at the locations of pseudo-integer surfaces. Thus, even though toroidal simulations with adiabatic electrons did not display flow shear layers at pseudo-integer surfaces, they do exhibit spikes in parallel correlation. We see that increasing $L_{y}$ increases the spacing between pseudo-integer surfaces, in accordance with (3.1). However, we also notice that neither the width, nor the height of the spikes changes with $L_{y}$. Thus, as we increase $L_{y}$ the regions affected by pseudo-integer surfaces occupy a smaller and smaller fraction of the domain. This indicates one way to eliminate the effect of pseudo-rational surfaces – increase $L_{y}$ until convergence is achieved. If $L_{y}$ is sufficiently large, the pseudo-rational surfaces will have a negligible impact on all volume-averaged quantities. Effectively we are ‘diluting’ away their influence.

Figure 7.

The parallel correlation between inboard and outboard locations following a field line at a given minor radial location within the flux tube. The base simulation has kinetic electrons with $N_{\text{pol}}=1$ and $L_{y}=250\unicode[STIX]{x1D70C}_{i}$ (black thick solid). The other three curves are produced by changing the base simulation to either have adiabatic electrons (grey thick solid), $L_{y}=500\unicode[STIX]{x1D70C}_{i}$ (red thin solid) or $N_{\text{pol}}=2$ (blue thick dashed). Note the $L_{y}=500\unicode[STIX]{x1D70C}_{i}$ simulation has only a single pseudo-integer surface, which occurs at $x=0$. The $N_{\text{pol}}=2$ simulation has no pseudo-integer surfaces, but has pseudo-half-integer surfaces at $x=\{0,\pm 25,\pm 50,\pm 75\}\unicode[STIX]{x1D70C}_{i}$.

Figure 7 also shows a second strategy to eliminate pseudo-integer surfaces – increase $N_{\text{pol}}$ until convergence is achieved. If the flux tube is longer, a turbulent eddy has to remain correlated over a longer distance in order to ‘bite its own tail’. In other words, the $N_{\text{pol}}=2$ has no pseudo-integer surfaces, only pseudo-half-integer surfaces. Accordingly, we see that, compared to the $N_{\text{pol}}=1$ simulation, the $N_{\text{pol}}=2$ case has much smaller peaks at its lowest order pseudo-rational surfaces, which are more closely spaced. This suggests that, given the parameters used for these simulations, a small number of poloidal turns should be sufficient to eliminate the effect of pseudo-rational surfaces and achieve convergence. This is consistent with the original convergence study of Beer et al. (Reference Beer, Cowley and Hammett1995), which found that two or three poloidal turns was sufficient to achieve convergence in gyrofluid simulations with adiabatic electrons. However, it was not obvious that kinetic and adiabatic electron simulations would behave similarly because, in linear simulations, kinetic electrons enable modes that are very extended along the magnetic field line (i.e. ‘giant tails’ in the ballooning envelope) (Hallatschek & Dorland Reference Hallatschek and Dorland2005).

4 Resolution study

We will now use the gyrokinetic code GENE (Jenko et al. Reference Jenko, Dorland, Kotschenreuther and Rogers2000; Görler et al. Reference Görler, Lapillonne, Brunner, Dannert, Jenko, Merz and Told2011) to perform a resolution study in $L_{y}$ and $N_{\text{pol}}$, the two strategies to eliminate pseudo-integer surfaces. Figures 8(a) and 8(b) show the results using standard Cyclone Base Case parameters (Dimits et al. Reference Dimits, Bateman, Beer, Cohen, Dorland, Hammett, Kim, Kinsey, Kotschenreuther and Kritz2000) with adiabatic electrons and kinetic electrons respectively (see table 2). The simulations with adiabatic electrons are considerably less expensive, so we are able to perform a more complete parameter scan and more rigorously ensure adequate resolution (see table 1). Regardless, both display qualitatively similar behaviour.

Figure 8.

Resolution studies in $L_{y}$ and $N_{\text{pol}}$ with the parameters of tables 1 and 2 using adiabatic electrons (a) and kinetic electrons (b). The open black circles indicate simulations with $N_{\text{pol}}=1$ that individually double the grid resolution or the domain width of the three other coordinates (i.e. $x$, $v_{\Vert }$ and $\unicode[STIX]{x1D707}$) as well as the grid resolutions in $y$ and $\unicode[STIX]{x1D712}$. The simulation at $L_{y}=2000\unicode[STIX]{x1D70C}_{i}$ was forced to have double the nominal $x$ domain width (and thus double the number of $x$ grid points) due to the aspect ratio quantization of (2.7). The top axis shows the value of $\unicode[STIX]{x1D70C}_{\ast }$ for which a given $L_{y}$ corresponds to the full flux surface.

For adiabatic electrons, figure 8(a) shows that the fully converged value of the ion heat flux appears to be a bit less than $7$ in gyroBohm units ($Q_{gB}=\unicode[STIX]{x1D70C}_{\ast }^{2}v_{r}n_{i}T_{r}$, where $n_{i}$ is the ion number density). Such convergence can be achieved by simply increasing $L_{y}$ to very large values while maintaining $N_{\text{pol}}=1$. The same result can also be achieved by increasing $N_{\text{pol}}$, which somewhat alleviates our concerns about the validity of using $N_{\text{pol}}>1$ (see § 2). However, this is only true when we are already using a sufficiently large value of $L_{y}$. If $L_{y}$ is not large enough, the simulation will not be fully converged even if $N_{\text{pol}}\rightarrow \infty$. This is for the same reason that we need a sufficiently large domain in the $x$ direction – self-interaction in the perpendicular direction. Completely independent of parallel self-interaction, it is still necessary to ensure that the turbulent correlation length in the $y$ direction is much smaller than the domain width in the $y$ direction (as explained in § 1). If this is not satisfied, convergence will be achieved as $N_{\text{pol}}\rightarrow \infty$ (because the parallel self-interaction will vanish), but it will not converge to the same value as when the binormal self-interaction is also eliminated. Note in figure 8 that the simulations that are not fully resolved in either $L_{y}$ or $N_{\text{pol}}$ have heat fluxes lower than the converged value. This indicates that both binormal and parallel self-interaction tend to reduce the ability of turbulence to transport energy. This may be because the correlations imposed by self-interaction reduce the degrees of freedom accessible to the turbulence, preventing it from behaving in ways that would more effectively drive transport.

Figure 8(a) shows that, using $L_{y}\approx 125\unicode[STIX]{x1D70C}_{i}$ and $N_{\text{pol}}=1$, as is common in the community (Watanabe et al. Reference Watanabe, Sugama, Ishizawa and Nunami2015), under-predicts the ion heat flux by approximately 25 % in adiabatic electron simulations. While this is already concerning, figure 8(b) shows that the effect of self-interaction increases when using kinetic electrons. Instead of a 25 % error, we see roughly 40 %. This is consistent with figure 7, which showed that the kinetic electron simulation had more pronounced spikes in the parallel correlation function. Moreover, the effect of self-interaction can be even bigger for other physical parameters. For example, Watanabe et al. (Reference Watanabe, Sugama, Ishizawa and Nunami2015) shows that taking Cyclone Base Case parameters and lowering the magnetic shear from ${\hat{s}}=0.8$ to $0.2$ increases the error in adiabatic electron simulations from 25 % to at least 50 %. Going further, we simulated Cyclone Base Case parameters with adiabatic electrons and ${\hat{s}}=0$ and found that changing $N_{\text{pol}}=1$ to $N_{\text{pol}}=2$ increased the heat flux by 70 %. Simulations of stellarators with ${\hat{s}}=0$ have even found that turbulence can be stable when using $N_{\text{pol}}=1$, but unstable when $N_{\text{pol}}=2$ (see figure 11 of Faber et al. (Reference Faber, Pueschel, Terry, Hegna and Roman2018)). Clearly, self-interaction through the parallel boundary condition can be a significant effect in local gyrokinetic simulations.

An important practical result is the relative computational cost of the two routes – increasing $N_{\text{pol}}$ or increasing $L_{y}$. This is briefly considered by Beer et al. (Reference Beer, Cowley and Hammett1995) in the context of gyrofluid simulations of a TFTR experimental shot, concluding that ‘it appears that faster convergence is obtained by allowing the domain to be longer than a parallel correlation length’. However, our conclusion from the above Cyclone Base Case simulations is the opposite – increasing $L_{y}$ while maintaining $N_{\text{pol}}=1$ was the cheapest way to reach convergence. For example, the $N_{\text{pol}}=2$, $L_{y}=250\unicode[STIX]{x1D70C}_{i}$ simulation required roughly double the computational resources of the $N_{\text{pol}}=1$, $L_{y}=1000\unicode[STIX]{x1D70C}_{i}$ simulation. This can be understood by studying how the coordinate system grids change in response to increasing $L_{y}$ versus $N_{\text{pol}}$ (see table 1).

We see that as $L_{y}$ is increased the number of $y$ grid points must be increased proportionally. Additionally, after $L_{y}$ reaches a certain value, equation (2.7) forces us to increase $L_{x}$ (and thus the number of $x$ grid points) proportionally as well. However, for the parameters of tables 1 and 2, this was not necessary until $L_{y}>1000\unicode[STIX]{x1D70C}_{i}$. For simulations with lower values of ${\hat{s}}$, for which parallel self-interaction is expected to be more of a concern, this constraint would kick in at lower values of $L_{y}$. On the other hand, when increasing $N_{\text{pol}}$ we clearly must increase the number of parallel grid points. Moreover, as explained in detail in Watanabe et al. (Reference Watanabe, Sugama, Ishizawa and Nunami2015), properly resolving a longer parallel domain turns out to require more radial grid points, due to the fact that longer domains are twisted more by magnetic shear. Adding more radial modes decreases the time step of explicit codes like GENE because we are allowing smaller radial scales, which makes the CFL limit more constraining. In light of this, it is not surprising that simulations with $N_{\text{pol}}>1$ are so costly – it is necessary to increase the parallel, radial and temporal resolutions proportionally with $N_{\text{pol}}$. Additionally, pushing to extreme values of $L_{y}$ has the added advantage of ensuring that self-interaction in the $y$ direction is completely eliminated.

While increasing $L_{y}$ was the better way to reach convergence in our simulations, there are a few important caveats that make generalization difficult. First, at lower values of ${\hat{s}}$ the domain aspect ratio quantization condition (i.e. (2.7)) kicks in at lower values of $L_{y}$, making the $L_{y}$ route more costly. Second, codes that partially use implicit methods to advance in time may observe different computational costs than what we have seen in GENE, which uses only explicit methods, as they are not limited by the CFL condition. Third, adding perpendicular grid points (specifically when increasing $L_{x}$ and $L_{y}$ at fixed resolution) is fundamentally different than adding parallel grid points because the perpendicular modes are all coupled through a Fourier transform in the nonlinear term, while only some of the parallel grid points are connected via the parallel streaming term. Thus, the relative computational and communication costs may depend strongly on the computer architecture, memory and parallelization scheme. Finally, Scott (Reference Scott2001) and Watanabe et al. (Reference Watanabe, Sugama, Ishizawa and Nunami2015) each present a clever way of potentially getting around the need to increase the radial and temporal resolutions with $N_{\text{pol}}$. While neither have been implemented in GENE, they both could make $N_{\text{pol}}>1$ simulations much more affordable, potentially more affordable than large $L_{y}$ simulations. The first approach, the ‘flux tube train’ (Watanabe et al. Reference Watanabe, Sugama, Ishizawa and Nunami2015), has been implemented in the local gyrokinetic codes GKV (Watanabe & Sugama Reference Watanabe and Sugama2006) and stella (Barnes, Parra & Landreman Reference Barnes, Parra and Landreman2019). The other approach, the ‘shifted metric’ coordinate system (Scott Reference Scott2001), has not been implemented in any local gyrokinetic code and it is currently unclear if it is possible to implement in a radially periodic domain (Told & Jenko Reference Told and Jenko2010).

As discussed in § 2, the parallel self-interaction occurring in local simulations is physical when $L_{y}$ corresponds to a full flux surface. The same is true of self-interaction in the $y$ direction – if we are modelling a full flux surface, the toroidal periodicity of a real device will lead to self-interaction if the toroidal correlation length is comparable to the toroidal circumference. Thus, we can also view the $N_{\text{pol}}=1$ points in figure 8 as a scan showing how the physical effect of self-interaction weakens with decreasing $\unicode[STIX]{x1D70C}_{\ast }$Footnote ⁵. This is indicated on the top horizontal axes in figure 8. For context, the DIII-D tokamak (Luxon Reference Luxon2002) has $\unicode[STIX]{x1D70C}_{\ast }\approx 1/300$ and a full flux surface simulation corresponds to

(4.1)$$\begin{eqnarray}L_{y}=\frac{2\unicode[STIX]{x03C0}}{q_{0}}\frac{1}{\unicode[STIX]{x1D70C}_{\ast }}\frac{x_{0}}{a}\unicode[STIX]{x1D70C}_{i}\approx 700\unicode[STIX]{x1D70C}_{i}.\end{eqnarray}$$

Therefore, figure 8 shows that self-interaction is a finite $\unicode[STIX]{x1D70C}_{\ast }$ effect that we are usually justified in eliminating when modelling present-day large tokamaks. However, it might still play an important role in smaller machines or specific parameter regimesFootnote ⁶.

Lastly, we performed a temperature gradient scan to determine how self-interaction affects the critical gradient and profile stiffness in Cyclone Base Case simulations. Figure 9 shows two data sets – one with a large value of $L_{y}$, where self-interaction has been mostly diluted away, and one with a small value of $L_{y}$, where self-interaction is still significant. We see that the difference in the ion heat flux between the two values of $L_{y}$ is generally reduced as the temperature gradient decreases. This leads to similar estimates of the critical gradient for both sets of simulations, giving evidence that large $L_{y}$ domain widths may not be necessary for simulations near the critical gradient. However, the self-interaction present in the simulations with smaller $L_{y}$ reduces the profile stiffness, more significantly when using kinetic electrons.

Figure 9.

A critical gradient study using adiabatic electrons (a) or kinetic electrons (b) for a domain width that mostly eliminates self-interaction (grey triangles) and another for a domain width that allows significant self-interaction (black circles). All simulations use the parameters of tables 1 and 2 and best-fit linear trendlines are also shown.