1 Introduction
Convective heat loss resulting from micro-turbulent fluctuations in a fusion reactor limits the ability of such a reactor to confine heat to the degree required to achieve fusion. Above some critical temperature gradient, this heat loss rises extremely rapidly, limiting the central temperature and hence the fusion reaction rate. Empirically, one way to improve fusion performance despite this obstacle is to build larger devices.Footnote 1 Current planning for the first demonstration power plant (Federici et al. Reference Federici, Kemp, Ward, Bachmann, Franke, Gonzalez, Lowry, Gadomska, Harman and Meszaros2014) envisions a device more than 20 times larger by volume than the world’s largest operating test reactor, the Joint European Torus (JET). Such a tactic remains the surest route to achieving fusion with our current understanding. However, the approach is not without difficulties. First, the individual cost of such a large plant is very high (although the power output scales accordingly), putting the construction of such reactors, as well as necessary precursor experiments and test reactors, out of the reach of all but major governments. Second, the larger designs place much greater demands upon the construction materials (see e.g. the discussion in Stork et al. (Reference Stork, Agostini, Boutard, Buckthorpe, Diegele, Dudarev, English, Federici, Gilbert and Gonzalez2014)).
The reactor design effort is constantly finding innovative ways to mitigate these problems. For example, one approach to achieving more compact and cost-effective designs is to use new high-temperature superconductors which have lower cost and greater tensile strength than the superconductors used today, and which can be constructed so as to allow easy disassembly for maintenance (Sorbom et al. Reference Sorbom, Ball, Palmer, Mangiarotti, Sierchio, Bonoli, Kasten, Sutherland, Barnard and Haakonsen2015).
However, while such approaches have frequently used state-of-the-art modelling for almost all of the construction of the reactor, none of these studies have incorporated the turbulent heat loss using nonlinear models which accurately capture the essential properties of the turbulence. Instead, the vast majority use scaling laws, either inspired by dimensional arguments or extrapolated from current experiments (Uckan Reference Uckan1990; Galambos et al. Reference Galambos, Perkins, Haney and Mandrekas1995; Luce et al. Reference Luce, Challis, Ide, Joffrin, Kamada, Politzer, Schweinzer, Sips, Stober and Giruzzi2014; Sorbom et al. Reference Sorbom, Ball, Palmer, Mangiarotti, Sierchio, Bonoli, Kasten, Sutherland, Barnard and Haakonsen2015) (see also discussion in Zohm et al. (Reference Zohm, Angioni, Fable, Federici, Gantenbein, Hartmann, Lackner, Poli, Porte and Sauter2013)). It is also possible to use more sophisticated ‘quasilinear’ models (Staebler, Kinsey & Waltz Reference Staebler, Kinsey and Waltz2007; Bourdelle et al. Reference Bourdelle, Citrin, Baiocchi, Casati, Cottier, Garbet and Imbeaux2016) which calculate the linear properties of the instabilities which drive the turbulence and use phenomenological rules, and calibration to a pre-determined subset of nonlinear simulations, to calculate the turbulent heat loss from these linear properties. Use of these quasilinear models includes both full reactor design studies (e.g. Jardin et al. (Reference Jardin, Kessel, Mau, Miller, Najmabadi, Chan, Chu, Lahaye, Lao and Petrie2006), Staebler & John (Reference Staebler and John2006), Wenninger et al. (Reference Wenninger, Arbeiter, Aubert, Aho-Mantila, Albanese, Ambrosino, Angioni, Artaud, Bernert and Fable2015)) and efforts to predict and optimise performance for particular devices (e.g. Mukhovatov et al. (Reference Mukhovatov, Shimomura, Polevoi, Shimada, Sugihara, Bateman, Cordey, Kardaun, Pereverzev and Voitsekhovitch2003), Budny (Reference Budny2009), Kinsey et al. (Reference Kinsey, Staebler, Candy, Waltz and Budny2011), Parail et al. (Reference Parail, Albanese, Ambrosino, Artaud, Besseghir, Cavinato, Corrigan, Garcia, Garzotti and Gribov2013), Meneghini et al. (Reference Meneghini, Snyder, Smith, Candy, Staebler, Belli, Lao, Park, Green and Elwasif2016)). Additionally, there have been advances in using neural-network approaches to allow such methods to fully replace the simple scaling laws by facilitating extremely fast simulation based on the quasilinear models (Citrin et al. Reference Citrin, Breton, Felici, Imbeaux, Aniel, Artaud, Baiocchi, Bourdelle, Camenen and Garcia2015a ). Nonetheless, nonlinear turbulence calculations, while used routinely for investigation of individual experiments and verification of quasilinear models, have been conspicuously absent from reactor design.Footnote 2
The primary reason for the absence of nonlinear models is their cost and complexity. Effectively, it would have been impossible until recently to have completed a design study including these models on any reasonable time scale. This absence represents a missed opportunity. Over the last 20 years there has been a great increase in understanding, derived from both experiment and theoretical inquiry, of the circumstances in which the levels of turbulence can be greatly reduced, without reducing the temperature gradient (Burrell Reference Burrell1997; Synakowski Reference Synakowski1999; Highcock et al. Reference Highcock, Barnes, Schekochihin, Parra, Roach and Cowley2010; Barnes et al. Reference Barnes, Parra, Highcock, Schekochihin, Cowley and Roach2011a ; Highcock et al. Reference Highcock, Schekochihin, Cowley, Barnes, Parra, Roach and Dorland2012; Citrin et al. Reference Citrin, Garcia, Görler, Jenko, Mantica, Told, Bourdelle, Hatch, Hogeweij and Johnson2015b ). Such effects can be, and have been, included in quasilinear and other models a posteriori. However, even where designs are based on extrapolation from the best performing experimental configurations, using the most complete reduced models available at the time (e.g. Jardin et al. (Reference Jardin, Kessel, Mau, Miller, Najmabadi, Chan, Chu, Lahaye, Lao and Petrie2006)), the designs can know nothing of potential new phenomena, leading to possible further gains, that might exist in the vast configuration space at their disposal (more than two decades ago, Galambos et al. (Reference Galambos, Perkins, Haney and Mandrekas1995) showed that in theory, where full control could be maintained over turbulent transport, it was possible to reduce the capital cost of the reactor plant by 30 %, that is, 10 billion (1995) US dollars, compared to the most advanced designs of the day). This still largely uncharted configuration space is best explored in conjunction with nonlinear models, since any exploration by reduced models could conceivably miss a crucial physics effect not yet included within the modelFootnote 3 (cf. work in the stellarator community which clearly demonstrates the benefit of tight integration of stellarator optimisation with nonlinear analysis; Xanthopoulos et al. Reference Xanthopoulos, Mynick, Helander, Turkin, Plunk, Jenko, Görler, Told, Bird and Proll2014).
Every fusion experiment costs a vast sum to build, and owing to the hard limits placed by all components of the reactor, can only explore a certain small region of the high- (or infinite-) dimensional design space. It is beholden upon theory and numerical study to search the design space, at several-orders-of-magnitude smaller expense, to identify the most promising regions. However, in order for this search to be meaningful, we must have confidence in the predictions. The simplified models used must either be pessimistic, or be used only in regimes where they have been benchmarked against more complete models or experiment, and thus preclude the possibility of entirely new design paradigms. By contrast, careful comparison with experiment over recent years has provided a high degree of confidence that nonlinear models (specifically gyrokinetic models, which use a five-dimensional reduction of the Vlasov equation valid in the conditions of a fusion reactor; for a review see Abel et al. (Reference Abel, Plunk, Wang, Barnes, Cowley, Dorland and Schekochihin2013)), without any provision of tuning or fitting parameters, are able to accurately predict what the properties of turbulence will be like in a given situation (White et al. Reference White, Howard, Greenwald, Reinke, Sung, Baek, Barnes, Candy, Dominguez and Ernst2013; Citrin et al. Reference Citrin, Jenko, Mantica, Told, Bourdelle, Dumont, Garcia, Haverkort, Hogeweij and Johnson2014; van Wyk et al. Reference van Wyk, Highcock, Schekochihin, Roach, Field and Dorland2016).
2 A first-principles model
Here we present the first case where a first-principles nonlinear model of turbulence (specifically a novel hybrid gyrofluid/gyrokinetic model, described below, which produces excellent agreement with gyrokinetic models) has been used, as part of a multi-scale transport analysis, to calculate the performance of the core of a given configuration ab initio and then seek for a higher performing solution. In this particular study, the optimisation algorithm chosen achieved an improvement in the fusion power per unit volume of 91 %. However, this first study is envisioned as a proof of concept, as the stepping stone for a much larger effort in which many more dimensions of parameter space are explored.
An illustration of the nested magnetic flux surfaces that confine a fusion plasma. The shape of these fluxes corresponds to the initial shape used for this study.
A magnetic fusion reactor is composed of a plasma confined by a series of nested toroidal magnetic flux surfaces (figure 1). The plasma is composed of ionised deuterium and tritium, which fuse together to produce helium, neutrons and a large amount of energy, provided sufficient pressure is achieved over sufficient volume in the centre of the plasma. The magnetic field is provided by a series of large external field coils and a current which flows through the plasma itself (as well as by various small coils responsible in addition for the stability of the field). The plasma is heated in a number of different ways: by the current itself, by the injection of a beam of high energy neutral particles and by high power electromagnetic radiation of various frequencies. Fuel is injected by means of puffs of gas, frozen pellets and neutral particle beams. Heat is lost via neutrons, radiation and by particle loss to the wall (which is concentrated in a special target region called the divertor). In this investigation we consider the core of the reactor; that is we consider the inner
${\sim}90\,\%$
(by minor radius) of the flux surfaces. We assume the existence of an edge transport barrier, that is, that there is good confinement with a steep pressure gradient in the remaining
${\sim}10\,\%$
of the plasma (a standard assumption since the discovery of the ‘high-confinement mode’; Wagner et al.
Reference Wagner, Fussmann, Grave, Keilhacker, Kornherr, Lackner, McCormick, Müller, Stäbler and Becker1984). We also assume that the barrier is independent of elongation and triangularity (this assumption, necessary because we are considering only core transport, differs from current experimental observations, and is discussed further below).
Schematic of the CoRFu framework.
The software used in this study to optimise such a reactor is the CodeRunner Fusion Optimisation Framework (CoRFu). This is a sizeable edifice of software comprising many independent parts, as illustrated in figure 2. The equations underlying the framework are given in appendices A–D. The calculation proceeds as follows.
To begin with, the choice of optimisation control parameters, and all other requisite input information, is provided to the CoRFu framework, which acts as the optimisation driver. From this input, an initial configuration is assembled, comprising parameterisations of:
(i) the shape of the flux surface at the edge of the plasma core;
(ii) the profile of the toroidal current; and
(iii) the initial pressure profile, including a fixed finite pressure at the edge of the plasma core.
These are used (link 1 in figure 2) by the Chease code (Lütjens, Bondeson & Sauter Reference Lütjens, Bondeson and Sauter1996) to calculate the shape of the magnetic flux surfaces (i.e. the equilibrium magnetic field).
For the initial configuration we consider a purely hypothetical tokamak. We do, however, choose machine parameters (see table 1) that are comparable in magnitude to those of a large tokamak – the Joint European Torus (JET) – and profiles that are plausible for such a tokamak (figure 5; see e.g. Barnes et al. (Reference Barnes, Parra, Highcock, Schekochihin, Cowley and Roach2011a ) or the tokamak profile database described in Roach et al. (Reference Roach, Walters, Budny, Imbeaux, Fredian, Greenwald, Stillerman, Alexander, Carlsson and Cary2008) for a comparison).
The initial configuration, including the magnetic equilibrium solution, details of external heat sources and the profiles of particle density (densities and external sources are held fixed) are then used (links 2 and 3) by the Trinity transport solver (Barnes et al. Reference Barnes, Abel, Dorland, Görler, Hammett and Jenko2010) to find the pressure (and thus the fusion power generated), across the whole of the plasma, with the additional assumption that the temperatures of the ions and electrons are kept equal by collisional equilibration.Footnote 4 This is done by allowing the pressure profile to evolve until the heat losses, due to turbulence (discussed below) and other effects (in this study, neoclassical transport and radiative losses, see appendix A), match the heat inputs, both external and that generated by fusion. However, since the magnetic equilibrium itself is affected by the pressure, the result of the Trinity calculation must be fed back to Chease (link 4), a new equilibrium generated and Trinity rerun, until the cycle converges, the pressure stops changing and a steady-state solution is found.
An illustration of a Trinity simulation showing above-threshold turbulence for the optimal configuration. The elongated ribbons to the left of the main picture are flux tubes (Beer, Cowley & Hammett Reference Beer, Cowley and Hammett1995), the domains in which the turbulence is calculated, using the local approximation which assumes the turbulent fluctuations are small-scale perpendicular to the magnetic field. Assuming statistical periodicity (which again results from the small scale of the fluctuations), these flux tubes can be repeated to fill the four flux surfaces being simulated (which are shown cut away to the right of the main picture). The picture to the right shows the cross-section of the turbulence.
Global properties of the initial, intermediate and optimal configurations.
To calculate the heat loss due to turbulence (links 5 and 8), Trinity uses several simultaneous copies of the GPU-based hybrid gyrofluid/gyrokinetic code Gryfx (see appendices B–D). Each copy of Gryfx calculates the heat flux resulting from the turbulence at a particular location in the plasma, that is, on a particular magnetic flux surface, given the pressure profile and magnetic equilibrium. A mean-field, multi-scale approach is used. The turbulent heat flux is calculated for the current pressure profile; the heat flux is then used to evolve the pressure profile, at which point the turbulent heat flux is recalculated. An illustration of the complete system is given in figure 3.
The Gryfx code divides the calculation of the heat flux in two (links 6 and 7) using a novel algorithm: the evolution of the smaller scale drift waves is calculated using the gyrofluid model (Dorland & Hammett Reference Dorland and Hammett1993; Beer & Hammett Reference Beer and Hammett1996) by Gryfx, and the evolution of the large-scale zonal flows which are generated by the turbulence is calculated by the gyrokinetic code Gs2 (Kotschenreuther, Rewoldt & Tang Reference Kotschenreuther, Rewoldt and Tang1995; Dorland et al. Reference Dorland, Jenko, Kotschenreuther and Rogers2000), with the nonlinear interaction between the two scales (and among the drift waves of various scales) being calculated in Gryfx (see appendix D). The simulations are electrostatic with a Boltzmann response for electrons.
The transport calculation within Trinity used only five radial grid points. This was to reduce the overall cost of this first study, and would be increased in future work. Nonetheless, we note that the 5-point derivative stencil used radially by Trinity, and the fact that the transport calculations evolved to steady state (where power balance was satisfied on each flux surface), will reduce the impact of the sparse radial grid on the results of this study.
3 Finding the optimal solution
When the pressure, magnetic equilibrium and turbulent heat fluxes have ceased to evolve, a steady-state solution has been found. At this point, the figure of merit for this solution (the fusion power per unit volume) is used by the optimisation driver CodeRunner (link 9 in figure 2), to generate a new set of control parameters according to the particular optimisation algorithm being used. In this study, the Nelder–Mead simplex algorithm was selectedFootnote
5
(not to be confused with the simplex algorithm in linear programming). This is a slow but robust algorithm that constructs a simplex with
$N+1$
vertices in an
$N$
-dimensional configuration space and then at each stage takes the ‘worst’ vertex and moves it to a more optimal place. Visually, one can perceive the simplex ‘crawling uphill’ (figure 4).
(a) Diagram illustrating the search path taken by the simplex optimisation algorithm. The suspended triangles represent the simplex at each iteration, with the colour and height of each triangle equal to the average of the values of
$P/V$
at each of the vertices. The curve in the base plane indicates the actual sequence of function calls, that is, the actual sequence of (
$\unicode[STIX]{x1D705}$
,
$\unicode[STIX]{x1D6FF}$
) values evaluated by Trinity. (b)
$P/V$
as a function of
$\unicode[STIX]{x1D705}$
and
$\unicode[STIX]{x1D6FF}$
, both the function values (squares) and an interpolated surface, with the colour being the value of
$P/V$
. The highest (optimal) value is at
$\unicode[STIX]{x1D705}=1.60$
,
$\unicode[STIX]{x1D6FF}=-0.625$
. An intermediate solution, with a triangularity within the bounds of what has currently been realised, is chosen at
$\unicode[STIX]{x1D705}=1.525$
,
$\unicode[STIX]{x1D6FF}=-0.25$
.
Using the new values of the control parameters provided by the simplex algorithm (i.e. the next ‘guess’) the whole of the previous cycle is repeated, determining the figure of merit for this next configuration. The calculation proceeds likewise until it has converged on a maximum to within some specified tolerance (or until a prescribed number of iterations have taken place).
(a) Profiles of ion temperature (
$T_{i}$
), electron density (
$n_{e}$
), safety factor (
$q$
) and heat input (
$S_{i}$
), showing the values of each quantity versus normalised radius
$\unicode[STIX]{x1D70C}$
.
$\unicode[STIX]{x1D70C}$
is a dimensionless quantity which labels the flux surface. In this work we examine core transport; thus at the outer flux surface (
$\unicode[STIX]{x1D70C}=1$
) there is a finite temperature of 3.2 keV and a finite electron density of
$0.6\times 10^{20}~\text{m}^{-3}$
. (b) Outer flux surface shapes of the initial, intermediate and optimal solutions.
The two control parameters chosen were the elongation and triangularity of the outer flux surface (as defined in Lütjens et al. (Reference Lütjens, Bondeson and Sauter1996)). The shape of the magnetic field has long been known to have a marked effect on the turbulence, and the choice of these two shaping parameters in particular was motivated by pioneering experimental work on the Tokamak à Configuration Variable (TCV), which was designed to allow large variation in both – demonstrating a consequent large variation in performance (Weisen et al. Reference Weisen, Alberti, Berry, Behn, Blanchard, Bosshard, Bühlmann, Chavan, Coda and Deschenaux1999).
The results are displayed at the bottom of figure 4. In this figure one can observe the evolution of the control parameters, starting with three initial guesses of
$(\unicode[STIX]{x1D705},\unicode[STIX]{x1D6FF})=(1.3,0.2)$
, (1.3, 0.1) and (1.4, 0.1). The lowest initial vertex is (1.3,0.2) with a power per unit volume (
$P/V$
) of
$0.0313~\text{MW}\,\text{m}^{-3}$
. Initially the algorithm moves to higher elongation and lower triangularity, before proceeding to keep elongation roughly constant and to continue to decrease triangularity. The iteration was terminated when values of triangularity moved significantly beyond the limits of what has been seen in experiment (
${\sim}-0.65$
(Pochelon et al.
Reference Pochelon, Goodman, Henderson, Angioni, Behn, Coda, Hofmann, Hogge, Kirneva and Martynov1999)). The final, optimal value of
$P/V$
of
$0.0598~\text{MW}\,\text{m}^{-3}$
with
$(\unicode[STIX]{x1D705},\unicode[STIX]{x1D6FF})=(1.6,-0.625)$
. Thus, an improvement of 91 % was discovered over the course of the optimisation. Since the final value of triangularity was somewhat extreme, we also provide data for an intermediate result with
$(\unicode[STIX]{x1D705},\unicode[STIX]{x1D6FF})=(1.525,-0.25)$
, for which
$P/V$
is
$0.0481~\text{MW}\,\text{m}^{-3}$
.
It is interesting also to consider the evolution of the confinement time, defined as the total amount of stored kinetic energy divided by the rate of energy injection, is displayed in figure 6. We see that optimising the confinement time, which reflects the ability of the system to confine heat, would have produced different results to optimising the fusion power per unit volume; in particular, the improvement at negative triangularity is less marked, and moving to higher elongation at positive triangularity may have produced a similar improvement in the confinement time.
Evolution of the confinement time in seconds, given as a function of elongation
$\unicode[STIX]{x1D705}$
and triangularity
$\unicode[STIX]{x1D6FF}$
. Squares represent data points, shown on top of an interpolated surface for illustration of the overall trend.
4 Discussion and outlook
It has been shown that, without any tuning of model parameters, it has been possible to generate a plausible solution for reactor performance and then optimise it. The final fusion power in this study is small compared to the input heating power (on the ions) of 20 MW, but that is because the parameters in this study (table 1) are not reactor-like, but rather of the same order of magnitude as those of JET, the largest existing device: a choice made to allow assessment of the credibility of the initial solution. Even more encouragingly, the qualitative result of this optimisation – that core confinement improved at lower triangularity and higher elongation – is directly corroborated by experimental measurements on the TCV tokamak (Weisen et al. Reference Weisen, Alberti, Berry, Behn, Blanchard, Bosshard, Bühlmann, Chavan, Coda and Deschenaux1999) and theoretical studies (Marinoni et al. Reference Marinoni, Brunner, Camenen, Coda, Graves, Lapillonne, Pochelon, Sauter and Villard2009), although as TCV is a much smaller device than that we consider here, with important effects arising from kinetic electron modes and finite gyroradius, quantitative comparisons cannot be justified.
The challenges inherent to this study, and the obstacles to such a study being carried out in the past, proceeded principally from the calculation of the turbulent fluxes. All in all a total of 8680 converged nonlinear turbulence calculations at a spatial resolution of
$192\times 96\times 20$
(radial
$\times$
poloidal
$\times$
parallel) were required, at a total cost of approximately 3000 GPU-accelerated node hours. That such an exercise was possible was owing to two developments. The first was the design of the implicit algorithm in the Trinity code which reduced the number of turbulence calculations required (Barnes et al.
Reference Barnes, Abel, Dorland, Görler, Hammett and Jenko2010). The second was the development of Gryfx, which by using a gyrokinetic response for the zonal flows and advanced nonlinear closures was able to overcome the shortcomings of gyrofluid codes in the 1990s (Dimits et al.
Reference Dimits, Bateman, Beer, Cohen, Dorland, Hammett, Kim, Kinsey, Kotschenreuther and Kritz2000) while still maintaining their huge speed advantage over gyrokinetic codes (typically a GRYFX run will be 20 times faster and 100 times less expensive than, for example, a Gs2 run; Mandell & Dorland Reference Mandell and Dorland2014).
It is only right to point out in this discussion that this study focuses solely on the plasma core, and not on the reactor components or the plasma edge. It should be noted that though edge transport barriers can be achieved in combination with negative triangularity (Pochelon et al. Reference Pochelon, Angelino, Behn, Brunner, Coda, Kirneva, Medvedev, Reimerdes, Rossel and Sauter2012), the evolution of the temperature at the edge of the domain, here held constant, would be strongly altered if the edge of the plasma was modelled in addition to the core, and thus the problem of optimising the performance of the plasma as whole is distinct from optimising that of the core, as is discussed further in Pochelon et al. (Reference Pochelon, Angelino, Behn, Brunner, Coda, Kirneva, Medvedev, Reimerdes, Rossel and Sauter2012). In fact it is likely that the edge confinement would degrade as the triangularity becomes negative, reducing the favourability of negative triangularity regimes (see e.g. Merle, Sauter & Yu Medvedev (Reference Merle, Sauter and Yu Medvedev2017)). Thus further work is needed to understand how these results would be modified if the plasma edge was modelled self-consistently with the core. In light of these limitations, it would be most desirable in the future to incorporate this methodology in a holistic model of the device (including checks on magneto-hydrodynamic stability); in particular it would be desirable to control the safety factor profile to maintain values greater than one in the core to avoid the sawtooth instability (the small increase in the safety factor required is not expected to significantly affect the results presented here). However, while ultimately all parts of a tokamak must be considered, by finding ways to improve core performance the demands placed upon the reactor design, and the requirement for high plasma confinement in the edge, can be lessened.
Having demonstrated the first successful optimisation of (core) confinement using a nonlinear model of turbulence, there are several directions that are immediately attractive to follow. The first is to use this technique to find optimal configurations for JET ahead of its first run using an active (deuterium–tritium) fuel mix in twenty years (using parameters that match JET far more closely than in this study). The second is to consider, in a similar fashion, ways in which the performance of ITER, the new, global fusion experiment being constructed, can be improved. It would also be important to extend Gryfx to consider electromagnetic effects and kinetic electron physics, and to include the effects of rotation, which are not considered here. Other priorities would be to include additional transport channels: electron and ion heat and particles, momentum, impurity transport and so on. The complex interplay between these channels can be the crucial factor in achieving higher performance (see e.g. Highcock et al. (Reference Highcock, Barnes, Parra, Schekochihin, Roach and Cowley2011)), and robust operation, and is an ideal area to study with a nonlinear framework such as CoRFu. The eventual goal is to switch to using optimisation algorithms which parallelise easily (Lewis Reference Lewis2004) and can search for global maxima, and run a vast parallel blue-skies search for a dramatically optimised fusion reactor.
Acknowledgements
The authors of this paper are indebted to A. Schekochihin, T. Fülöp, G. Hammett, J. Ball, J. Citrin, M. Landremann, I. Abel, F. van Wyk, F. Parra, C. Roach, S. Newton, I. Pusztai and J. Omotani for helpful discussions, ideas and support. This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014–2018 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. This work has also been supported by the Framework grant for Strategic Energy Research (Dnr. 2014–5392) from Vetenskapsrådet. N.R.M. is supported by the US Department of Energy via the DOE CSGF program, provided under grant DE-FG02-97ER25308. This work used the TACC Stampede supercomputer and the ARCHER UK National Supercomputing Service (http://www.archer.ac.uk).
Appendix A. Equations underlying the optimisation framework
The equations that govern the system are presented in Abel et al. (Reference Abel, Plunk, Wang, Barnes, Cowley, Dorland and Schekochihin2013), which builds upon an already large body of work (significantly Frieman & Chen Reference Frieman and Chen1982; Sugama & Horton Reference Sugama and Horton1998).
They follow a multi-scale approach, with a clear separation in time between the evolution of the safety factor (the twist of the magnetic field lines, which evolves on the resistive time scale; Abel & Cowley Reference Abel and Cowley2013), the evolution of the pressure (which evolves on the confinement time scale) and the evolution of the turbulence (which occurs on the time scale of the linear micro-instabilities which drive it). The equations governing all three are
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}q}{\unicode[STIX]{x2202}t}=\frac{c}{4\unicode[STIX]{x03C0}^{2}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}V^{\prime }\left\langle \boldsymbol{E}\boldsymbol{\cdot }\boldsymbol{B}\right\rangle _{\unicode[STIX]{x1D713}} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{3}{2}\frac{1}{V^{\prime }}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}V^{\prime }\left\langle n_{s}\right\rangle _{\unicode[STIX]{x1D713}}T_{s}+\frac{1}{V^{\prime }}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}V^{\prime }\left\langle Q_{s}\right\rangle _{\unicode[STIX]{x1D713}}=S_{s} & \displaystyle\end{eqnarray}$$
and
$$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle \frac{\unicode[STIX]{x2202}h_{s}}{\unicode[STIX]{x2202}t}+\left(v_{\Vert }\hat{\boldsymbol{b}}+\boldsymbol{V}_{D_{s}}+\frac{c}{B}\hat{\boldsymbol{b}}\times \unicode[STIX]{x1D735}\left\langle \unicode[STIX]{x1D711}\right\rangle _{\mathbf{R}}\right)\boldsymbol{\cdot }\unicode[STIX]{x1D735}h_{s}\nonumber\\ \displaystyle & & \displaystyle \displaystyle \quad =\left\langle C\left[h_{s}\right]\right\rangle _{\mathbf{R}}+\frac{Z_{s}eF_{s}}{T_{s}}\frac{\unicode[STIX]{x2202}\left\langle \unicode[STIX]{x1D711}\right\rangle _{\mathbf{R}}}{\unicode[STIX]{x2202}t}-\frac{\unicode[STIX]{x2202}F_{s}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}\left(\hat{\boldsymbol{b}}\times \unicode[STIX]{x1D735}\left\langle \unicode[STIX]{x1D711}\right\rangle _{\mathbf{R}}\right)\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713},\end{eqnarray}$$
respectively. In these equations,
$q$
is the magnetic safety factor,
$\unicode[STIX]{x1D713}$
is the poloidal magnetic flux which is contained within a flux surface,
$V$
is the volume of the flux surface,
$V^{\prime }=\text{d}V/\text{d}\unicode[STIX]{x1D713}$
is the incremental volume element (loosely, the flux surface area),
$\boldsymbol{E}$
is the electric field and
$\boldsymbol{B}$
the magnetic field, and
$\left\langle \right\rangle _{\unicode[STIX]{x1D713}}$
denotes an average over the flux surface. The index
$s$
labels the charged particle species (e.g. electron, deuterium ion, tritium ion etc.),
$T_{s}$
and
$n_{s}$
are the species temperature and density,
$Q_{s}$
is the flux of heat across a flux surface, and
$S_{s}$
is a volumetric heat source. The quantity
$h_{s}$
is the fluctuating (turbulent) part of the particle distribution function,
$v_{\Vert }$
is the velocity along the field line,
$\hat{\boldsymbol{b}}=\boldsymbol{B}/\left|\boldsymbol{B}\right|$
,
$\boldsymbol{V}_{D_{s}}$
represents the effect of the magnetic field inhomogeneities on particle motion,
$\unicode[STIX]{x1D711}$
is the perturbed electric potential,
$C$
represents the effects of collisions between particles,
$Z_{s}e$
is the species charge and
$F_{s}$
is the equilibrium (slowly varying) component of the species distribution function. The operator
$\left\langle \right\rangle _{\mathbf{R}}$
denotes an average over gyrophase at fixed gyrocentre location
$\boldsymbol{R}$
. The turbulent component of the heat flux
$Q_{s}$
can be trivially obtained from the turbulent distribution function
$h_{s}$
$$\begin{eqnarray}Q_{s}=\left\langle \int \text{d}^{3}\boldsymbol{v}\frac{m_{s}v^{2}}{2}\left[\left(\frac{c}{B}\hat{\boldsymbol{b}}\times \unicode[STIX]{x1D735}\left\langle \unicode[STIX]{x1D711}\right\rangle _{\mathbf{R}}\right)\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\right]\left\langle h_{s}\right\rangle _{\mathbf{r}}\right\rangle _{\unicode[STIX]{x1D713},t}\end{eqnarray}$$
where
$\left\langle \right\rangle _{\mathbf{r}}$
denotes a gyroaverage at constant particle position. The (usually) subdominant neoclassical component of the heat flux is modelled using analytical approximations (Chang & Hinton Reference Chang and Hinton1982), as are radiative losses due to bremsstrahlung (Glasstone & Lovberg Reference Glasstone and Lovberg1960).
Since we are seeking a steady-state solution, equation (A 1) is not evolved during this study. Instead we prescribe a fixed profile of surface averaged toroidal current. This causes
$q$
to change with the pressure gradient, not on the resistive time scale, during the evolution towards steady state, but once steady state is reached
$q$
ceases to evolve and the solution is self-consistent. Using the prescribed profile of toroidal current, a prescribed outer flux surface and the pressure profile, the poloidal flux (and hence the shape of the magnetic flux surfaces) can be determined using the Grad–Shafranov equation:
$$\begin{eqnarray}\unicode[STIX]{x1D6E5}^{\ast }\unicode[STIX]{x1D713}=-4\unicode[STIX]{x03C0}R^{2}\mathop{\sum }_{s}n_{s}T_{s}\left\{\frac{\text{d}\ln n_{s}}{\text{d}\unicode[STIX]{x1D713}}+\frac{\text{d}\ln T_{s}}{\text{d}\unicode[STIX]{x1D713}}\right\}-I(\unicode[STIX]{x1D713})\frac{\text{d}I}{\text{d}\unicode[STIX]{x1D713}},\end{eqnarray}$$
where
$\unicode[STIX]{x1D6E5}^{\ast }$
is the Grad–Shafranov operator
$$\begin{eqnarray}\unicode[STIX]{x1D6E5}^{\ast }\unicode[STIX]{x1D713}=\left(\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}R^{2}}-\frac{1}{R}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}R}+\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}Z^{2}}\right)\unicode[STIX]{x1D713}.\end{eqnarray}$$
Note that when using the Chease code (Lütjens et al.
Reference Lütjens, Bondeson and Sauter1996) to solve the Grad–Shafranov equation, the usual function
$I\unicode[STIX]{x2202}I/\text{d}\unicode[STIX]{x1D713}$
can be replaced as an input by the surface averaged toroidal current density, which is what is specified in this study.
The equation for the pressure (A 2) is evolved by Trinity (Barnes et al.
Reference Barnes, Abel, Dorland, Görler, Hammett and Jenko2010; Trinity is also capable of evolving the density and rotation, which are kept fixed for this study). This equation determines how the evolution of the pressure is governed by the total flux of heat, given sources and the shape of the magnetic flux surfaces (that is, the solution of the Grad–Shafranov equation). Since the solution of the Grad–Shafranov equation itself varies on the same time scale as the pressure, the solution is periodically updated with the new pressure profile until steady state is reached. The pressure equation can be evolved separately for each charged species; in this study, to reduce expense, it is solved for the deuterium ions only. It is then assumed that collisional processes will rapidly equilibrate the temperatures of all species, and hence that the temperatures of the electrons and tritium ions is the same as that of the deuterium ions. A fixed pressure is set at the outer magnetic flux surface (
$\unicode[STIX]{x1D70C}=1$
); this is set to a pressure to be expected at the edge of the core, that is, at the top of an edge transport barrier.
The turbulent distribution function, from which can be calculated the turbulent fluxes, is determined by the gyrokinetic equation (A 3). This can be rigorously derived from first principles using assumptions which are valid in a fusion reactor. However, it remains too expensive to solve for the purposes of this study.
Instead, velocity space integrals are taken of this equation to create a hierarchy of fluid moments. In the early 1990s, a set of closures for this hierarchy were developed which accurately captured the linear response for drift waves, as well as finite Larmor radius effects (Hammett & Perkins Reference Hammett and Perkins1990; Dorland & Hammett Reference Dorland and Hammett1993; Beer & Hammett Reference Beer and Hammett1996). Unfortunately, while successful in many cases, these ‘gyrofluid’ models were unable to capture correctly two key properties of the turbulence, namely, the excitation of large-scale zonal flows (Dimits et al.
Reference Dimits, Bateman, Beer, Cohen, Dorland, Hammett, Kim, Kinsey, Kotschenreuther and Kritz2000) and the phenomenon of perpendicular nonlinear phase mixing (Tatsuno et al.
Reference Tatsuno, Dorland, Schekochihin, Plunk, Barnes, Cowley and Howes2009). This resulted in the overprediction of heat fluxes by gyrofluid models. However, in recent years a new hybrid gyrofluid/gyrokinetic code has been developed, Gryfx (see appendices B–D), which overcomes these weaknesses, producing excellent agreement with codes that solve the gyrokinetic equation (figure 7) whilst still taking orders of magnitude less time. A principle advance has been the use of a gyrokinetic solver for the linear zonal flow response (appendix D). This enables Gryfx to capture the characteristic ‘Dimits shift’ (Dimits et al.
Reference Dimits, Bateman, Beer, Cohen, Dorland, Hammett, Kim, Kinsey, Kotschenreuther and Kritz2000), where, for the parameters used in figure 7, the turbulent transport is suppressed to negligible levels in the range
$4\leqslant (1/T_{i})\,\text{d}T_{i}/\text{d}\unicode[STIX]{x1D70C}\lesssim 5$
even though the temperature gradient is above the linear stability threshold. This nonlinear upshift in the threshold has long been attributed to zonal flow dynamics, and our results are consistent with this theory. At large temperature gradients, the turbulence is stronger and the zonal flows are not as effective at suppressing the turbulence. In this regime, new closures that model the effects of nonlinear phase mixing (appendix C) are the dominant effect in reducing the gyrofluid heat flux predictions to gyrokinetic levels. It is also important to note that recent theoretical and numerical work has showed that higher velocity space moments are energetically subdominant in the regimes of interest, producing renewed confidence that a closure with a sufficient number of moments can capture the important dynamics (Parker et al.
Reference Parker, Highcock, Schekochihin and Dellar2016; Schekochihin et al.
Reference Schekochihin, Parker, Highcock, Dellar, Dorland and Hammett2016). In addition, since Gryfx includes the full quadratic nonlinearity (the fourth term on the left-hand side of (A 3)), it is expected to capture important new phenomena such as subcritical turbulence.
Comparison between gyrokinetic simulations (Gs2), the original gyrofluid model (Beer & Hammett Reference Beer and Hammett1996) and Gryfx. The turbulent heat flux is shown as a function of the ion temperature gradient for Cyclone Base Case (Dimits et al. Reference Dimits, Bateman, Beer, Cohen, Dorland, Hammett, Kim, Kinsey, Kotschenreuther and Kritz2000) parameters.
When using Gryfx to calculate the turbulent fluxes within Trinity, it is essential to use sufficient resolution to resolve the turbulent phenomena, to check for any pathologies and to ensure that the turbulent calculation has reached convergence. There are many turbulent phenomena that can make this challenging, particularly near the threshold for the onset of turbulence, including long time scale oscillations and large amplitudes of the zonal flows. This is covered further in appendix E.
Appendix B. Gyrofluid equations
The gyrofluid model in Gryfx is based on the
$4+2$
toroidal gyrofluid model of Beer & Hammett (Reference Beer and Hammett1996), which describes the time evolution of six guiding centre moments of the gyrokinetic equation, (A 3): density (
$n$
), parallel velocity (
$u_{\Vert }$
), parallel and perpendicular temperature (
$T_{\Vert }$
and
$T_{\bot }$
) and the parallel fluxes of parallel and perpendicular heat (
$q_{\Vert }$
and
$q_{\bot }$
). These moments can be defined via the following velocity space integrals of the gyroaveraged fluctuating distribution function
$g_{s}=h_{s}-(Z_{s}e/T_{s})\left\langle \unicode[STIX]{x1D711}\right\rangle _{\mathbf{R}}F_{s}$
:
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D6FF}n=\int g\,\text{d}^{3}v\quad n_{0}\unicode[STIX]{x1D6FF}u_{\Vert }=\int g~v_{\Vert }\,\text{d}^{3}v\\ \displaystyle \unicode[STIX]{x1D6FF}p_{\Vert }=T_{s}\unicode[STIX]{x1D6FF}n+n_{s}\unicode[STIX]{x1D6FF}T_{\Vert }=m_{s}\!\int g~v_{\Vert }^{2}\,\text{d}^{3}v~~\unicode[STIX]{x1D6FF}p_{\bot }=T_{s}\unicode[STIX]{x1D6FF}n+n_{s}\unicode[STIX]{x1D6FF}T_{\bot }=(m_{s}/2)\!\int \!gv_{\bot }^{2}\,\text{d}^{3}v\\ \displaystyle \unicode[STIX]{x1D6FF}q_{\Vert }=-3n_{s}T_{s}\unicode[STIX]{x1D6FF}u_{\Vert }+m_{s}\!\int g~v_{\Vert }^{3}\,\text{d}^{3}v\quad \unicode[STIX]{x1D6FF}q_{\bot }=-n_{s}T_{s}\unicode[STIX]{x1D6FF}u_{\Vert }+(m_{s}/2)\!\int gv_{\Vert }v_{\bot }^{2}\,\text{d}^{3}v.\end{array}\right\} & & \displaystyle\end{eqnarray}$$
These fluctuating quantities (along with the electrostatic potential
$\unicode[STIX]{x1D711}$
) are normalised as
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{a}{\unicode[STIX]{x1D70C}_{i}}\left(\frac{\unicode[STIX]{x1D6FF}n}{n_{s}},\frac{\unicode[STIX]{x1D6FF}u_{\Vert }}{v_{t,s}},\frac{\unicode[STIX]{x1D6FF}T}{T_{s}},\frac{\unicode[STIX]{x1D6FF}q}{n_{s}m_{s}v_{t,s}^{3}},\frac{e\unicode[STIX]{x1D711}}{T_{i}}\right)=\left({\tilde{n}},\tilde{u} _{\Vert },\tilde{T},\tilde{q},\tilde{\unicode[STIX]{x1D711}}\right), & \displaystyle\end{eqnarray}$$
where
$v_{t,s}=\sqrt{T_{s}/m_{s}}$
is the ion thermal speed,
$n_{s}$
and
$T_{s}$
are the equilibrium ion density and temperature respectively,
$\unicode[STIX]{x1D70C}_{s}$
is the ion gyroradius and
$a$
is the normalising length defined to be half the diameter of the last closed flux surface at the midplane. Here the subscript
$i$
denotes the reference ion species. These moment definitions follow Snyder & Hammett (Reference Snyder and Hammett2001) and are consistent with Beer’s original definitions. Note that we will evolve
$\tilde{T}_{\Vert }$
and
$\tilde{T}_{\bot }$
, whereas Beer used
$\tilde{p}_{\Vert }={\tilde{n}}+\tilde{T}_{\Vert }$
and
$\tilde{p}_{\bot }={\tilde{n}}+\tilde{T}_{\bot }$
. Hereafter quantities will be normalised and non-dimensionalised unless otherwise specified, so we will drop the tildes unless they are necessary for clarity.Footnote
6
Further, note that these normalised moments are equivalent to the first several Laguerre–Hermite velocity moments of
$g$
, as described in Mandell et al. (Reference Mandell, Dorland and Landreman2018):
$$\begin{eqnarray}\displaystyle \left(G_{0,0},G_{0,1},\sqrt{2}G_{0,2},\sqrt{6}G_{0,3},G_{1,0},G_{1,1}\right)=\left(n,u_{\Vert },T_{\Vert },q_{\Vert },T_{\bot },q_{\bot }\right), & & \displaystyle\end{eqnarray}$$
where the Laguerre–Hermite moments are defined by (again in non-dimensional form)
$$\begin{eqnarray}G_{\ell ,m}=2\unicode[STIX]{x03C0}\int _{-\infty }^{\infty }\text{d}v_{\Vert }\int _{0}^{\infty }\text{d}\unicode[STIX]{x1D707}\,B~\frac{(-1)^{\ell }}{\sqrt{m!}}\text{L}_{\ell }(\unicode[STIX]{x1D707}B)\text{He}_{m}(v_{\Vert })g\end{eqnarray}$$
with
$\text{L}_{\ell }$
the Laguerre polynomials and
$\text{He}_{m}$
the probabilist’s Hermite polynomials.
The gyrofluid equations that we solve in Gryfx can then be written as
$$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle \frac{\unicode[STIX]{x2202}n}{\unicode[STIX]{x2202}t}+{\mathcal{N}}_{n}+B\unicode[STIX]{x1D735}_{\Vert }\frac{u_{\Vert }}{B}-\left(\frac{1}{L_{ns}}+\frac{1}{2L_{Ts}}\hat{\unicode[STIX]{x1D6FB}}_{\bot }^{2}\right)\text{i}\unicode[STIX]{x1D714}_{\star }\unicode[STIX]{x1D6F7}\nonumber\\ \displaystyle & & \displaystyle \displaystyle \quad +\left(2+\frac{1}{2}\hat{\unicode[STIX]{x1D6FB}}_{\bot }^{2}\right)\text{i}\unicode[STIX]{x1D714}_{d}\unicode[STIX]{x1D6F7}+\text{i}\unicode[STIX]{x1D714}_{d}\left(T_{\Vert }+T_{\bot }+2n\right)=0\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle \frac{\unicode[STIX]{x2202}u_{\Vert }}{\unicode[STIX]{x2202}t}+{\mathcal{N}}_{u_{\Vert }}+B\unicode[STIX]{x1D735}_{\Vert }\frac{n+T_{\Vert }}{B}+\unicode[STIX]{x1D735}_{\Vert }\unicode[STIX]{x1D6F7}+\left(T_{\bot }+n+\frac{1}{2}\hat{\unicode[STIX]{x1D6FB}}_{\bot }^{2}\unicode[STIX]{x1D6F7}\right)\unicode[STIX]{x1D735}_{\Vert }\ln B\nonumber\\ \displaystyle & & \displaystyle \displaystyle \quad +\,\text{i}\unicode[STIX]{x1D714}_{d}\left(q_{\Vert }+q_{\bot }+4u_{\Vert }\right)=0\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle \frac{\unicode[STIX]{x2202}T_{\Vert }}{\unicode[STIX]{x2202}t}+{\mathcal{N}}_{T_{\Vert }}+B\unicode[STIX]{x1D735}_{\Vert }\frac{q_{\Vert }+2u_{\Vert }}{B}+2\left(q_{\bot }+u_{\Vert }\right)\unicode[STIX]{x1D735}_{\Vert }\ln B-\frac{1}{L_{Ts}}\text{i}\unicode[STIX]{x1D714}_{\star }\unicode[STIX]{x1D6F7}+2\text{i}\unicode[STIX]{x1D714}_{d}\unicode[STIX]{x1D6F7}\nonumber\\ \displaystyle & & \displaystyle \displaystyle \quad +\,\text{i}\unicode[STIX]{x1D714}_{d}\left(6T_{\Vert }+2n\right)+2\left|\unicode[STIX]{x1D714}_{d}\right|\left(\unicode[STIX]{x1D708}_{1}T_{\Vert }+\unicode[STIX]{x1D708}_{2}T_{\bot }\right)=-\frac{2}{3}\unicode[STIX]{x1D708}_{ss}\left(T_{\Vert }-T_{\bot }\right)\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle \frac{\unicode[STIX]{x2202}T_{\bot }}{\unicode[STIX]{x2202}t}+{\mathcal{N}}_{T_{\bot }}-B\unicode[STIX]{x1D735}_{\Vert }\frac{u_{\Vert }}{B}+B^{2}\unicode[STIX]{x1D735}_{\Vert }\frac{q_{\bot }+u_{\Vert }}{B^{2}}-\left[\frac{1}{2L_{ns}}\hat{\unicode[STIX]{x1D6FB}}_{\bot }^{2}+\frac{1}{L_{Ts}}\left(1+\hat{\hat{\unicode[STIX]{x1D6FB}}}_{\bot }^{2}\right)\right]\text{i}\unicode[STIX]{x1D714}_{\star }\unicode[STIX]{x1D6F7}\nonumber\\ \displaystyle & & \displaystyle \displaystyle \quad +\left(1+\hat{\unicode[STIX]{x1D6FB}}_{\bot }^{2}+\hat{\hat{\unicode[STIX]{x1D6FB}}}_{\bot }^{2}\right)\text{i}\unicode[STIX]{x1D714}_{d}\unicode[STIX]{x1D6F7}+\text{i}\unicode[STIX]{x1D714}_{d}\left(4T_{\bot }+n\right)+2\left|\unicode[STIX]{x1D714}_{d}\right|\left(\unicode[STIX]{x1D708}_{3}T_{\Vert }+\unicode[STIX]{x1D708}_{4}T_{\bot }\right)\nonumber\\ \displaystyle & & \displaystyle \displaystyle \quad =\frac{1}{3}\unicode[STIX]{x1D708}_{ss}\left(T_{\Vert }-T_{\bot }\right)\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle \frac{\unicode[STIX]{x2202}q_{\Vert }}{\unicode[STIX]{x2202}t}+{\mathcal{N}}_{q_{\Vert }}+\left(3+\unicode[STIX]{x1D6FD}_{\Vert }\right)\unicode[STIX]{x1D735}_{\Vert }T_{\Vert }+\sqrt{2}D_{\Vert }\left|k_{\Vert }\right|q_{\Vert }+\text{i}\unicode[STIX]{x1D714}_{d}\left(-3q_{\Vert }-3q_{\bot }+6u_{\Vert }\right)\nonumber\\ \displaystyle & & \displaystyle \displaystyle \quad +\left|\unicode[STIX]{x1D714}_{d}\right|\left(\unicode[STIX]{x1D708}_{5}u_{\Vert }+\unicode[STIX]{x1D708}_{6}q_{\Vert }+\unicode[STIX]{x1D708}_{7}q_{\bot }\right)=-\unicode[STIX]{x1D708}_{ss}q_{\Vert }\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle \frac{\unicode[STIX]{x2202}q_{\bot }}{\unicode[STIX]{x2202}t}+{\mathcal{N}}_{q_{\bot }}+\unicode[STIX]{x1D735}_{\Vert }\left(T_{\bot }+\frac{1}{2}\hat{\unicode[STIX]{x1D6FB}}_{\bot }^{2}\unicode[STIX]{x1D6F7}\right)+\sqrt{2}D_{\bot }\left|k_{\Vert }\right|q_{\bot }\nonumber\\ \displaystyle & & \displaystyle \quad +\left(T_{\bot }-T_{\Vert }+\hat{\hat{\unicode[STIX]{x1D6FB}}}_{\bot }^{2}\unicode[STIX]{x1D6F7}-\frac{1}{2}\hat{\unicode[STIX]{x1D6FB}}_{\bot }^{2}\unicode[STIX]{x1D6F7}\right)\unicode[STIX]{x1D735}_{\Vert }\ln B+\text{i}\unicode[STIX]{x1D714}_{d}\left(-q_{\Vert }-q_{\bot }+u_{\Vert }\right)\nonumber\\ \displaystyle & & \displaystyle \displaystyle \quad +\left|\unicode[STIX]{x1D714}_{d}\right|\left(\unicode[STIX]{x1D708}_{8}u_{\Vert }+\unicode[STIX]{x1D708}_{9}q_{\Vert }+\unicode[STIX]{x1D708}_{10}q_{\bot }\right)=-\unicode[STIX]{x1D708}_{ss}q_{\bot }\end{eqnarray}$$
where
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D735}_{\Vert }=v_{ts}\hat{\boldsymbol{b}}\boldsymbol{\cdot }\unicode[STIX]{x1D735},\quad b=k_{\bot }^{2}\unicode[STIX]{x1D70C}_{s}^{2},\quad \unicode[STIX]{x1D6F7}=\unicode[STIX]{x1D6E4}_{0}^{1/2}(b)\unicode[STIX]{x1D711},\quad \boldsymbol{v}_{\unicode[STIX]{x1D6F7}}=\hat{\boldsymbol{b}}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D6F7},\\[5.0pt] \displaystyle \frac{1}{2}\hat{\unicode[STIX]{x1D6FB}}_{\bot }^{2}\unicode[STIX]{x1D6F7}=b\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E4}_{0}^{1/2}}{\unicode[STIX]{x2202}b}\unicode[STIX]{x1D711},\quad \hat{\hat{\unicode[STIX]{x1D6FB}}}_{\bot }^{2}\unicode[STIX]{x1D6F7}=b\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}b^{2}}(b\unicode[STIX]{x1D6E4}_{0}^{1/2})\unicode[STIX]{x1D711},\\[7.0pt] \displaystyle \text{i}\unicode[STIX]{x1D714}_{\ast }=-\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\boldsymbol{\cdot }\hat{\boldsymbol{b}}\times \unicode[STIX]{x1D735},\quad \text{i}\unicode[STIX]{x1D714}_{d}=\frac{\unicode[STIX]{x1D70F}_{s}}{Z_{s}B^{2}}\hat{\boldsymbol{b}}\times \unicode[STIX]{x1D735}B\boldsymbol{\cdot }\unicode[STIX]{x1D735},\end{array}\right\} & & \displaystyle\end{eqnarray}$$
and the species thermal velocity
$v_{ts}$
, gyroradius
$\unicode[STIX]{x1D70C}_{s}$
, temperature
$\unicode[STIX]{x1D70F}_{s}$
, charge
$Z_{s}$
, equilibrium density and temperature scale lengths
$L_{ns}$
and
$L_{Ts}$
, and collision frequency
$\unicode[STIX]{x1D708}_{ss}$
have all been non-dimensionalised. The nonlinear terms are denoted by
${\mathcal{N}}_{}$
and will be addressed in detail in appendix C; the final form of these terms is given in (C 10)–(C 15). The quasineutrality constraint for a single ion species is
$$\begin{eqnarray}n_{e}=\frac{n}{1+b/2}-\frac{bT_{\bot }}{2(1+b/2)^{2}}+(\unicode[STIX]{x1D6E4}_{0}-1)\unicode[STIX]{x1D711},\end{eqnarray}$$
where
$\unicode[STIX]{x1D6E4}_{n}(b)=\text{I}_{n}(b)\text{e}^{-b}$
and
$\text{I}_{n}(b)=\text{i}^{-n}\text{J}_{n}(\text{i}b)$
is the modified Bessel function. When electrons are assumed to be adiabatic, which is the case for all results in this paper, we have
$$\begin{eqnarray}n_{e}=\frac{T_{i}}{T_{e}}(\unicode[STIX]{x1D711}-\left\langle \unicode[STIX]{x1D711}\right\rangle _{\unicode[STIX]{x1D713}}),\end{eqnarray}$$
where
$\left\langle \unicode[STIX]{x1D711}\right\rangle _{\unicode[STIX]{x1D713}}$
is a flux surface average.
The coefficients
$D_{\Vert }$
,
$D_{\bot }$
,
$\unicode[STIX]{x1D6FD}_{\Vert }$
and
$\unicode[STIX]{x1D708}_{1}\ldots \unicode[STIX]{x1D708}_{10}$
are set by the closures and taken to be the same as in Beer & Hammett (Reference Beer and Hammett1996). These closure approximations are carefully chosen to capture important kinetic effects, notably Landau damping, phase mixing from toroidal
$\unicode[STIX]{x1D735}B$
and curvature drifts, and finite Larmor radius (FLR) effects. The resulting gyrofluid model can reproduce the gyrokinetic linear dispersion relation quite accurately.
Appendix C. Closures for nonlinear FLR phase mixing
Phase mixing processes, like Landau damping, are fundamentally caused by the fact that particles have a distribution of velocities. For the case of Landau damping, the spread in parallel velocities of particles freely streaming along field lines causes neighbouring particles to move apart. This (linear) phase mixing process smears away spatial perturbations, even in the collisionless limit, and drives the formation of small-scale structures in
$f(v_{\Vert })$
, with
$\unicode[STIX]{x1D6FF}v_{\Vert }\sim 1/k_{\Vert }t$
.
A similar phase mixing process is associated with the nonlinear term in the gyrokinetic equation. This term represents random mixing by gyroaveraged
$\boldsymbol{E}\times \boldsymbol{B}$
flows, and thus produces small scale structure in space. There is a spread in the gyroaveraged
$\boldsymbol{E}\times \boldsymbol{B}$
velocities of particles, as higher energy particles with larger gyroradii average over more fluctuations in the potential and thereby have a slower
$\boldsymbol{E}\times \boldsymbol{B}$
drift than lower energy particles with smaller gyroradii. This leads to phase mixing perpendicular to the field and drives the formation of structure in
$f(v_{\bot })$
. Thus the nonlinear term simultaneously produces small-scale structures in both physical space and perpendicular velocity space.
This process was first recognised in the context of gyrofluid closures by Dorland & Hammett (Reference Dorland and Hammett1993). Later, Schekochihin et al. (Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009) identified the existence of a kinetic cascade due to nonlinear phase mixing. They predicted the key properties of this cascade, and placed it in the context of the broader concept of entropy cascades. Notably, nonlinear phase mixing was found to be the dominant method of generating small-scale structure in velocity space, outpacing Landau damping. Tatsuno et al. (Reference Tatsuno, Plunk, Barnes, Dorland, Howes and Numata2012) studied nonlinear phase mixing in the context of freely decaying turbulence, identifying three regimes of importance, from collisional to collisionless. Howes et al. (Reference Howes, TenBarge, Dorland, Quataert, Schekochihin, Numata and Tatsuno2011) found numerical evidence supporting the existence of nonlinear phase mixing and the entropy cascade at small scales in electromagnetic, kinetic Alfvén wave turbulence. The effects of nonlinear phase mixing have also been observed experimentally, in laboratory magnetised plasmas (Kawamori Reference Kawamori2013) and in the solar wind (Chen et al. Reference Chen, Wicks, Horbury and Schekochihin2010).
We will now derive a gyrofluid closure to model the effects of nonlinear FLR phase mixing. We start with a simple kinetic problem from which we derive a damped kinetic response. We then choose a dissipative gyrofluid closure and fit closure coefficients so that the gyrofluid response closely matches the kinetic response.
C.1 Simple kinetic problem
To illustrate the essence of the nonlinear phase mixing process, we follow Dorland & Hammett (Reference Dorland and Hammett1993) and start with a kinetic picture in sheared slab geometry. We consider a zonal potential that varies sinusoidally in only the
$x$
direction with wavenumber
$k_{x}:\unicode[STIX]{x1D711}=\unicode[STIX]{x1D711}_{ZF}(x)=\unicode[STIX]{x1D719}_{ZF}\sin (k_{x}x)$
. The zonal
$\boldsymbol{E}\times \boldsymbol{B}$
flow is then
$\boldsymbol{v}_{E}=v_{E}(x){\hat{y}}=(\unicode[STIX]{x2202}\unicode[STIX]{x1D711}_{ZF}/\unicode[STIX]{x2202}x){\hat{y}}=k_{x}\unicode[STIX]{x1D719}_{ZF}\cos (k_{x}x){\hat{y}}$
. Assuming no gradients in the equilibrium
$F_{0}$
and no parallel gradients, the gyrokinetic equation reduces to a one-dimensional advection equation involving the nonlinear term,
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}g}{\unicode[STIX]{x2202}t}+J_{0}\left(\frac{k_{x}v_{\bot }}{\unicode[STIX]{x1D6FA}}\right)v_{E}\frac{\unicode[STIX]{x2202}g}{\unicode[STIX]{x2202}y}=0. & & \displaystyle\end{eqnarray}$$
Taking a Maxwellian initial perturbation with a single mode in
$y$
with wavenumber
$k_{y}$
,
$g(t=0)=\text{e}^{\text{i}k_{y}y}F_{M}$
, the perturbation then evolves as
$$\begin{eqnarray}\displaystyle g(t)=F_{M}\text{e}^{\text{i}k_{y}[y-J_{0}(k_{x}v_{\bot }/\unicode[STIX]{x1D6FA})v_{E}t]}\simeq F_{M}\text{e}^{\text{i}k_{y}(y-v_{E}t)}\text{e}^{\text{i}k_{y}bv_{E}tv_{\bot }^{2}/4v_{t}^{2}}, & & \displaystyle\end{eqnarray}$$
where on the right we have expanded to first order in small
$b$
by taking
$J_{0}\simeq 1-(b/4)v_{\bot }^{2}$
, and here
$b=k_{x}^{2}v_{t}^{2}/\unicode[STIX]{x1D6FA}^{2}=k_{x}^{2}\unicode[STIX]{x1D70C}^{2}$
. In this small-
$k_{\bot }\unicode[STIX]{x1D70C}$
limit we can analytically calculate the kinetic perturbed density response, given by
$$\begin{eqnarray}\displaystyle n_{\text{kin}}(t)=\frac{1}{n_{0}}\int \text{d}^{3}v\,g\simeq \text{e}^{\text{i}k_{y}(y-v_{E}t)}\frac{1}{1-\text{i}k_{y}bv_{E}t/2}. & & \displaystyle\end{eqnarray}$$
Thus we see that the density response decays in time with a long tail that goes like
$1/t$
. This is the behaviour that we will want to capture with our fluid closure. Note that we can also numerically integrate the exact kinetic solution in (C 2) to capture the full
$J_{0}$
effects; for
$k_{\bot }\unicode[STIX]{x1D70C}\leqslant 1$
the exact response is nearly identical to the small-
$k_{\bot }\unicode[STIX]{x1D70C}$
response given in (C 3).
C.2 Fluid picture
If we take Laguerre moments in
$\unicode[STIX]{x1D707}B=v_{\bot }^{2}/2$
of our simple one-dimensional gyrokinetic equation (C 1), we will see that each Laguerre moment is coupled via the nonlinear term, which is a manifestation of the phase mixing process (Mandell et al.
Reference Mandell, Dorland and Landreman2018). Examining the equations for the first two Laguerre moments,
$G_{0,0}=n$
and
$G_{1,0}=T_{\bot }$
, and again taking the small-
$k_{\bot }\unicode[STIX]{x1D70C}$
limit, we have
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}n}{\unicode[STIX]{x2202}t}+\left(1-\frac{b}{2}\right)v_{E}\frac{\unicode[STIX]{x2202}n}{\unicode[STIX]{x2202}y}-\frac{b}{2}v_{E}\frac{\unicode[STIX]{x2202}T_{\bot }}{\unicode[STIX]{x2202}y}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}T_{\bot }}{\unicode[STIX]{x2202}t}+\left(1-\frac{3b}{2}\right)v_{E}\frac{\unicode[STIX]{x2202}T_{\bot }}{\unicode[STIX]{x2202}y}-\frac{b}{2}v_{E}\frac{\unicode[STIX]{x2202}n}{\unicode[STIX]{x2202}y}-\frac{b}{2}v_{E}\frac{\unicode[STIX]{x2202}G_{2,0}}{\unicode[STIX]{x2202}y}=0, & \displaystyle\end{eqnarray}$$
where the
$G_{2,0}$
Laguerre moment is not evolved but requires closure. These equations are identical to the Beer equations in the one-dimensional low
$k_{\bot }\unicode[STIX]{x1D70C}$
limit, with the exception of the last
$G_{2,0}$
term in (C 5). Now that we have identified this extra term, we can generalise these equations to three dimensions, and to higher
$k_{\bot }\unicode[STIX]{x1D70C}$
using the full FLR expressions (Dorland & Hammett Reference Dorland and Hammett1993):
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}n}{\unicode[STIX]{x2202}t}+\boldsymbol{v}_{\unicode[STIX]{x1D6F7}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}n+\left[\frac{1}{2}\hat{\unicode[STIX]{x1D6FB}}_{\bot }^{2}\boldsymbol{v}_{\unicode[STIX]{x1D6F7}}\right]\boldsymbol{\cdot }\unicode[STIX]{x1D735}T_{\bot }=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \hspace{-5.0pt}\frac{\unicode[STIX]{x2202}T_{\bot }}{\unicode[STIX]{x2202}t}+\boldsymbol{v}_{\unicode[STIX]{x1D6F7}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}T_{\bot }+\left[\frac{1}{2}\hat{\unicode[STIX]{x1D6FB}}_{\bot }^{2}\boldsymbol{v}_{\unicode[STIX]{x1D6F7}}\right]\boldsymbol{\cdot }\unicode[STIX]{x1D735}n+\left[\hat{\hat{\unicode[STIX]{x1D6FB}}}_{\bot }^{2}\boldsymbol{v}_{\unicode[STIX]{x1D6F7}}\right]\boldsymbol{\cdot }\unicode[STIX]{x1D735}T_{\bot }+\left[\frac{1}{2}\hat{\unicode[STIX]{x1D6FB}}_{\bot }^{2}\boldsymbol{v}_{\unicode[STIX]{x1D6F7}}\right]\boldsymbol{\cdot }\unicode[STIX]{x1D735}G_{2,0}=0. & \displaystyle\end{eqnarray}$$
C.3 Fluid closure
Now we must find a closure expression for the extra term involving
$G_{2,0}$
in (C 7). If we simply set
$G_{2,0}=0$
, as Beer did, we get an oscillatory (undamped) solution for
$n$
and
$T_{\bot }$
. In the spirit of the closures pioneered by Hammett and Perkins to model Landau damping (Hammett & Perkins Reference Hammett and Perkins1990), we can instead use dissipative closures to produce fluid density and temperature responses that are damped like the kinetic responses. Thus we will choose a closure of the form
$$\begin{eqnarray}\displaystyle G_{2,0}=\frac{\left|[\frac{1}{2}\hat{\unicode[STIX]{x1D6FB}}_{\bot }^{2}\boldsymbol{v}_{\unicode[STIX]{x1D6F7}}]\boldsymbol{\cdot }\unicode[STIX]{x1D735}\right|}{[\frac{1}{2}\hat{\unicode[STIX]{x1D6FB}}_{\bot }^{2}\boldsymbol{v}_{\unicode[STIX]{x1D6F7}}]\boldsymbol{\cdot }\unicode[STIX]{x1D735}}\left(\unicode[STIX]{x1D707}_{1}n+\unicode[STIX]{x1D707}_{2}T_{\bot }\right), & & \displaystyle\end{eqnarray}$$
where we follow Beer & Hammett (Reference Beer and Hammett1996) and allow each coefficient
$\unicode[STIX]{x1D707}$
to have a dissipative and a reactive (non-dissipative) piece, given by
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D707}=\unicode[STIX]{x1D707}_{r}+\unicode[STIX]{x1D707}_{i}\frac{\left|[\frac{1}{2}\hat{\unicode[STIX]{x1D6FB}}_{\bot }^{2}\boldsymbol{v}_{\unicode[STIX]{x1D6F7}}]\boldsymbol{\cdot }\unicode[STIX]{x1D735}\right|}{[\frac{1}{2}\hat{\unicode[STIX]{x1D6FB}}_{\bot }^{2}\boldsymbol{v}_{\unicode[STIX]{x1D6F7}}]\boldsymbol{\cdot }\unicode[STIX]{x1D735}}=(\unicode[STIX]{x1D707}_{r},\unicode[STIX]{x1D707}_{i}). & & \displaystyle\end{eqnarray}$$
We set the
$\unicode[STIX]{x1D707}_{1}$
and
$\unicode[STIX]{x1D707}_{2}$
coefficients by numerically minimising the difference between the kinetic density response (C 3) and the fluid density response found by evolving the
$n$
and
$T_{\bot }$
equations after inserting the closures. We have found that
$\unicode[STIX]{x1D707}_{1}=(0.747,-0.078)$
and
$\unicode[STIX]{x1D707}_{2}=(1.368,-2.023)$
produces a fluid response that fits the kinetic response reasonably well for
$k_{\bot }\unicode[STIX]{x1D70C}\lesssim 1$
.
C.4 Extension to other moment equations
One can follow a similar procedure to develop nonlinear phase mixing closures for the remaining gyrofluid equations in the
$4+2$
model. The
$u_{\Vert }$
and
$q_{\bot }$
equations form a coupled system identical to the system we have studied above, so we can use the same 2-moment closure and simply replace
$n$
and
$T_{\bot }$
with
$u_{\Vert }$
and
$q_{\bot }$
, respectively. The
$T_{\Vert }$
and
$q_{\Vert }$
equations are also identical; each of these is not coupled to other equations, so the closures that appear in these equations are 1-moment closures. The closure coefficient is found with a similar procedure to the one used above, by fitting the fluid response for
$T_{\Vert }$
to the corresponding kinetic response.
This completes our derivation of a gyrofluid closure to model nonlinear FLR phase mixing. Adding these new terms to the original nonlinear terms, the final gyrofluid nonlinear terms used in Gryfx are given by
$$\begin{eqnarray}\displaystyle {\mathcal{N}}_{n} & = & \displaystyle \boldsymbol{v}_{\unicode[STIX]{x1D6F7}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}n+\left[{\textstyle \frac{1}{2}}\hat{\unicode[STIX]{x1D6FB}}_{\bot }^{2}\boldsymbol{v}_{\unicode[STIX]{x1D6F7}}\right]\boldsymbol{\cdot }\unicode[STIX]{x1D735}T_{\bot },\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\mathcal{N}}_{u_{\Vert }} & = & \displaystyle \boldsymbol{v}_{\unicode[STIX]{x1D6F7}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}u_{\Vert }+\left[{\textstyle \frac{1}{2}}\hat{\unicode[STIX]{x1D6FB}}_{\bot }^{2}\boldsymbol{v}_{\unicode[STIX]{x1D6F7}}\right]\boldsymbol{\cdot }\unicode[STIX]{x1D735}q_{\bot },\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\mathcal{N}}_{T_{\Vert }} & = & \displaystyle \boldsymbol{v}_{\unicode[STIX]{x1D6F7}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}T_{\Vert }+\left|\left[{\textstyle \frac{1}{2}}\hat{\unicode[STIX]{x1D6FB}}_{\bot }^{2}\boldsymbol{v}_{\unicode[STIX]{x1D6F7}}\right]\boldsymbol{\cdot }\unicode[STIX]{x1D735}\right|\unicode[STIX]{x1D707}_{3}T_{\Vert },\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\mathcal{N}}_{T_{\bot }} & = & \displaystyle \boldsymbol{v}_{\unicode[STIX]{x1D6F7}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}T_{\bot }+\left[{\textstyle \frac{1}{2}}\hat{\unicode[STIX]{x1D6FB}}_{\bot }^{2}\boldsymbol{v}_{\unicode[STIX]{x1D6F7}}\right]\boldsymbol{\cdot }\unicode[STIX]{x1D735}n+\left[\hat{\hat{\unicode[STIX]{x1D6FB}}}_{\bot }^{2}\boldsymbol{v}_{\unicode[STIX]{x1D6F7}}\right]\boldsymbol{\cdot }\unicode[STIX]{x1D735}T_{\bot }\nonumber\\ \displaystyle & & \displaystyle +\left|\left[{\textstyle \frac{1}{2}}\hat{\unicode[STIX]{x1D6FB}}_{\bot }^{2}\boldsymbol{v}_{\unicode[STIX]{x1D6F7}}\right]\boldsymbol{\cdot }\unicode[STIX]{x1D735}\right|\left(\unicode[STIX]{x1D707}_{1}n+\unicode[STIX]{x1D707}_{2}T_{\bot }\right),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\mathcal{N}}_{q_{\Vert }} & = & \displaystyle \boldsymbol{v}_{\unicode[STIX]{x1D6F7}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}q_{\Vert }+\left|\left[{\textstyle \frac{1}{2}}\hat{\unicode[STIX]{x1D6FB}}_{\bot }^{2}\boldsymbol{v}_{\unicode[STIX]{x1D6F7}}\right]\boldsymbol{\cdot }\unicode[STIX]{x1D735}\right|\unicode[STIX]{x1D707}_{3}q_{\Vert },\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\mathcal{N}}_{q_{\bot }} & = & \displaystyle \boldsymbol{v}_{\unicode[STIX]{x1D6F7}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}q_{\bot }+\left[{\textstyle \frac{1}{2}}\hat{\unicode[STIX]{x1D6FB}}_{\bot }^{2}\boldsymbol{v}_{\unicode[STIX]{x1D6F7}}\right]\boldsymbol{\cdot }\unicode[STIX]{x1D735}u_{\Vert }+\left[\hat{\hat{\unicode[STIX]{x1D6FB}}}_{\bot }^{2}\boldsymbol{v}_{\unicode[STIX]{x1D6F7}}\right]\boldsymbol{\cdot }\unicode[STIX]{x1D735}q_{\bot }\nonumber\\ \displaystyle & & \displaystyle +\left|\left[{\textstyle \frac{1}{2}}\hat{\unicode[STIX]{x1D6FB}}_{\bot }^{2}\boldsymbol{v}_{\unicode[STIX]{x1D6F7}}\right]\boldsymbol{\cdot }\unicode[STIX]{x1D735}\right|\left(\unicode[STIX]{x1D707}_{1}u_{\Vert }+\unicode[STIX]{x1D707}_{2}q_{\bot }\right),\end{eqnarray}$$
where
${\mathcal{N}}_{m}$
represents the nonlinear terms in the
$m$
moment equation, and the absolute value terms comprise our new closure terms with closure coefficients
$\unicode[STIX]{x1D707}_{1}=(0.747,-0.078),\unicode[STIX]{x1D707}_{2}=(1.368,-2.023)$
and
$\unicode[STIX]{x1D707}_{3}=(0.456,-0.724)$
. This set of closures is similar to those presented in Dorland & Hammett (Reference Dorland and Hammett1993), but that model required higher order terms of order
$k_{\bot }^{4}\unicode[STIX]{x1D70C}^{4}$
in the
$n$
and
$u_{\Vert }$
equations, which is beyond the order of accuracy of the usual FLR terms. We avoid this by including closure terms proportional to
$n$
and
$u_{\Vert }$
in the
$T_{\bot }$
and
$q_{\bot }$
equations, respectively.
In the derivation of these new closures, we only considered nonlinear FLR phase mixing from a static, one-dimensional potential. In an evolving three-dimensional system, the phase mixing process is much more complicated, and our rough model may overestimate or underestimate the amount of phase mixing. Nonetheless, our model will capture the correct scaling of the mixing, producing physically motivated damping at large
$\unicode[STIX]{x1D711}$
and high
$k_{\bot }$
. In this way our new terms can be interpreted as a type of hyperviscosity, but one that damps at the
$\unicode[STIX]{x1D70C}$
scale as opposed to conventional hyperviscosity models that damp at the grid scale.
Finally we must address how we evaluate and implement a term of the form
${\mathcal{P}}=|\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}|M$
. Since we use a Fourier spectral representation for the equations, we are interested in the Fourier transform of
${\mathcal{P}}$
. Denoting the Fourier transform of a quantity
$x$
as
$\hat{x}_{k}$
, and defining
${\mathcal{N}}_{}=\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}M$
, we evaluate the closure term as
$$\begin{eqnarray}\displaystyle \hat{{\mathcal{P}}}_{k}=\widehat{(|\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}|M)}_{k}=|\hat{{\mathcal{N}}}_{k}|\frac{\hat{M}_{k}}{|\hat{M}_{k}|}, & & \displaystyle\end{eqnarray}$$
where all operations are performed in Fourier space, and
$\hat{{\mathcal{N}}}_{k}$
can be calculated pseudospectrally in the same manner as the usual nonlinear terms.
Appendix D. Hybrid gyrokinetic–gyrofluid zonal flow model
A major drawback of the Beer gyrofluid model is the inability to accurately model zonal flows. These nonlinearly driven sheared poloidal
$\boldsymbol{E}\times \boldsymbol{B}$
flows have been shown to play a key role in determining the turbulence saturation level. Therefore inaccuracy in zonal flow dynamics has been one of the main sources of disagreement between gyrofluid and gyrokinetic turbulence models.Footnote
7
Attempts were made, with limited success, to modify the gyrofluid closures (Beer & Hammett Reference Beer and Hammett1998) to capture the linearly undamped component of the flows derived by Rosenbluth & Hinton (Reference Rosenbluth and Hinton1998). These closure modifications had limited success (Dimits et al.
Reference Dimits, Bateman, Beer, Cohen, Dorland, Hammett, Kim, Kinsey, Kotschenreuther and Kritz2000), and also relied on simplifying assumptions about the magnetic geometry, which would not be able to capture the dependence of zonal flow dynamics on shaping (Xiao, Catto & Dorland Reference Xiao, Catto and Dorland2007). Further, accurately modelling zonal flows involves resolving sharp features in the distribution function in the
$v_{\Vert }$
direction resulting from trapped particle dynamics, which in a moment approach would require much higher Hermite moments in
$v_{\Vert }$
than used in the Beer model (Mandell et al.
Reference Mandell, Dorland and Landreman2018).
Instead of seeking more complicated gyrofluid closure modifications to improve zonal flow accuracy, we avoid this issue by employing a hybrid approach in Gryfx: we evolve the zonal flow modes with a fully gyrokinetic model (using the gyrokinetic code Gs2), while continuing to evolve the non-zonal modes with the gyrofluid model given by (B 5)–(B 10) above. Because we use a Fourier spectral representation for the perpendicular configuration space discretisation in Gryfx, and because the Fourier modes only interact via the nonlinearity, we can easily choose an alternative algorithm for the linear evolution of the zonal (
$k_{y}=0$
) modes. In order to nonlinearly couple the zonal and non-zonal modes, we must first be able to transform the gyrokinetic distribution function into gyrofluid moments, and vice versa. The transformation from the gyrokinetic distribution function to gyrofluid moments is simply the process of taking velocity moments, which can also be interpreted as projecting the distribution function onto a Laguerre–Hermite basis as shown in (B 4). For the inverse transformation, from gyrofluid moments to the gyrokinetic distribution function, we can expand the distribution function in the Laguerre–Hermite basis as (Mandell et al.
Reference Mandell, Dorland and Landreman2018)
$$\begin{eqnarray}\displaystyle \frac{g}{F_{0}} & = & \displaystyle \displaystyle \mathop{\sum }_{\ell =0}\mathop{\sum }_{m=0}\frac{(-1)^{\ell }}{\sqrt{m!}}~\text{L}_{\ell }(\unicode[STIX]{x1D707}B)~\text{He}_{m}(v_{\Vert })~G_{\ell ,m}\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & {\approx} & \displaystyle n+v_{\Vert }u_{\Vert }+\frac{1}{2}\left(v_{\Vert }^{2}-1\right)T_{\Vert }+\left(\unicode[STIX]{x1D707}B-1\right)T_{\bot }\nonumber\\ \displaystyle & & \displaystyle \quad +\,\frac{1}{2}\left(\frac{v_{\Vert }^{3}}{3}-v_{\Vert }\right)q_{\Vert }+v_{\Vert }\left(\unicode[STIX]{x1D707}B-1\right)q_{\bot },\end{eqnarray}$$
where in the second line we have explicitly expressed the expansion in terms of the Beer gyrofluid moments and truncated.
Following the same logic, we can expand the gyrokinetic nonlinear term in terms of the gyrofluid nonlinear terms given in (C 10)–(C 15):
$$\begin{eqnarray}\displaystyle {\mathcal{N}}_{gk} & = & \displaystyle \hat{\boldsymbol{b}}\times \unicode[STIX]{x1D735}\left\langle \unicode[STIX]{x1D711}\right\rangle _{\mathbf{R}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}h=\hat{\boldsymbol{b}}\times \unicode[STIX]{x1D735}\left\langle \unicode[STIX]{x1D711}\right\rangle _{\mathbf{R}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}g\nonumber\\ \displaystyle & {\approx} & \displaystyle {\mathcal{N}}_{n}+v_{\Vert }{\mathcal{N}}_{u_{\Vert }}+\frac{1}{2}\left(v_{\Vert }^{2}-1\right){\mathcal{N}}_{T_{\Vert }}+\left(\unicode[STIX]{x1D707}B-1\right){\mathcal{N}}_{T_{\bot }}\nonumber\\ \displaystyle & & \displaystyle \quad +\,\frac{1}{2}\left(\frac{v_{\Vert }^{3}}{3}-v_{\Vert }\right){\mathcal{N}}_{q_{\Vert }}+v_{\Vert }\left(\unicode[STIX]{x1D707}B-1\right){\mathcal{N}}_{q_{\bot }}.\end{eqnarray}$$
Thus we can construct the gyrokinetic nonlinear term for the
$k_{y}=0$
zonal modes from the
$k_{y}=0$
component of the six gyrofluid nonlinear terms. The full hybrid algorithm then proceeds as follows.
(i) Calculate moments of the
$k_{y}=0$
component of the gyrokinetic distribution function (e.g. via (B 4)).(ii) Evaluate the six gyrofluid nonlinear terms, (C 10)–(C 15), for all Fourier modes.
(iii) Evaluate the
$k_{y}=0$
component of the gyrokinetic nonlinear term via (D 3).(v) Evolve the
$k_{y}=0$
modes with equation (A 3) using the nonlinear term from (iii).
As mentioned above, we couple to the gyrokinetic code Gs2 to evolve the gyrokinetic equation in step (v). Further, we do not include the nonlinear phase mixing closure terms in the
$k_{y}=0$
component of the gyrofluid nonlinear terms. This is in part due to the fact that introducing damping from closure models can suppress the zonal flow residual (which is precisely the reason we moved to gyrokinetic evolution of the zonal modes). Thus in our model, the nonlinear drive for the zonal flows comes only from the lowest several moments of the distribution function; this approximation is justified by the results of Rogers, Dorland & Kotschenreuther (Reference Rogers, Dorland and Kotschenreuther2000).
Finally, note that steps (iv) and (v) can be executed in parallel. In Gryfx we take advantage of GPU–CPU concurrency by simultaneously evolving the non-zonal gyrofluid equations on the GPU and the zonal gyrokinetic equation on the CPU. Using a single GPU and
${\sim}16$
CPU cores (a common supercomputer node configuration), the two steps take roughly the same amount of wall clock time. This means that in our scheme the additional cost (in terms of wall clock time) of using a fully gyrokinetic model for the zonal flows is minimal.
Appendix E. Practical considerations
The presentation of equations and the methodology given above has, for reasons of clarity and brevity, skipped over some of the most thorny issues encountered during the current study. However, we believe that a discussion of them should be present in this work, and we hope that by including it here as an appendix, we can save at least some time and effort for future researchers. In this discussion it helps to think of the operation of CoRFu as a series of six nested loops, as illustrated in figure 8. Within each of these loops issues can arise. The principle practical challenges encountered in the construction and use of the CoRFu framework (independent of the challenges in constructing its individual components such as Trinity or GRYFX) were:
(i) incorporating evolution of the magnetic equilibrium;
(ii) determining whether a particular turbulence calculation had reached steady state;
(iii) encountering failures of either the turbulence code or the Grad–Shafranov (GS) code;
(iv) dealing with turbulent thresholds;
(v) making efficient use of resources; and
(vi) choosing resolution parameters for the turbulence code.
In the following sections we discuss each of these challenges.
Schematic of CoRFu operation showing the hierarchy of decision loops which must be walked through.
E.1 Incorporating evolution of the magnetic equilibrium
The Trinity transport solver has the ability to evolve the magnetic equilibrium internally. However, at present it is not able to treat the evolution of the magnetic equilibrium implicitly, as it does the profiles of density, temperature and flow. Therefore it would be necessary to use a much smaller transport step to avoid numerical instability.
The advantage of this study is that we are only seeking steady-state transport and magnetic equilibrium solutions. Thus, provided Trinity discovers a solution in which both the profiles and the equilibrium do not vary in time, this solution will be self-consistent and thus acceptable. This allows us to separate the evolution of the Grad–Shafranov solution from Trinity, as is illustrated in figure 8. In effect, Trinity evolves using a fixed magnetic equilibrium until the profiles stop changing. At this point the new pressure profile is passed to Chease which recalculates the equilibrium. Trinity is then re-run with new equilibrium and the cycle is repeated until the optimisation criterion converges to within a specified tolerance.
E.2 Determining whether a particular turbulence calculation has reached steady state
This is an extensive subject fraught with complication. Not only is it hard to automate what enters into human judgement in these cases (in many cases the best way of determining whether a simulation has reached saturation is to ‘eyeball’ the time traces) but there are cases when even the best of classification systems would fail (e.g. the case of slow zonal flow amplitude rise in electron-temperature-gradient (ETG) driven turbulence; Colyer et al. Reference Colyer, Schekochihin, Parra, Roach, Barnes, Ghim and Dorland2017).
The difficulty proceeds from having no a priori knowledge of the time scales of a given situation: in particular, those of the linear growth, time to saturation and the longest time scales in the saturated state (usually related to zonal flows). A typical dilemma might be: is the heat flux no longer changing because
(i) the simulation has reached a saturated nonlinear state; or
(ii) the system was stable and all modes died away to a constant noise level; or
(iii) the growth rate was so small that the change in amplitude was undetectable on the time scale allotted?
The authors of this paper have experimented with many techniques, including rolling time averages, moving averages and frequency analysis. In this work, a simpler approach was used which represented a balance between the competing objectives of certainty, efficiency and robustness.
Each turbulence calculation (that is, for each radial location, each element of the Jacobian and each transport iteration) was broken up into a number of stages, with all fields cached between stages so that the calculation could resume instantly from the same state. The number of time steps in each stage was chosen in combination with the time step so that each stage lasted
$O(10)$
turnover times (
$a/v_{\text{thi}}$
). The turbulence code driver then ran the calculation for two initial stages, and then proceeded to run each calculation according to the decision tree presented in figure 9. This decision tree is extremely restrictive about allowing a simulation to be declared converged. This choice was motivated by the observation that a bad value for the heat flux incurred such a great time penalty – by disrupting the sensitive implicit Newton iteration – that running Gryfx for a little longer than was typically necessary saved time in the long run.
Schematic of the CoRFu turbulence code driver showing the logic used to determine if a turbulence calculation is saturated (
$Q_{o}$
is the heat flux from the previous iteration,
$Q_{p}$
is the heat flux from the previous stage of the turbulence calculation and
$Q_{c}$
is the heat flux from the current stage).
E.3 Encountering failures of either Newton iteration, the turbulence code or the GS solver
It is of course possible for any component of CoRFu to experience a failure during the very long operation typically required. This could be due to hardware failures, file corruptions, and so on. However, the three most common failures that were encountered were
(i) failure of the Newton iteration to converge due to overlarge time step or a bad flux value;
(ii) numerical instability from an overly large time step in Gryfx; and
(iii) failure of Chease to converge on a solution, usually as a result of a non-monotonic pressure profile but sometimes from straying to an extreme location in parameter space.
In the case of a failure of the Newton iteration, one likely cause was the inherent noise in the turbulent fluxes. In particular, this was a problem if a large change in a gradient (e.g. the temperature gradient
$\unicode[STIX]{x1D705}_{T}=R/L_{T}$
) produced only a small change in a transport coefficient (e.g. the heat flux
$Q$
). Since we are inverting a Jacobian which essentially depends on quantities such as
$\unicode[STIX]{x0394}Q/\unicode[STIX]{x0394}\unicode[STIX]{x1D705}_{T}$
, if the true value of
$\unicode[STIX]{x0394}Q$
is small, and there is noise, of the same magnitude, which contrives to give
$\unicode[STIX]{x0394}Q$
close to zero, there is an artificial stationary point created in the search domain which will lead the Newton iteration to take a large jump incorrectly. In the future, it is planned to implement a Levenberg–Marquardt mixed search which will be more resistant to this; for the current study it was necessary to be very stringent with the turbulence convergence conditions, and to use a large, fixed value of
$\unicode[STIX]{x0394}\unicode[STIX]{x1D705}_{T}=0.6$
. The reasoning behind this was that since
$Q\propto \unicode[STIX]{x1D705}_{T}^{3}$
(Barnes, Parra & Schekochihin Reference Barnes, Parra and Schekochihin2011b
), a constant value of
$\unicode[STIX]{x0394}\unicode[STIX]{x1D705}_{T}$
would continue to produce a respectable
$\unicode[STIX]{x0394}Q$
as
$Q$
tended to 0, but would not produce an excessively large
$\unicode[STIX]{x0394}Q$
far from the threshold.
With this implemented, a remaining likely cause of a failure was a real stationary point in the search domain, and the solution was to make the transport time step smaller. In the next section we will discuss this in more detail with regard to efficiency. However, there was another consideration when resetting the time step: what to do about the state of the turbulence calculation.
If the Newton iteration had deviated very far from the previous profiles it might be reasonable to assume that since the turbulent state at the end of the failed iteration was very far from that at the beginning, it would be better to start the turbulence calculation from noise to avoid the turbulence displaying violent transients (particularly if the end of the failed iteration was in a strongly driven regime). In fact, since (because of stiff transport) the steady-state profiles were never very far from marginal (with low growth rates), such a strategy would incur a large time penalty and run the risk of mistaking a low-growing turbulent case for a stable case.
In the case of Gryfx, a time step was chosen that was sufficiently small for almost all cases. The turbulence calculation driver would typically start a Gryfx simulation stage using the state of the previous simulation stage (since either it was continuing a calculation till convergence, or it was starting a new calculation with similar parameters to the previous); however this behaviour was changed if necessary (note that resetting the initial condition close to threshold when the growth rate was low incurred a significant time penalty). Otherwise, the following steps were taken to mitigate problems:
(i) If the previous Trinity Newton iteration did not converge, values of fields and moments were kept but time averages of fluxes were reset.
(ii) If the previous fluxes were very small (
${<}10^{-7}$
), all values were reset.(iii) If the previous fluxes were either very large (
${>}1e10$
) or had NaN values, all values were reset. Typically this was because a quirk of the initial condition pushed the turbulence amplitudes to values too large for the time step and could be fixed by resetting, but if this occurred repeatedly the optimisation would be halted and require intervention.(iv) If fewer than 3 turbulence stages had been run the moving averages would be reset to eliminate initial transients.
(v) If greater than 3 turbulence stages had been run, the time scale of the moving average would gradually increase in successive stages (proportional to the total turbulence time elapsed in these successive stages) to force eventual convergence of the heat flux. This was a necessary to ensure that the transport calculation never became stuck, and had no impact on the results, once again because only steady-state solutions were required.
In the case of the Grad–Shafranov solver Chease, since failures were often caused by an unconverged pressure profile, CoRFu first attempted to run Trinity for an extra period of time using the previous equilibrium. Occasionally in preliminary studies (though not in the study reported here) the GS solver failed because the control parameters strayed into a region of parameter space where a solution could not be found. In such cases a highly unfavourable value was returned to prevent the algorithm pursuing that search direction.
E.4 Dealing with turbulent thresholds
To a user of the CoRFu framework, it is an inconvenient reality that optimal solutions typically lie close to the zero-turbulence manifold. This means that the transport solver Trinity will spend most of its operation hovering around the turbulent threshold. This is a problem because of the extra time required to calculate the turbulent transport. If the previous transport iteration was below the threshold and the subsequent iteration is likewise, the turbulence calculation will end rapidly because any initial noise will die away. If the previous iteration was above threshold, and the subsequent one is too, then the turbulence will rapidly adjust to the new drive parameters. However, if the previous iteration was below threshold and the subsequent one above, (i.e. the transport solver crosses the zero-turbulence manifold from below) the turbulence calculation will not converge until amplitudes have grown linearly to saturation amplitudes. Since we are close to threshold with small growth rates, this may take a very long time. Equally unpleasant is when the zero-turbulence manifold is crossed from above: in this case the free energy of the previous vigorous turbulence must be absorbed by only weakly damped modes.
The strategy for dealing with the zero-turbulence manifold was to adjust the way that Trinity calculated the derivatives of fluxes (e.g.
$Q$
) with respect to profile gradients (e.g.
$\unicode[STIX]{x1D705}_{T}$
). Effectively we needed to make sure that, as far as possible, the lower part of the derivative stencil stays below the threshold, and the higher part stays above the threshold. This is achieved firstly by ensuring that the stencil is always fairly broad (at the cost of reducing the convergence rate of the Newton iteration and hence the size of the time step). Fortunately, this requirement coincided with the requirement described in § E.3 for a large
$\unicode[STIX]{x0394}Q$
. The second way to achieve this is to set a tight convergence criterion for the Newton iteration to make sure that it rarely artificially (that is, at odds with the continuous time evolution) crossed the threshold (in terms of Trinity parameters, errflr was set to 0.02 and flrfac was set to 4.0).
E.5 Making efficient use of resources
Although Gryfx is a very fast and efficient code, it is very important, given the number of turbulence calculations carried out, to strictly control the cost of the optimisation study. For many reasons, as given above, it is necessary to be overly cautious and run the turbulence simulations until a strict set of criteria have been satisfied. Nonetheless, it will always be the case that during every iteration of the transport calculation some flux calculations will converge before others. Therefore the flux driver is set up so that the optimisation can run on an arbitrary number of GPUs, and that if the number of GPUs is less than the number of flux calculations carried out at a time, then at each stage of the flux calculation the GPUs will cycle through each of the flux calculations requested by the driver. As the number of stages increases, more and more of the turbulence calculations will have converged and so the GPUs will cycle more quickly through the remaining ones. In the current study, with eight flux calculations per transport iteration, running on four GPUs, this arrangement reduced the simulation cost by close to 50 %.
An equally important consideration with regard to cost is the choices made for iteration and time stepping parameters within Trinity. When seeking a steady-state solution, a certain amount of time will need to elapse while the profiles adjust to the new parameters provided by the optimisation driver. A first assumption might be that increasing the size of the time step will reduce the cost of evolving the profiles for that time (remembering that the Trinity algorithm is implicit and therefore stable for large time steps). However, the larger the time step, the greater the difference between the profiles at the end and the profiles at the beginning of the time step, and thus the greater the chance that the Newton search will not be able to land upon the future profiles. If this occurs, the time step must be reduced and the step repeated, effectively doubling the cost of that time step. Increasing the maximum number of steps in the Newton iteration would increase the chances of it reaching a solution; however, it would also increase the cost should the iteration still fail.
By trial and error, the following choices were found to enable efficient and robust operation.
(i) The threshold for rejecting a Newton iteration was set high (errtol=0.1) to avoid the cost of repeating an iteration (this value was still low enough to avoid wacky iterations triggered by a stationary point).
(ii) The threshold for saying that the Newton iteration was converged and therefore fewer than the allowed number of steps were needed was set low (errflr=0.02) so that the iteration was always as accurate as possible, leading to an accurate time evolution.
(iii) The threshold for increasing the time step was set very low (flrfac=4 so that the threshold was 0.02/4) so that the time step was only increased if the iterations were proceeding very smoothly.
(iv) The factor by which the time step was reduced after a failed iteration was set high (deltadj=4) to avoid repeated failed iterations. For example, when crossing a turbulent threshold it was often necessary to rapidly reduce the transport time step: if deltadj was small, it would take a long time before the time step was reduced sufficiently.
(v) The maximum allowed number of steps in the Newton iteration was set to 4. This represented a balance between giving the search the greatest chance of succeeding and reducing the cost if the search eventually failed.
E.6 Choosing resolution parameters for the turbulence code
The resolution required for the turbulence calculation, including perpendicular box size, perpendicular and parallel resolutions, and time step, vary significantly with the values of the driving parameters and the magnetic geometry.
In the present case, the GPUs that were available for the calculation were sufficiently powerful to allow the increase of the Gryfx resolution to a size that was more than adequate for the majority of simulations (as determined by manual inspection of a sample of turbulence simulations), without an unmanageable increase in cost.
There are a posteriori ways of checking that resolutions were sufficient, such as examining turbulent spectra, which would allow one to avoid such a blanket cost increase. It is also, of course, possible to check that the output quantities such as heat flux do not change with increased resolution. However, the automation of such investigations, which are usually carried out by hand for individual turbulence calculations, has not yet been implemented in the CoRFu framework, and would constitute significant additional research.
