Eliminating turbulent self-interaction through the parallel boundary condition in local gyrokinetic simulations
Abstract
In this work, we highlight an issue that may reduce the accuracy of many local nonlinear gyrokinetic simulations — turbulent self-interaction through the parallel boundary condition. Given a sufficiently long parallel correlation length, individual turbulent eddies can span the full domain and “bite their own tails,” thereby altering their statistical properties. Such self-interaction is only modeled accurately when the simulation domain corresponds to a full flux surface, otherwise it is artificially strong. For Cyclone Base Case parameters and typical domain sizes, we find that this mechanism modifies the heat flux by roughly 40% and it can be even more important. The effect is largest when using kinetic electrons, low magnetic shear, and strong turbulence drive (i.e. steep background gradients). It is found that parallel self-interaction can be eliminated by increasing the parallel length and/or the binormal width of the simulation domain until convergence is achieved.
pacs:
52.30.Gz, 52.35.Ra, 52.55.Fa, 52.65.Tt1 Introduction
Local nonlinear gyrokinetic simulations are one of the most commonly used tools to assess turbulent transport in tokamaks. They solve the gyrokinetic model [1, 2, 3, 4], a nonlinear system of integro-differential equations that have been derived to model plasma turbulence as accurately as possible. Indeed, the primary approximation of gyrokinetics, (i.e. the tokamak minor radius is much larger than the ion gyroradius ), is satisfied in many existing tokamaks by factors of several hundred. Unfortunately, even with its approximations, such a high-fidelity model is very computationally expensive. For this reason, it has been crucial to streamline the calculations as much as possible.
To this end, it is helpful to minimize the volume of the simulation domain and only include the minimal number of turbulent eddies needed to obtain a statistically-relevant representation of turbulence. Instead of modeling full magnetic flux surfaces across most of the plasma minor radius (referred to as “global” simulations), a calculation can get by with a much smaller domain (referred to as “local” simulations). Since turbulence in tokamaks is aligned with the magnetic field and is very anisotropic, it is important for the domain of a local simulation to have similar characteristics. Thus, the domain is long in the direction parallel to the magnetic field and quite narrow in the two perpendicular directions. Full flux surfaces are not usually modeled, nor is a large fraction of the minor radius. Typically, the domain has a rectangular cross-section on the outboard midplane and is deformed into a parallelogram at other poloidal locations due to the effect of magnetic shear. Such a domain is called a “flux-tube” [5] (see figure 1). It can enable the computational cost of local simulations to be orders of magnitude lower than global simulations, particularly when modeling large devices.
An important difference between local and global simulations concerns boundary conditions. In global simulations, the boundaries within the flux surface are straightforward — because entire flux surfaces are modeled, the physical periodicity in the toroidal and poloidal directions can be implemented directly. However, the radial boundary condition for global simulations is less obvious. Typically Dirichlet boundary conditions are applied together with “buffer regions,” which are radial regions that contain sources of particles and energy [6, 7] ^{1}^{1}1Additionally, an artificial Krook damping operator is often used to prevent profile relaxation in the simulation region.. Hence, it is important to monitor global simulations to ensure that the results are not sensitive to the particularities of this complex and artificial boundary condition.
In contrast, all of the boundary conditions for local simulations are elegant [5], though they still require care (as is emphasized by this work). Since the domain of local simulations is narrow in the radial direction, the background plasma parameters (e.g. density, temperature) are well approximated by a constant value with a linear gradient to drive turbulence. This naturally reflects the inherent scale separation between fluctuations and background gradients, which is a fundamental assumption of the gyrokinetic model [4]. One consequence is that the turbulence experiences the same drive and background plasma conditions at opposing sides of the domain, so it should be statistically identical. This is also true in the binormal direction (i.e. the direction within the flux surface and perpendicular to the magnetic field). For simplicity, instead of generating statistically identical turbulence for these boundaries, it is substituted with exactly identical turbulence by using periodic boundary conditions. Using exactly identical turbulence instead of statistically identical turbulence does have the potential to introduce unphysical correlations into the system. However, as long as the correlation length of the turbulence is much smaller than the domain, this will not occur. In practice, the validity of this can be determined by testing for convergence in the widths (radial and binormal) of the domain.
For nonlinear simulations, the boundary condition in the parallel direction is similar, but with one additional complicating factor — the effect of magnetic shear. In this work we will study the parallel boundary condition and its effect on convergence in the parallel and binormal sizes of the simulation domain. We will find that the domain length typically used by the community, one poloidal turn, can introduce unphysical turbulent correlations, directly affecting the accuracy of the results. Individual turbulent eddies can remain correlated across the entire parallel length of the domain and interact with themselves, which is unphysical unless the domain corresponds to the full flux surface. This was first discussed in the original flux-tube model paper [5] in the context of gyrofluid simulations. Subsequent gyrokinetic studies [8, 9, 10] investigated the ability of this mechanism to generate localized, steady corrugations in the background profiles. More recent work considered the self-interaction of linear eigenmodes and revealed that their parallel self-interaction can even interfere with the numerical convergence of nonlinear simulations with respect to the domain size [11, 12] ^{2}^{2}2Note that the self-interaction considered in these previous works (i.e. the nonlinear interaction of a linear eigenmode with itself) is closely related to the self-interaction considered here (i.e. the interaction of a nonlinear turbulent eddy with itself), but the two are distinct. Most notably, in the absence of magnetic shear, linear eigenmodes would not exhibit self-interaction, but nonlinear eddies still could.. Additionally, it has been explored in gyrokinetic simulations of low magnetic shear stellarators [13, 14]. This paper will show that parallel self-interaction can be eliminated by lengthening the simulation domain in the parallel direction to multiple poloidal turns or making the binormal width of the domain large.
2 The parallel boundary condition
The parallel boundary condition generally used by nonlinear flux-tube simulations is called the “twist-and-shift” condition [5]. 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 , though, due to axisymmetry, they are still free to be at different toroidal angles . This constraint requires the length of the domain to be an integer number of poloidal circuits, which we represent by .
Coordinate | Grid range | Number of gridpoints | ||
Adiabatic | Kinetic | Adiabatic | Kinetic | |
timestep CFL limit |
Parameter | Value | Parameter | Value |
Minor radius of flux-tube, | 0.54 | Major radius, | |
Safety factor, | Magnetic shear, | ||
Temperature gradient, | Density gradient, | ||
Ion-e^{-} mass ratio, | Ion-e^{-} temperature ratio, | ||
Effective ion charge, | Collision frequency, | ||
4^{th} order hyperdiffusion [15], | 4^{th} order hyperdiffusion [15], |
While is free to be any positive integer, it has become overwhelmingly standard in the fusion community to use (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 . For example, such simulations are often not “globally consistent” [17]. This means that, when , the narrow flux-tube can permit Fourier modes that do not exist on the full, doubly-periodic flux surface. However, in the true 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 (see figure 2), where is the toroidal width of the flux-tube. Since is a limited to a finite number of gyroradii in local simulations, the distance separating the poloidal turns is an asymptotically large number of gyroradii. Thus, the particles will perceive the perpendicular direction within the flux surface as infinite to lowest order in the 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 to model a surface), violating global perpendicular periodicity will only introduce an error that is small in . This is intuitive — global consistency should be unimportant in local simulations. Another concern is that flux-tubes with 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 ). 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.
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.
(1) |
where is the safety factor, is the minor radial coordinate within the domain, is the magnetic shear, and the subscript indicates the quantity is evaluated at the center 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 centered on one another [5]. 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 shear^{3}^{3}3For a mathematical explanation of the parallel boundary condition, see reference [5]..
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. Given that we have chosen to center the parallelograms with respect to one another, the flux surface in the center 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 because these surfaces will exhibit the strongest parallel self-interaction. Now if , 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 artifact of the boundary condition and don’t 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 occur midway between the pseudo-integer surfaces.
These pseudo-rational surfaces are also involved in an important constraint that arises from the combination of the twist-and-shift and radial boundary conditions — the aspect ratio of the domain becomes discretized (i.e. the radial domain width divided by the binormal domain width ). To see this, 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 has a matching field line at with which it is equivalent. These field lines are the same and must remain matched as you move in the parallel direction. Otherwise, a particle moving purely along a field line could find itself transported it 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 will always connect them to two other field lines that also form a matching pair. In general, this will only be achieved if the shift in toroidal angle across the domain caused by magnetic shear is equal to the integer times the toroidal domain width, i.e.
(2) |
For the example shown in figures 4 and 5, we see that and 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 , the field lines on the radial boundary would no longer close on themselves, but instead to a different matching pair offset in the direction by .
To see how equation (2) discretizes the domain aspect ratio, we invert the definition of the binormal coordinate
(3) |
(where is a straight-field line poloidal angle) and use equation (1) to calculate the toroidal shift across the sheared domain to be
(4) | ||||
(5) |
Since
(6) |
follows from evaluating equation (3) at constant and , we can use equations (2) and (5) to see that the possible values for the radial domain width are
(7) |
Thus, we find that 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 you cannot decrease below . This causes simulations with very low, but finite magnetic shear to become expensive.
Lastly, it is important to note that the pseudo-rational surfaces created by the parallel boundary condition can be physical [9, 11, 12]. For example, a surface should be modeled using and a surface using . However, for the overall simulation to be physical, the spacing between pseudo-rational surfaces must also be correct. For example, the domain shown in figure 5 has 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 device^{4}^{4}4Or more precisely, twice the distance between the integer surfaces in the device if the safety factor profile was given by equation (1).. 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. To see this, take the change in safety factor across the simulation domain and set it equal to the number of integer surfaces in the simulation . Then, substituting equations (6) and (7) with shows that (i.e. the toroidal domain size is equal to the entire toroidal domain). Thus, we see that modeling rational surfaces accurately requires full flux surface simulations with radial domain widths on the scale of the tokamak minor radius . Boxes of this size offer little computational savings compared to global simulations, but if we try to get by with a smaller domain of , it will have a radial density of pseudo-integer surfaces that is 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 isn’t 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 equation (7), we see that the spacing is given by
(8) |
so we can vary either or . Additionally, we can specifically check the low-order pseudo-rational surfaces to see if they are exhibiting any unusual behavior. 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 flow shear layers [8, 9, 10, 11, 12] 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 [10], 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 behavior has been seen in toroidal ion-scale simulations. Furthermore, figure 6 shows that similar structures can be found at electron-scales in slab Electron Temperature Gradient (ETG) simulations.
These 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 flows, the tokamak will set up a parallel flow in order to convert the flow layers into a purely toroidal flow [4]. Such flow represents intrinsic rotation, which is generally prohibited by an underlying symmetry of the gyrokinetic model [18, 19, 20]. Thus, there must be a symmetry-breaking mechanism present. One plausible mechanism is variation in turbulence characteristics [21]. This could be introduced by the existence of pseudo-rational surfaces (e.g. the turbulence on integer surfaces is different than on the neighboring irrational surfaces). Normally, this symmetry-breaking effect is small in 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 . There are only two physical effects that cause variation in turbulence characteristics to break the symmetry (see Section 9 of reference [21]): 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
(9) |
where the subscript signifies the non-zonal portion of the quantity and indicates an average over any coordinate . The quantity indicates the degree of correlation between two points and on the same field line (i.e. at constant and ). Figure 7 shows the -averaged correlation between adjacent outboard and inboard midplanes for various toroidal simulations. Note that all simulations with 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 a pseudo-integer surfaces, they do exhibit spikes in parallel correlation. We see that increasing increases the spacing between pseudo-integer surfaces, in accordance with equation (8). However, we also notice that neither the width, nor the height of the spikes changes with . Thus, as we increase 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 until convergence is achieved. If 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 also shows a second strategy to eliminate pseudo-integer surfaces — increase 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 has no pseudo-integer surfaces, only pseudo-half-integer surfaces. Accordingly, we see that, compared to the simulation, the 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 reference [5], 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 simulation, kinetic electrons enable modes that are very extended along the magnetic field line (i.e. “giant tails” in the ballooning envelope) [22].
4 Resolution study
We will now perform a resolution study in and , the two strategies to eliminate pseudo-integer surfaces. Figures 8(a) and 8(b) show the results using standard Cyclone Base Case parameters [16] 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 behavior.
For adiabatic electrons, figure 8(a) shows that the fully converged value of the ion heat flux appears to be a bit less than in gyroBohm units (, where is the ion number density). You can achieve such convergence by simply increasing to very large values while maintaining . The same result can also be achieved by increasing , which somewhat alleviates our concerns about the validity of using (see Section 2). However, this is only true when you are already using a sufficiently large value of . If isn’t large enough, the simulation won’t be fully converged even if . This is for the same reason that we need a sufficiently large domain in the direction — self-interaction in the perpendicular direction. Completely independent of parallel self-interaction, you still need to ensure that the turbulent correlation length in the direction is much smaller than the domain width in the direction (as explained in section 1). If this is not satisfied, you will get convergence as (because the parallel self-interaction will vanish), but it won’t converge to the same value as when the binormal self-interaction is also eliminated. Note that both binormal and parallel self-interaction is observed 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 and , as is common in the community [23], under-predicts the ion heat flux by about 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 error, we see roughly . 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, reference [23] shows that taking Cyclone Base Case parameters and lowering the magnetic shear from to increases the error in adiabatic electron simulations from to at least . Going further, we simulated Cyclone Base Case parameters with adiabatic electrons and and found that changing to increased the heat flux by . Simulations of stellarators with have even found that turbulence can be stable when using , but unstable when (see figure 11 of reference [14]). 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 or increasing . This is briefly considered by reference [5] 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 while maintaining was the cheapest way to reach convergence. For example, the , simulation required roughly double the computational resources of the , simulation. This can be understood by studying how the coordinate system grids change in response to increasing versus (see table 1).
We see that as is increased the number of gridpoints must be increased proportionally. Additionally, after reaches a certain value, equation (7) forces us to increase (and thus the number of gridpoints) proportionally as well. However, for the parameters of tables 1 and 2, this was not necessary until . For simulations with lower values of , for which parallel self-interaction is expected to be more of a concern, this constraint would kick in at lower values of . On the other hand, when increasing we clearly must increase the number of parallel gridpoints. Moreover, as explained in detail in reference [23], properly resolving a longer parallel domain turns out to require more radial gridpoints, due to the fact that longer domains are twisted more by magnetic shear. Adding more radial modes decreases the timestep 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 are so costly — you must increase the parallel, radial, and temporal resolutions proportionally with . Additionally, pushing to extreme values of has the added advantage of ensuring that self-interaction in the direction is completely eliminated.
Thus, it seems that increasing is generally the better way to reach convergence. However, there are two important caveats. First, at lower values of the domain aspect ratio quantization condition (i.e. equation (7)) kicks in at lower values of , making the route more costly. Second, references [23] and [24] each present a clever way of potentially getting around the need to increase the radial and temporal resolutions with . While neither have been implemented in GENE, they both could make simulations much more affordable. The first approach, the “flux-tube train” [23], has been implemented in the local gyrokinetic codes GKV [25] and stella [26]. The other approach, the “shifted metric” coordinate system [24], 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 [27].
As discussed in section 2, the parallel self-interaction occurring in local simulations is physical when corresponds to a full flux surface. The same is true of self-interaction in the direction — if you are modeling 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, you can also view the points in figure 8 as a scan showing how the physical effect of self-interaction weakens with decreasing ^{5}^{5}5This convergence looks quantitatively similar to the system size investigations in references [28, 29].. This is indicated on the top horizontal axes in figure 8. For context, the DIII-D tokamak [30] has and a full flux surface simulation corresponds to
(10) |
Therefore, figure 8 shows that self-interaction is a finite effect that we are usually justified in eliminating when modeling present-day large tokamaks. However, it might still play an important role in smaller machines or specific parameter regimes ^{6}^{6}6Note that there are many other finite effects apart from the self-interaction discussed here. Some of these are contained in current global gyrokinetic codes (e.g. profile shearing), while many others have only been derived recently [31] and are not in any code. Thus, while going to a full flux surface simulation will allow you to properly treat self-interaction, it is no guarantee that your overall simulation will be accurate for a small machines. On the contrary, the fact that self-interaction is important suggests that many of the other finite effects may be important too..
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 , where self-interaction has been mostly diluted away, and one with a small value of , where self-interaction is still significant. We see that the difference in the ion heat flux between the two values of 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 domain widths may not be necessary for simulations near the critical gradient. However, the self-interaction present in the simulations with smaller reduces the profile stiffness, more significantly when using kinetic electrons.
5 Conclusions
In this work we have shown that turbulent self-interaction through the parallel boundary condition can significantly affect the results of local gyrokinetic simulations. Such self-interaction can be physical, but only when the simulation domain corresponds to a full flux surface for the device being modeled. Using a narrow flux-tube to model a large device will artificially strengthen self-interaction, which can reduce the heat flux. This implies that any simulation with significant parallel self-interaction is only physically accurate for a single finite value of . In other words, self-interaction is a finite effect.
To achieve the true limit assumed in deriving gyrokinetics, one would like to completely eliminate parallel self-interaction. We have shown in figure 8 that this can be achieved by increasing the width of the simulation domain in the binormal direction and/or lengthening the simulation domain using . Currently available results suggest that self-interaction is stronger for kinetic electrons, low values of magnetic shear, and strong turbulence drive (i.e. steep background gradients). This may be because these parameters tend to increase the parallel correlation length of the turbulence. To verify that self-interaction has been eliminated, one should pay special attention to convergence in the binormal domain width and check for spikes in the parallel correlations function (e.g. figure 7). Additionally, one should take care when using a resolution study done with adiabatic electrons to justify domain resolutions for simulations using kinetic electrons.
References
References
- [1] P.J. Catto. Linearized gyro-kinetics. Plasma Phys., 20(7):719, 1978.
- [2] E.A. Frieman and L. Chen. Nonlinear gyrokinetic equations for low-frequency electromagnetic waves in general plasma equilibria. Phys. Fluids, 25:502, 1982.
- [3] A.J. Brizard and T.S. Hahm. Foundations of nonlinear gyrokinetic theory. Rev. Mod. Phys., 79(2):421, 2007.
- [4] I.G. Abel, G.G. Plunk, E. Wang, M.A. Barnes, S.C. Cowley, W. Dorland, and A.A. Schekochihin. Multiscale gyrokinetics for rotating tokamak plasmas: Fluctuations, transport, and energy flows. Rep. Prog. Phys, 76:116201, 2013.
- [5] MA Beer, SC Cowley, and GW Hammett. Field-aligned coordinates for nonlinear simulations of tokamak turbulence. Phys. Plasmas, 2(7):2687, 1995.
- [6] J. Candy and R.E. Waltz. An Eulerian gyrokinetic-Maxwell solver. J. Comput. Phys., 186(2):545, 2003.
- [7] T. Görler, X. Lapillonne, S. Brunner, T. Dannert, F. Jenko, F. Merz, and D. Told. The global version of the gyrokinetic turbulence code GENE. Journal of Computational Physics, 230(18):7053, 2011.
- [8] J. Candy. Beta scaling of transport in microturbulence simulations. Phys. Plasmas, 12(7):072307, 2005.
- [9] R.E. Waltz, M.E. Austin, K.H. Burrell, and J. Candy. Gyrokinetic simulations of off-axis minimum-q profile corrugations. Phys. Plasmas, 13(5):052301, 2006.
- [10] J. Dominski, S. Brunner, T. Goerler, F. Jenko, D. Told, and L. Villard. How non-adiabatic passing electron layers of linear microinstabilities affect turbulent transport. Phys. Plasmas, 22(6):062303, 2015.
- [11] A. Weikl, A.G. Peeters, F. Rath, F. Seiferling, R. Buchholz, S.R. Grosshauser, and D. Strintzi. The occurrence of staircases in ITG turbulence with kinetic electrons and the zonal flow drive through self-interaction. Phys. Plasmas, 25(7):072305, 2018.
- [12] Ajay C. J., S. Brunner, B. McMillan, J. Dominski, and J. Ball. How the self-interaction mechanism affects zonal flow drive and convergence of flux-tube simulations of turbulent transport. in preparation, 2019.
- [13] M.F. Martin, M. Landreman, P. Xanthopoulos, N.R. Mandell, and W. Dorland. The parallel boundary condition for turbulence simulations in low magnetic shear devices. Plasma Phys. Control. Fusion, 60(9):095008, 2018.
- [14] B.J. Faber, M.J. Pueschel, P.W. Terry, C.C. Hegna, and J.E. Roman. Stellarator microinstabilities and turbulence at low magnetic shear. Journal of Plasma Physics, 84(5), 2018.
- [15] M.J. Pueschel, T. Dannert, and F. Jenko. On the role of numerical dissipation in gyrokinetic vlasov simulations of plasma microturbulence. Computer Physics Communications, 181(8):1428–1437, 2010.
- [16] A.M. Dimits, G. Bateman, M.A. Beer, B.I. Cohen, W. Dorland, G.W. Hammett, C. Kim, J.E. Kinsey, M. Kotschenreuther, A.H. Kritz, and others. Comparisons and physics basis of tokamak transport models and turbulence simulations. Phys. Plasmas, 7:969, 2000.
- [17] B. Scott. Global consistency for thin flux tube treatments of toroidal geometry. Phys. Plasmas, 5(6):2334–2339, 1998.
- [18] A.G. Peeters, C. Angioni, and others. Linear gyrokinetic calculations of toroidal momentum transport in a tokamak due to the ion temperature gradient mode. Phys. Plasmas, 12(7):072515, 2005.
- [19] F.I. Parra, M. Barnes, and A.G. Peeters. Up-down symmetry of the turbulent transport of toroidal angular momentum in tokamaks. Phys. Plasmas, 18(6):062501, 2011.
- [20] H. Sugama, T.H. Watanabe, M. Nunami, and S. Nishimura. Momentum balance and radial electric fields in axisymmetric and nonaxisymmetric toroidal plasmas. Plasma Phys. Control. Fusion, 53(2):024004, 2011.
- [21] F.I. Parra and M. Barnes. Intrinsic rotation in tokamaks: Theory. Plasma Phys. Control. Fusion, 57(4):045002, 2015.
- [22] K. Hallatschek and W. Dorland. Giant electron tails and passing electron pinch effects in tokamak-core turbulence. Phys. Rev. Lett., 95:055002, Jul 2005.
- [23] T.-H. Watanabe, H. Sugama, A. Ishizawa, and M. Nunami. Flux tube train model for local turbulence simulation of toroidal plasmas. Phys. Plasmas, 22(2):022507, 2015.
- [24] B. Scott. Shifted metric procedure for flux tube treatments of toroidal geometry: Avoiding grid deformation. Phys. Plasmas, 8(2):447–458, 2001.
- [25] T.-H. Watanabe and H. Sugama. Velocity–space structures of distribution function in toroidal ion temperature gradient turbulence. Nucl. Fusion, 46(1):24, 2006.
- [26] M. Barnes, F.I. Parra, and M. Landreman. stella: An operator-split, implicit–explicit f-gyrokinetic code for general magnetic field configurations. J. Comput. Phys., 391:365, 2019.
- [27] D. Told and F. Jenko. Applicability of different geometry approaches to simulations of turbulence in highly sheared magnetic fields. Phys. Plasmas, 17(4):042302, 2010.
- [28] B.F. McMillan, X. Lapillonne, S. Brunner, L. Villard, S. Jolliet, A. Bottino, T. Görler, and F. Jenko. System size effects on gyrokinetic turbulence. Phys. Rev. Lett., 105(15):155001, 2010.
- [29] Z. Lin, S. Ethier, T.S. Hahm, and W.M. Tang. Verification of gyrokinetic particle simulation of device size scaling of turbulent transport. Plasma Sci. Technol., 14(12):1125, 2012.
- [30] J.L. Luxon. A design retrospective of the DIII-D tokamak. Nucl. Fusion, 42(5):614, 2002.
- [31] F.I. Parra and I. Calvo. Phase-space Lagrangian derivation of electrostatic gyrokinetics in general geometry. Plasma Phys. Control. Fusion, 53(4):045001, 2011.