Optimisation of confinement in a fusion reactor using a nonlinear turbulence model
UUnder consideration for publication in J. Plasma Phys. Optimisation of confinement in a fusionreactor using a nonlinear turbulence model
E. G. Highcock † , N. R. Mandell , M. Barnes , and W. Dorland Department of Physics, Chalmers University of Technology, Gothenburg, Sweden Rudolph Peierls Centre for Theoretical Physics, 1 Keble Road, Oxford, UK Culham Centre for Fusion Energy, Culham Science Centre, Abingdon, UK Department of Astrophysical Sciences, Princeton University, Princeton, New Jersey, USA Department of Physics, University of Maryland, College Park, Maryland, USA(Received xx; revised xx; accepted xx)
The confinement of heat in the core of a magnetic fusion reactor is optimised using amultidimensional optimisation algorithm. For the first time in such a study, the lossof heat due to turbulence is modelled at every stage using first-principles nonlinearsimulations which accurately capture the turbulent cascade and large-scale zonal flows.The simulations utilise a novel approach, with gyrofluid treatment of the small-scale driftwaves and gyrokinetic treatment of the large-scale zonal flows. A simple near-circularequilibrium with standard parameters is chosen as the initial condition. The figure ofmerit, fusion power per unit volume, is calculated, and then two control parameters,the elongation and triangularity of the outer flux surface, are varied, with the algorithmseeking to optimise the chosen figure of merit. A two-fold increase in the plasma powerper unit volume is achieved by moving to higher elongation and strongly negativetriangularity.
1. Introduction
Convective heat loss resulting from micro-turbulent fluctuations in a fusion reactorlimits the ability of such a reactor to confine heat to the degree required to achievefusion. Above some critical temperature gradient, this heat loss rises extremely rapidly,limiting the central temperature and hence the fusion reaction rate. Empirically, one wayto improve fusion performance despite this obstacle is to build larger devices. † Currentplanning for the first demonstration power plant (Federici et al. et al. (2014)).The reactor design effort is constantly finding innovative ways to mitigate these † Email address for correspondence: [email protected] † This is based on the observation that the global confinement of heat improves with size (EJDoyle et al. et al. a r X i v : . [ phy s i c s . p l a s m - ph ] A p r E. G. Highcock, N. R. Mandell, M. Barnes and W. Dorland problems. For example, one approach to achieving more compact and cost-effectivedesigns is to use new high-temperature superconductors which have lower cost and greatertensile strength than the superconductors used today, and which can be constructed soas to allow easy disassembly for maintenance (Sorbom et al. et al. et al. et al. et al. (2013)). It is also possible to use more sophisticated “quasilinear” models (Staebler et al. et al. et al. (2015); Jardin et al. (2006))and efforts to predict and optimise performance for particular devices (e.g. Mukhovatov et al. (2003); Budny (2009); Kinsey et al. (2011); Parail et al. (2013); Meneghini et al. (2016)). Additionally, there have been advances in using neural-network approaches toallow such methods to fully replace the simple scaling laws by facilitating extremely fastsimulation based on the quasilinear models (Citrin et al. a ). Nonetheless, nonlinearturbulence calculations, while used routinely for investigation of individual experimentsand verification of quasilinear models, have been conspicuously absent from reactordesign. † 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 studyincluding these models on any reasonable timescale. This absence represents a missedopportunity. 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 thelevels of turbulence can be greatly reduced, without reducing the temperature gradient(Burrell 1997; Synakowski 1999; Highcock et al. et al. a ; Highcock et al. et al. b ). Such effects can be, and have been, included in quasilinearand other models a posteriori . However, even where designs are based on extrapolationfrom the best performing experimental configurations, using the most complete reducedmodels available at the time (e.g. Jardin et al. (2006)), the designs can know nothing ofpotential new phenomena, leading to possible further gains, that might exist in the vastconfiguration space at their disposal (more than two decades ago, Galambos et al. (1995)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 effectnot yet included within the model ‡ (c.f work in the stellarator community which clearly † In stating this we have not overlooked the important work across the entire field (which alack of space prevents us from surveying) studying the performance of and possible improvementsto existing designs such as ITER and DEMO using nonlinear simulations. It is the success ofsuch work that motivates inclusion of nonlinear analysis in the design process itself. ‡ As discussed in more detail below, this proof-of-principle study includes only anelectrostatic, adiabatic-electron, hybrid gyrofluid-gyrokinetic nonlinear model, but the ptimisation of fusion confinement Figure 1.
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. demonstrates the benefit of tight integration of stellarator optimisation with nonlinearanalysis; Xanthopoulos et al. (2014)).Every fusion experiment costs a vast sum to build, and owing to the hard limits placedby 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 studyto search the design space, at several-orders-of-magnitude smaller expense, to identifythe most promising regions. However, in order for this search to be meaningful, we musthave 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 completemodels or experiment and thus preclude the possibility of entirely new design paradigms.By contrast, careful comparison with experiment over recent years has provided a highdegree of confidence that nonlinear models (specifically gyrokinetic models, which usea five-dimensional reduction of the Vlasov equation valid in the conditions of a fusionreactor; for a review see Abel et al. (2013)), without any provision of tuning or fittingparameters, are able to accurately predict what the properties of turbulence will be likein a given situation (White et al. et al. et al.
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 producesexcellent agreement with gyrokinetic models) has been used, as part of a multi-scaletransport 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, theoptimisation algorithm chosen achieved an improvement in the fusion power per unitvolume of 91%. However, this first study is envisioned as a proof of concept, as thestepping stone for a much larger effort in which many more dimensions of parameterspace are explored.A magnetic fusion reactor is composed of a plasma confined by a series of nestedtoroidal magnetic flux surfaces (Fig. 1). The plasma is composed of ionised deuteriumand tritium, which fuse together to produce helium, neutrons, and a large amount ofenergy, provided sufficient pressure is achieved over sufficient volume in the centre of methodology demonstrated is already capable of being extended to include a full gyrokineticelectromagnetic multi-species calculation, albeit at increased cost.
E. G. Highcock, N. R. Mandell, M. Barnes and W. Dorland the plasma. The magnetic field is provided by a series of large external field coils and acurrent which flows through the plasma itself (as well as by various small coils responsiblein addition for the stability of the field). The plasma is heated in a number of differentways: 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 bymeans 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 regioncalled the divertor). In this investigation we consider the core of the reactor; that is weconsider the inner ∼
90% (by minor radius) of the flux surfaces. We assume the existenceof an edge transport barrier, that is, that there is good confinement with a steep pressuregradient in the remaining ∼
10% of the plasma (a standard assumption since the discoveryof the “high-confinement mode”; Wagner et al. (1984)). We also assume that the barrieris independent of elongation and triangularity (this assumption, necessary because weare considering only core transport, differs from current experimental observations, andis discussed further below).The software used in this study to optimise such a reactor is the
CodeRunner
FusionOptimisation Framework (
CoRFu ). This is a sizeable edifice of software comprising manyindependent parts, as illustrated in Fig. 2. The equations underlying the framework aregiven in Appendices A–D. The calculation proceeds as follows.To begin with, the choice of optimisation control parameters, and all other requisiteinput information, is provided to the
CoRFu framework, which acts as the optimisationdriver. From this input, an initial configuration is assembled, comprising parameterisa-tions of: • the shape of the flux surface at the edge of the plasma core; • the profile of the toroidal current, and • the initial pressure profile, including a fixed finite pressure at the edge of the plasmacore.These are used (Link 1 in Fig. 2) by the Chease code (L¨utjens et al. purely hypothetical tokamak. We do, however,choose machine parameters (see Table 1) that are comparable in magnitude to those ofa large tokamak—the Joint European Torus (JET)—and profiles that are plausible forsuch a tokamak (Fig. 5; see e.g. Barnes et al. (2011 a ) or the tokamak profile databasedescribed in Roach et al. (2008) for a comparison).The initial configuration, including the magnetic equilibrium solution, details of exter-nal heat sources and the profiles of particle density (densities and external sources areheld fixed) are then used (Links 2 and 3) by the Trinity transport solver (Barnes et al. † This is done by allowing the pressure profileto 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 theheat inputs, both external and that generated by fusion. However, since the magneticequilibrium itself is affected by the pressure, the result of the
Trinity calculation mustbe fed back to
Chease (Link 4), a new equilibrium generated, and
Trinity rerun, untilthe cycle converges, the pressure stops changing and a steady-state solution is found. † This limits the possibility of interplay between the electron and ion heating and losschannels, either to good or ill; we note again that this and other limitations of thisproof-of-principles study are discussed more fully below. ptimisation of fusion confinement Figure 2.
Schematic of the
CoRFu framework.
To calculate the heat loss due to turbulence (Links 5 and 8),
Trinity uses severalsimultaneous copies of the GPU-based hybrid gyrofluid/gyrokinetic code
Gryfx (seeAppendices B–D). Each copy of
Gryfx calculates the heat flux resulting from theturbulence at a particular location in the plasma, that is, on a particular magnetic fluxsurface, given the pressure profile and magnetic equilibrium. A mean-field, multi-scaleapproach 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 turbulentheat flux is recalculated. An illustration of the complete system is given in Fig. 3.The
Gryfx code divides the calculation of the heat flux in two (Links 6 and 7) usinga novel algorithm: the evolution of the smaller scale drift waves is calculated using thegyrofluid model (Dorland & Hammett 1993; Beer & Hammett 1996) by
Gryfx , and theevolution of the large-scale zonal flows which are generated by the turbulence is calculatedby the gyrokinetic code
Gs2 (Dorland et al. et al.
E. G. Highcock, N. R. Mandell, M. Barnes and W. Dorland
Figure 3.
An illustration of a
Trinity simulation showing above-threshold turbulence for theoptimal configuration. The elongated ribbons to the left of the main picture are flux tubes (Beer et al. being calculated in
Gryfx (see Appendix D). The simulations are electrostatic with aBoltzmann response for electrons.The transport calculation within
Trinity used only five radial grid points. This wasto 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 thefact that the transport calculations evolved to steady state (where power balance wassatisfied on each flux surface), will reduce the impact of the sparse radial grid on theresults of this study.
3. Finding the optimal solution
When the pressure, magnetic equilibrium, and turbulent heat fluxes have ceased toevolve, a steady-state solution has been found. At this point, the figure of merit forthis solution (the fusion power per unit volume) is used by the optimisation driver
CodeRunner (Link 9 in Fig. 2), to generate a new set of control parameters accordingto the particular optimisation algorithm being used. In this study, the Nelder-Meadsimplex algorithm was selected † (not to be confused with the simplex algorithm in linear † The simplex algorithm was chosen for its reliability in this initial study. However, the
CoRFu framework is expected to be used to explore a multi-dimensional design spacewith a plethora of local maxima. This next step will first require a careful study of thesuitability of various multi-dimensional optimisation algorithms both standard and novel,whether heuristics like simulated annealing or genetic algorithms, and whether gradient-freelike Nelder-Mead, or pseudo-gradient like BFGS and variants. Lewis (2004) provides a good ptimisation of fusion confinement et al. (1996)). The shape of the magnetic field haslong been known to have a marked effect on the turbulence, and the choice of these twoshaping parameters in particular was motivated by pioneering experimental work on theTCV tokamak, which was designed to allow large variation in both—demonstrating aconsequent large variation in performance (Weisen et al. κ , δ ) =(1.3, 0.2), (1.3, 0.1) and (1.4, 0.1). The lowest initial vertex is (1.3,0.2) with a powerper unit volume ( P/V ) of 0.0313 MW m − . Initially the algorithm moves to higherelongation and lower triangularity, before proceeding to keep elongation roughly constantand to continue to decrease triangularity. The iteration was terminated when values oftriangularity moved significantly beyond the limits of what has been seen in experiment( ∼ -0.65 Pochelon et al. (1999)). The final, optimal value of P/V of 0.0598 MW m − with( κ , δ ) = (1.6, -0.625). Thus, an improvement of 91% was discovered over the course ofthe optimisation. Since the final value of triangularity was somewhat extreme, we alsoprovide data for an intermediate result with ( κ , δ ) = (1.525, -0.25), for which P/V is0.0481 MW m − .It is interesting also to consider the evolution of the confinement time, defined as thetotal amount of stored kinetic energy divided by the rate of energy injection, is displayedin Fig. 6. We also see that optimising the confinement time, which reflects the ability ofthe system to confine heat, would have produced different results to optimising the fusionpower per unit volume; in particular, the improvement at negative triangularity is lessmarked, and moving to higher elongation at positive triangularity may have produced asimilar improvement in the confinement time.
4. Discussion and Outlook
It has been shown that, without any tuning of model parameters, it has been possibleto generate a plausible solution for reactor performance and then optimise it. Thefinal fusion power in this study is small compared to the input heating power (onthe ions) of 20MW, but that is because the parameters in this study (Tab. 1) arenot reactor-like, but rather of the same order of magnitude as those of JET, thelargest existing device: a choice made to allow assessment of the credibility of theinitial solution. Even more encouragingly, the qualitative result of this optimisation—that core confinement improved at lower triangularity and higher elongation—is directly overview. The field of nonlinear multi-dimensional optimisation is large and increasing, withboth commercial (ESTECO 2018) and freely available Abramson et al. (2006) software offerings.Wide, too, is the literature on the comparison of algorithms, whether individual examples(Manousopoulos & Michalopoulos 2009) or more abstract papers investigating the techniques ofalgorithm comparison, or metamodelling (Simpson et al.
E. G. Highcock, N. R. Mandell, M. Barnes and W. Dorland . . . . − . − . − . . . . . . . . . . P / V κδ P / V . . . . . . . . .
35 1 . .
45 1 . .
55 1 . . κ − . − . − . − . − . − . . . δ . . . . . . . P / V (b) Figure 4.
Top . Diagram illustrating the search path taken by the simplex optimisationalgorithm. The suspended triangles represent the simplex at each iteration, with the colourand 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 actualsequence of ( κ , δ ) values evaluated by Trinity . Bottom . P/V as a function of κ and δ , boththe function values (squares) and an interpolated surface, with the colour being the value of P/V . The highest (optimal) value is at κ = 1 . δ = − . κ = 1 . δ = − . corroborated by experimental measurements on the TCV tokamak (Weisen et al. et al. ptimisation of fusion confinement . . . . T i ( k e V ) InitialIntermediateOptimal n e ( m − ) q S i ( M W / m − ) ρ Initial Intermediate Optimal
Figure 5.
Top . Profiles of ion temperature ( T i ), electron density ( n e ), safety factor ( q ) and heatinput ( S i ), showing the values of each quantity versus normalised radius ρ . ρ is a dimensionlessquantity which labels the flux surface. In this work we examine core transport; thus at theouter flux surface ( ρ = 1) there is a finite temperature of 3.2 keV and a finite electron densityof 0.6 × m − . Bottom . Outer flux surface shapes of the initial, intermediate and optimalsolutions. all a total of 8680 converged nonlinear turbulence calculations at a spatial resolution of192 × ×
20 (radial × poloidal × parallel) were required, at a total cost of approximately3000 GPU-accelerated node hours. That such an exercise was possible was owing to twodevelopments. The first was the design of the implicit algorithm in the Trinity codewhich reduced the number of turbulence calculations required (Barnes et al.
E. G. Highcock, N. R. Mandell, M. Barnes and W. Dorland . .
35 1 . .
45 1 . .
55 1 . . κ − . − . − . − . − . − . . . δ . . . . . . . τ E / s Figure 6.
Evolution of the confinement time in seconds, given as a function of elongation κ and triangularity δ . Squares represent data points, shown on top of an interpolated surface forillustration of the overall trend. Initial Intermediate OptimalOn-axis Magnetic Field. (T) 2.91 2.84 2.81On-axis Major Radius (m) 3.12 3.17 3.21Ion input heat S i (MW) 20 20 20Plasma current I p (MA) 1.9 1.9 1.9Elongation κ δ V (m ) 59.3 67.1 71.8Confinement time (s) 0.179 0.221 0.237Fusion Power P (MW) 1.85 3.23 4.30 P/V (MW/m ) 0.0313 0.0481 0.0598 Table 1.
Global properties of the initial, intermediate and optimal configurations second was the development of
Gryfx , which by using a gyrokinetic response for thezonal flows and advanced nonlinear closures was able to overcome the shortcomings ofgyrofluid codes in the 1990s (Dimits et al.
Gs2 run; Mandell & Dorland (2014)).It is only right to point out in this discussion that this study focuses solely on theplasma core, and not on the reactor components or the plasma edge. It should benoted that though edge transport barriers can be achieved in combination with negativetriangularity (Pochelon et al. et al. et al. (2017)). Thus further work is needed tounderstand how these results would be modified if the plasma edge was modelled self- ptimisation of fusion confinement
Gryfx to consider electromagnetic effects and kinetic electron physics, and toinclude the effects of rotation, which are not considered here. Other priorities would beto include additional transport channels: electron and ion heat and particles, momentum,impurity transport and so on. The complex interplay between these channels can be thecrucial factor in achieving higher performance (see e.g. Highcock et al. (2011)), and robustoperation, and is an ideal area to study with a nonlinear framework such as
CoRFu . Theeventual goal is to switch to using optimisation algorithms which parallelize easily andcan search for global maxima, and run a vast parallel blue-skies search for a dramaticallyoptimised fusion reactor.
5. Acknowledgements
Appendix A. Equations underlying the optimisation framework
The equations that govern the system are presented in (Abel et al. (2013)), which buildsupon an already large body of work (significantly Frieman & Chen (1982); Sugama &Horton (1998)).They follow a multi-scale approach, with a clear separation in time between theevolution of the safety factor (the twist of the magnetic field lines, which evolves onthe resistive timescale; Abel & Cowley (2013)), the evolution of the pressure (whichevolves on the confinement timescale), and the evolution of the turbulence (which occurson the timescale of the linear micro-instabilities which drive it). The equations governing2
E. G. Highcock, N. R. Mandell, M. Barnes and W. Dorland all three are ∂q∂t = c π ∂∂ψ V (cid:48) (cid:104) E · B (cid:105) ψ (A 1)32 1 V (cid:48) ∂∂t V (cid:48) (cid:104) n s (cid:105) ψ T s + 1 V (cid:48) ∂∂ψ V (cid:48) (cid:104) Q s (cid:105) ψ = S s (A 2)and ∂h s ∂t + (cid:16) v (cid:107) ˆ b + V D s + cB ˆ b × ∇ (cid:104) ϕ (cid:105) R (cid:17) · ∇ h s = (cid:104) C [ h s ] (cid:105) R + Z s eF s T s ∂ (cid:104) ϕ (cid:105) R ∂t − ∂F s ∂ψ (cid:16) ˆ b × ∇ (cid:104) ϕ (cid:105) R (cid:17) · ∇ ψ, (A 3) respectively. In these equations, q is the magnetic safety factor, ψ is the poloidal magneticflux which is contained within a flux surface, V is the volume of the flux surface, V (cid:48) = dV /dψ is the incremental volume element (loosely, the flux surface area), E is the electricfield and B the magnetic field, and (cid:104)(cid:105) ψ 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) partof the particle distribution function, v (cid:107) is the velocity along the field line, ˆ b = B / | B | , V D s represents the effect of the magnetic field inhomogeneities on particle motion, ϕ isthe 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 thespecies distribution function. The operator (cid:104)(cid:105) R denotes an average over gyrophase atfixed gyrocentre location R . The turbulent component of the heat flux Q s can be triviallyobtained from the turbulent distribution function h s Q s = (cid:28)(cid:90) d v m s v (cid:104)(cid:16) cB ˆ b × ∇ (cid:104) ϕ (cid:105) R (cid:17) · ∇ ψ (cid:105) (cid:104) h s (cid:105) r (cid:29) ψ,t (A 4)where (cid:104)(cid:105) r denotes a gyroaverage at constant particle position. The (usually) sub-dominant neoclassical component of the heat flux is modelled using analytical approxi-mations (Chang & Hinton 1982), as are radiative losses due to Bremsstrahlung (Glasstone& Lovberg 1960).Since we are seeking a steady state solution, equation (A 1) is not evolved duringthis study. Instead we prescribe a fixed profile of surface averaged toroidal current. Thiscauses q to change with the pressure gradient, not on the resistive timescale, during theevolution towards steady state, but once steady state is reached q ceases to evolve and thesolution is self-consistent. Using the prescribed profile of toroidal current, a prescribedouter flux surface and the pressure profile, the poloidal flux (and hence the shape of themagnetic flux surfaces) can be determined using the Grad-Shafranov equation: ∆ ∗ ψ = − πR (cid:88) s n s T s (cid:26) d ln n s dψ + d ln T s dψ (cid:27) − I ( ψ ) dIdψ , (A 5)where ∆ ∗ is the Grad-Shafranov operator ∆ ∗ ψ = (cid:18) ∂ ∂R − R ∂∂R + ∂ ∂Z (cid:19) ψ. (A 6)Note that when using the Chease code (L¨utjens et al. (1996)) to solve the Grad-Shafranov equation, the usual function
I∂I/dψ can be replaced as an input by the surfaceaveraged toroidal current density, which is what is specified in this study. ptimisation of fusion confinement
Trinity (Barnes et al. (2010);
Trinity is also capable of evolving the density and rotation, which are kept fixed forthis study). This equation determines how the evolution of the pressure is governed bythe 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-Shafranovequation itself varies on the same timescale as the pressure, the solution is periodicallyupdated with the new pressure profile until steady-state is reached. The pressure equationcan be evolved separately for each charged species; in this study, to reduce expense,it is solved for the deuterium ions. It is then assumed that collisional processes willrapidly equilibrate the temperatures of all species, and hence that the temperatures ofthe electrons and tritium ions is the same as that of the deuterium ions. A fixed pressureis set at the outer magnetic flux surface ( ρ = 1); this is set to a pressure to be expectedat 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 firstprinciples using assumptions which are valid in a fusion reactor. However, it remains tooexpensive to solve for the purposes of this study.Instead, velocity space integrals are taken of this equation to create a hierarchy of fluidmoments. In the early 1990s, a set of closures for this hierarchy were developed whichaccurately captured the linear response for drift waves, as well as finite Larmor radiuseffects (Beer & Hammett 1996; Dorland & Hammett 1993; Hammett & Perkins 1990).Unfortunately, while successful in many cases, these “gyrofluid” models were unable tocapture correctly two key properties of the turbulence, namely, the excitation of large-scale zonal flows (Dimits et al. et al. Gryfx (see Appendices B–D), which overcomes these weaknesses,producing excellent agreement with codes that solve the gyrokinetic equation (Fig. 7)whilst still taking orders of magnitude less time. As is described in the main text, aprinciple 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. (cid:54) (1 /T i ) dT i /dρ (cid:46) et al. et al. Gryfx includesthe full quadratic nonlinearity (the fourth term on the left hand side of equation (A 3)),it is expected to capture important new phenomena such as subcritical turbulence.When using
Gryfx to calculate the turbulent fluxes within
Trinity , it is essential touse sufficient resolution to resolve the turbulent phenomena, to check for any pathologies,and to ensure that the turbulent calculation has reached convergence. There are manyturbulent phenomena that can make this challenging, particularly near the threshold for4
E. G. Highcock, N. R. Mandell, M. Barnes and W. Dorland Q i ( g y r o B o h m U n i t s ) Logarithmic Temperature Gradient (1 /T i ) dT i /dρ Gs2
Original Gyrofluid
Gryfx
Figure 7.
Comparison between gyrokinetic simulations (
Gs2 ), the original gyrofluid model(Beer & Hammett (1996)), and
Gryfx . The turbulent heat flux is shown as a function of theion temperature gradient for Cyclone Base Case (Dimits et al. the onset of turbulence, including long timescale oscillations and large amplitudes of thezonal 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 1996), which describes the time evolution of six guiding center moments of thegyrokinetic equation, (A 3): density ( n ), parallel velocity ( u (cid:107) ), parallel and perpendiculartemperature ( T (cid:107) and T ⊥ ), and the parallel fluxes of parallel and perpendicular heat ( q (cid:107) and q ⊥ ). These moments can be defined via the following velocity space integrals of thegyroaveraged fluctuating distribution function g s = h s − ( Z s e/T s ) (cid:104) ϕ (cid:105) R F s : δn = ∫ g d vδp (cid:107) = T s δn + n s δT (cid:107) = m s ∫ g v (cid:107) d vδq (cid:107) = − n s T s δu (cid:107) + m s ∫ g v (cid:107) d v n δu (cid:107) = ∫ g v (cid:107) d vδp ⊥ = T s δn + n s δT ⊥ = ( m s / ∫ g v ⊥ d vδq ⊥ = − n s T s δu (cid:107) + ( m s / ∫ g v (cid:107) v ⊥ d v These fluctuating quantities (along with the electrostatic potential ϕ ) are normalized as aρ i (cid:18) δnn s , δu (cid:107) v t,s , δTT s , δqn s m s v t,s , eϕT i (cid:19) = (cid:16) ˜ n, ˜ u (cid:107) , ˜ T , ˜ q, ˜ ϕ (cid:17) (B 1)where v t,s = (cid:112) T s /m s is the ion thermal speed, n s and T s are the equilibrium ion densityand temperature respectively, ρ s is the ion gyroradius, and a is the normalizing lengthdefined to be half the diameter of the last closed flux surface at the midplane. Here thesubscript i denotes the reference ion species. These moment definitions follow (Snyder& Hammett 2001) and are consistent with Beer’s original definitions. Note that we willevolve ˜ T (cid:107) and ˜ T ⊥ , whereas Beer used ˜ p (cid:107) = ˜ n + ˜ T (cid:107) and ˜ p ⊥ = ˜ n + ˜ T ⊥ . Hereafter quantitieswill be normalised and non-dimensionalised unless otherwise specified, so we will drop thetildes unless they are necessary for clarity. † Further, note that these normalised momentsare equivalent to the first several Laguerre-Hermite velocity moments of g , as described † Our non-dimensionalisation is chosen to be consistent with Appendix A of (Mandell et al. ptimisation of fusion confinement et al. (cid:16) G , , G , , √ G , , √ G , , G , , G , (cid:17) = (cid:0) n, u (cid:107) , T (cid:107) , q (cid:107) , T ⊥ , q ⊥ (cid:1) , (B 2)where the Laguerre-Hermite moments are defined by (again in non-dimensional form) G (cid:96),m = 2 π (cid:90) ∞−∞ dv (cid:107) (cid:90) ∞ dµB ( − (cid:96) √ m ! L (cid:96) ( µB )He m ( v (cid:107) ) g (B 3)with L (cid:96) the Laguerre polynomials and He m the probabilist’s Hermite polynomials.The gyrofluid equations that we solve in Gryfx can then be written as ∂n∂t + N n + B ∇ (cid:107) u (cid:107) B − (cid:18) L ns + 12 L T s ˆ ∇ ⊥ (cid:19) iω (cid:63) Φ + (cid:18) ∇ ⊥ (cid:19) iω d Φ + iω d (cid:0) T (cid:107) + T ⊥ + 2 n (cid:1) = 0 (B 4) ∂u (cid:107) ∂t + N u (cid:107) + B ∇ (cid:107) n + T (cid:107) B + ∇ (cid:107) Φ + (cid:18) T ⊥ + n + 12 ˆ ∇ ⊥ Φ (cid:19) ∇ (cid:107) ln B + iω d (cid:0) q (cid:107) + q ⊥ + 4 u (cid:107) (cid:1) = 0 (B 5) ∂T (cid:107) ∂t + N T (cid:107) + B ∇ (cid:107) q (cid:107) + 2 u (cid:107) B + 2 (cid:0) q ⊥ + u (cid:107) (cid:1) ∇ (cid:107) ln B − L T s iω (cid:63) Φ + 2 iω d Φ + iω d (cid:0) T (cid:107) + 2 n (cid:1) + 2 | ω d | (cid:0) ν T (cid:107) + ν T ⊥ (cid:1) = − ν ss (cid:0) T (cid:107) − T ⊥ (cid:1) (B 6) ∂T ⊥ ∂t + N T ⊥ − B ∇ (cid:107) u (cid:107) B + B ∇ (cid:107) q ⊥ + u (cid:107) B − (cid:20) L ns ˆ ∇ ⊥ + 1 L T s (cid:16) ∇ ⊥ (cid:17)(cid:21) iω (cid:63) Φ + (cid:16) ∇ ⊥ + ˆˆ ∇ ⊥ (cid:17) iω d Φ + iω d (4 T ⊥ + n ) + 2 | ω d | (cid:0) ν T (cid:107) + ν T ⊥ (cid:1) = 13 ν ss (cid:0) T (cid:107) − T ⊥ (cid:1) (B 7) ∂q (cid:107) ∂t + N q (cid:107) + (cid:0) β (cid:107) (cid:1) ∇ (cid:107) T (cid:107) + √ D (cid:107) (cid:12)(cid:12) k (cid:107) (cid:12)(cid:12) q (cid:107) + iω d (cid:0) − q (cid:107) − q ⊥ + 6 u (cid:107) (cid:1) + | ω d | (cid:0) ν u (cid:107) + ν q (cid:107) + ν q ⊥ (cid:1) = − ν ss q (cid:107) (B 8) ∂q ⊥ ∂t + N q ⊥ + ∇ (cid:107) (cid:18) T ⊥ + 12 ˆ ∇ ⊥ Φ (cid:19) + √ D ⊥ (cid:12)(cid:12) k (cid:107) (cid:12)(cid:12) q ⊥ + (cid:18) T ⊥ − T (cid:107) + ˆˆ ∇ ⊥ Φ −
12 ˆ ∇ ⊥ Φ (cid:19) ∇ (cid:107) ln B + iω d (cid:0) − q (cid:107) − q ⊥ + u (cid:107) (cid:1) + | ω d | (cid:0) ν u (cid:107) + ν q (cid:107) + ν q ⊥ (cid:1) = − ν ss q ⊥ (B 9)where ∇ (cid:107) = v ts ˆ b · ∇ , b = k ⊥ ρ s , Φ = Γ / ( b ) ϕ, v Φ = ˆ b × ∇ Φ,
12 ˆ ∇ ⊥ Φ = b ∂Γ / ∂b ϕ, ˆˆ ∇ ⊥ Φ = b ∂ ∂b ( bΓ / ) ϕ,iω ∗ = −∇ ψ · ˆ b × ∇ , iω d = τ s Z s B ˆ b × ∇ B · ∇ , and the species thermal velocity v ts , gyroradius ρ s , temperature τ s , charge Z s , equilibriumdensity and temperature scale lengths L ns and L T s , and collision frequency ν ss have allbeen non-dimensionalised. The nonlinear terms are denoted by N and will be addressedin detail in Appendix C; the final form of these terms is given in (C 10–C 15). Thequasineutrality constraint for a single ion species is n e = n b/ − bT ⊥ b/ + ( Γ − ϕ, (B 10)6 E. G. Highcock, N. R. Mandell, M. Barnes and W. Dorland where Γ n ( b ) = I n ( b ) e − b and I n ( b ) = i − n J n ( ib ) is the modified Bessel function. Whenelectrons are assumed to be adiabatic, which is the case for all results in this paper, wehave n e = T i T e ( ϕ − (cid:104) ϕ (cid:105) ψ ) , (B 11)where (cid:104) ϕ (cid:105) ψ is a flux surface average.The coefficients D (cid:107) , D ⊥ , β (cid:107) , and ν − ν are set by the closures and taken to be thesame as in (Beer & Hammett 1996). These closure approximations are carefully chosen tocapture important kinetic effects, notably Landau damping, phase mixing from toroidal ∇ B and curvature drifts, and finite Larmor radius (FLR) effects. The resulting gyrofluidmodel 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 factthat particles have a distribution of velocities. For the case of Landau damping, thespread in parallel velocities of particles freely streaming along field lines causes neigh-boring particles to move apart. This (linear) phase mixing process smears away spatialperturbations, even in the collisionless limit, and drives the formation of small-scalesstructures in f ( v (cid:107) ), with δv (cid:107) ∼ /k (cid:107) t .A similar phase mixing process is associated with the nonlinear term in the gyrokineticequation. This term represents random mixing by gyro-averaged E × B flows, and thusproduces small scale structure in space. There is a spread in the gyro-averaged E × B velocities of particles, as higher energy particles with larger gyroradii average over morefluctuations in the potential and thereby have a slower E × B drift than lower energyparticles with smaller gyroradii. This leads to phase mixing perpendicular to the fieldand drives the formation of structure in f ( v ⊥ ). Thus the nonlinear term simultaneouslyproduces small scale structures in both physical space and perpendicular velocity space.This process was first recognized in the context of gyrofluid closures by Dorland &Hammett (1993). Later, Schekochihin et al. (2009) identified the existence of a kineticcascade 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, nonlinearphase mixing was found to be the dominant method of generating small-scale structurein velocity space, outpacing Landau damping. Tatsuno et al. (2012) studied nonlinearphase mixing in the context of freely decaying turbulence, identifying three regimesof importance, from collisional to collisionless. Howes et al. (2011) found numericalevidence supporting the existence of nonlinear phase mixing and the entropy cascade atsmall scales in electromagnetic, kinetic Alfv´en wave turbulence. The effects of nonlinearphase mixing have also been observed experimentally, in laboratory magnetized plasmas(Kawamori 2013) and in the solar wind (Chen et al. Simple kinetic problem
To illustrate the essence of the nonlinear phase mixing process, we follow Dorland &Hammett (1993) and start with a kinetic picture in sheared slab geometry. We considera zonal potential that varies sinusoidally in only the x direction with wavenumber k x : ptimisation of fusion confinement ϕ = ϕ ZF ( x ) = φ ZF sin( k x x ). The zonal E × B flow is then v E = v E ( x )ˆ y = ∂ϕ ZF ∂x ˆ y = k x φ ZF cos( k x x )ˆ y . Assuming no gradients in the equilibrium F and no parallel gradients,the gyrokinetic equation reduces to a one-dimensional advection equation involving thenonlinear term, ∂g∂t + J (cid:18) k x v ⊥ Ω (cid:19) v E ∂g∂y = 0 . (C 1)Taking a Maxwellian initial perturbation with a single mode in y with wavenumber k y , g ( t = 0) = e ik y y F M , the perturbation then evolves as g ( t ) = F M e ik y [ y − J ( k x v ⊥ /Ω ) v E t ] (cid:39) F M e ik y ( y − v E t ) e ik y bv E tv ⊥ / v t , (C 2)where on the right we have expanded to first order in small b by taking J (cid:39) − b v ⊥ ,and here b = k x v t /Ω = k x ρ . In this small- k ⊥ ρ limit we can analytically calculate thekinetic perturbed density response, given by n kin ( t ) = 1 n (cid:90) d v g (cid:39) e ik y ( y − v E t ) − ik y bv E t/ . (C 3)Thus we see that the density response decays in time with a long tail that goes like 1 /t .This is the behavior that we will want to capture with our fluid closure. Note that wecan also numerically integrate the exact kinetic solution in (C 2) to capture the full J effects; for k ⊥ ρ (cid:54) k ⊥ ρ response givenin (C 3). C.2. Fluid picture
If we take Laguerre moments in µB = v ⊥ / et al. G , = n and G , = T ⊥ , and again taking the small- k ⊥ ρ limit, we have ∂n∂t + (cid:18) − b (cid:19) v E ∂n∂y − b v E ∂T ⊥ ∂y = 0 , (C 4) ∂T ⊥ ∂t + (cid:18) − b (cid:19) v E ∂T ⊥ ∂y − b v E ∂n∂y − b v E ∂G , ∂y = 0 , (C 5)where the G , Laguerre moment is not evolved but requires closure. These equations areidentical to the Beer equations in the 1D low k ⊥ ρ limit, with the exception of the last G , term in (C 5). Now that we have identified this extra term, we can generalize theseequations to 3D, and to higher k ⊥ ρ using the full FLR expressions (Dorland & Hammett1993): ∂n∂t + v Φ · ∇ n + (cid:20)
12 ˆ ∇ ⊥ v Φ (cid:21) · ∇ T ⊥ = 0 , (C 6) ∂T ⊥ ∂t + v Φ · ∇ T ⊥ + (cid:20)
12 ˆ ∇ ⊥ v Φ (cid:21) · ∇ n + (cid:104) ˆˆ ∇ ⊥ v Φ (cid:105) · ∇ T ⊥ + (cid:20)
12 ˆ ∇ ⊥ v Φ (cid:21) · ∇ G , = 0 . (C 7)C.3. Fluid closure
Now we must find a closure expression for the extra term involving G , in (C 7). If wesimply set G , = 0, as Beer did, we get an oscillatory (undamped) solution for n and T ⊥ . In the spirit of the closures pioneered by Hammett and Perkins to model Landau8 E. G. Highcock, N. R. Mandell, M. Barnes and W. Dorland damping (Hammett & Perkins 1990), we can instead use dissipative closures to producefluid density and temperature responses that are damped like the kinetic responses. Thuswe will choose a closure of the form G , = (cid:12)(cid:12)(cid:12) [ ˆ ∇ ⊥ v Φ ] · ∇ (cid:12)(cid:12)(cid:12) [ ˆ ∇ ⊥ v Φ ] · ∇ ( µ n + µ T ⊥ ) , (C 8)where we follow (Beer & Hammett 1996) and allow each coefficient µ to have a dissipativeand a reactive (non-dissipative) piece, given by µ = µ r + µ i (cid:12)(cid:12)(cid:12) [ ˆ ∇ ⊥ v Φ ] · ∇ (cid:12)(cid:12)(cid:12) [ ˆ ∇ ⊥ v Φ ] · ∇ = ( µ r , µ i ) . (C 9)We set the µ and µ coefficients by numerically minimizing the difference between thekinetic density response (C 3) and the fluid density response found by evolving the n and T ⊥ equations after inserting the closures. We have found that µ = (0 . , − . µ = (1 . , − . k ⊥ ρ (cid:46)
1. C.4.
Extension to other moment equations
One can follow a similar procedure to develop nonlinear phase mixing closures forthe remaining gyrofluid equations in the 4+2 model. The u (cid:107) and q ⊥ equations form acoupled system identical to the system we have studied above, so we can use the same2-moment closure and simply replace n and T ⊥ with u (cid:107) and q ⊥ , respectively. The T (cid:107) and q (cid:107) equations are also identical; each of these is not coupled to other equations, so theclosures that appear in these equations are 1-moment closures. The closure coefficient isfound with a similar procedure to the one used above, by fitting the fluid response for T (cid:107) to the corresponding kinetic response.This completes our derivation of a gyrofluid closure to model nonlinear FLR phasemixing. Adding these new terms to the original nonlinear terms, the final gyrofluidnonlinear terms used in Gryfx are given by N n = v Φ · ∇ n + (cid:20)
12 ˆ ∇ ⊥ v Φ (cid:21) · ∇ T ⊥ , (C 10) N u (cid:107) = v Φ · ∇ u (cid:107) + (cid:20)
12 ˆ ∇ ⊥ v Φ (cid:21) · ∇ q ⊥ , (C 11) N T (cid:107) = v Φ · ∇ T (cid:107) + (cid:12)(cid:12)(cid:12)(cid:12)(cid:20)
12 ˆ ∇ ⊥ v Φ (cid:21) · ∇ (cid:12)(cid:12)(cid:12)(cid:12) µ T (cid:107) , (C 12) N T ⊥ = v Φ · ∇ T ⊥ + (cid:20)
12 ˆ ∇ ⊥ v Φ (cid:21) · ∇ n + (cid:104) ˆˆ ∇ ⊥ v Φ (cid:105) · ∇ T ⊥ + (cid:12)(cid:12)(cid:12)(cid:12)(cid:20)
12 ˆ ∇ ⊥ v Φ (cid:21) · ∇ (cid:12)(cid:12)(cid:12)(cid:12) ( µ n + µ T ⊥ ) , (C 13) N q (cid:107) = v Φ · ∇ q (cid:107) + (cid:12)(cid:12)(cid:12)(cid:12)(cid:20)
12 ˆ ∇ ⊥ v Φ (cid:21) · ∇ (cid:12)(cid:12)(cid:12)(cid:12) µ q (cid:107) , (C 14) N q ⊥ = v Φ · ∇ q ⊥ + (cid:20)
12 ˆ ∇ ⊥ v Φ (cid:21) · ∇ u (cid:107) + (cid:104) ˆˆ ∇ ⊥ v Φ (cid:105) · ∇ q ⊥ + (cid:12)(cid:12)(cid:12)(cid:12)(cid:20)
12 ˆ ∇ ⊥ v Φ (cid:21) · ∇ (cid:12)(cid:12)(cid:12)(cid:12) (cid:0) µ u (cid:107) + µ q ⊥ (cid:1) , (C 15)where N m represents the nonlinear terms in the m moment equation, and theabsolute value terms comprise our new closure terms with closure coefficients ptimisation of fusion confinement µ = (0 . , − . , µ = (1 . , − . µ = (0 . , − . k ⊥ ρ in the n and u (cid:107) equations, which is beyondthe order of accuracy of the usual FLR terms. We avoid this by including closure termsproportional to n and u (cid:107) in the T ⊥ and q ⊥ equations, respectively.In the derivation of these new closures, we only considered nonlinear FLR phase mixingfrom a static, 1D potential. In an evolving 3D system, the phase mixing process is muchmore complicated, and our rough model may overestimate or underestimate the amountof phase mixing. Nonetheless, our model will capture the correct scaling of the mixing,producing physically motivated damping at large ϕ and high k ⊥ . In this way our newterms can be interpreted as a type of hyperviscosity, but one that damps at the ρ scaleas 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 P = | v ·∇| M . Since we use a Fourier spectral representation for the equations, we are interestedin the Fourier transform of P . Denoting the Fourier transform of a quantity x as ˆ x k , anddefining N = v · ∇ M , we evaluate the closure term asˆ P k = (cid:92) ( | v · ∇| M ) k = | ˆ N k | ˆ M k | ˆ M k | , (C 16)where all operations are performed in Fourier space, and ˆ N k can be calculated pseu-dospectrally 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 zonalflows. These nonlinearly-driven sheared poloidal E × B flows have been shown to playa key role in determining the turbulence saturation level. Therefore inaccuracy in zonalflow dynamics has been one of the main sources of disagreement between gyrofluid andgyrokinetic turbulence models. † Attempts were made, with limited success, to modify thegyrofluid closures (Beer & Hammett 1998) to capture the linearly undamped componentof the flows derived by Rosenbluth & Hinton (1998). These closure modifications hadlimited success (Dimits et al. et al. v (cid:107) direction resultingfrom trapped particle dynamics, which in a moment approach would require much higherHermite moments in v (cid:107) than used in the Beer model (Mandell et al. Gryfx : we evolvethe 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 byequations (B 4–B 9) above. Because we use a Fourier spectral representation for theperpendicular configuration space discretization in
Gryfx , and because the Fouriermodes only interact via the nonlinearity, we can easily choose an alternative algorithmfor the linear evolution of the zonal ( k y = 0) modes. In order to nonlinearly couple † This is not to say that zonal flow dynamics is the only source of disagreement. Our resultsindicate that even with accurate zonal flow dynamics, a model of nonlinear phase mixing isrequired to produce the agreement between the gyrofluid and gyrokinetic models shown inFigure 7. E. G. Highcock, N. R. Mandell, M. Barnes and W. Dorland the zonal and non-zonal modes, we must first be able to transform the gyrokineticdistribution function into gyrofluid moments, and vice versa. The transformation fromthe gyrokinetic distribution function to gyrofluid moments is simply the process of takingvelocity moments, which can also be interpreted as projecting the distribution functiononto a Laguerre-Hermite basis as shown in equation (B 3). For the inverse transformation,from gyrofluid moments to the gyrokinetic distribution function, we can expand thedistribution function in the Laguerre-Hermite basis as (Mandell et al. gF = (cid:88) (cid:96) =0 (cid:88) m =0 ( − (cid:96) √ m ! L (cid:96) ( µB ) He m ( v (cid:107) ) G (cid:96),m (D 1) ≈ n + v (cid:107) u (cid:107) + 12 (cid:16) v (cid:107) − (cid:17) T (cid:107) + ( µB − T ⊥ + 12 (cid:32) v (cid:107) − v (cid:107) (cid:33) q (cid:107) + v (cid:107) ( µB − q ⊥ (D 2)where in the second line we have explicitly expressed the expansion in terms of the Beergyrofluid moments and truncated.Following the same logic, we can expand the gyrokinetic nonlinear term in terms ofthe gyrofluid nonlinear terms given in (C 10–C 15): N gk = ˆ b × ∇ (cid:104) ϕ (cid:105) R · ∇ h = ˆ b × ∇ (cid:104) ϕ (cid:105) R · ∇ g ≈ N n + v (cid:107) N u (cid:107) + 12 (cid:16) v (cid:107) − (cid:17) N T (cid:107) + ( µB − N T ⊥ + 12 (cid:32) v (cid:107) − v (cid:107) (cid:33) N q (cid:107) + v (cid:107) ( µB − N q ⊥ , (D 3)Thus we can construct the gyrokinetic nonlinear term for the k y = 0 zonal modes fromthe k y = 0 component of the six gyrofluid nonlinear terms. The full hybrid algorithmthen proceeds as follows.(i) Calculate moments of the k y = 0 component of the gyrokinetic distributionfunction ( e.g. via (B 3)).(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).(iv) Evolve the k y (cid:54) = 0 modes with equations (B 4–B 9).(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 gyrokineticequation in step (v). Further, we do not include the nonlinear phase mixing closureterms in the k y = 0 component of the gyrofluid nonlinear terms. This is in part dueto the fact that introducing damping from closure models can suppress the zonal flowresidual (which is precisely the reason we moved to gyrokinetic evolution of the zonalmodes). Thus in our model, the nonlinear drive for the zonal flows comes only from thelowest several moments of the distribution function; this approximation is justified bythe results of (Rogers et al. Gryfx we takeadvantage of GPU-CPU concurrency by simultaneously evolving the non-zonal gyrofluidequations on the GPU and the zonal gyrokinetic equation on the CPU. Using a singleGPU and ∼
16 CPU cores (a common supercomputer node configuration), the two stepstake roughly the same amount of wall clock time. This means that in our scheme the ptimisation of fusion confinement
Appendix E. Practical Considerations
The presentation of equations and the methodology given above has, for reasons ofclarity and brevity, skipped over some of the most thorny issues encountered during thecurrent study. However, we believe that a discussion of them should be present in thiswork, and we hope that by including it here as an appendix, we can save at least sometime and effort for future researchers. In this discussion it helps to think of the operationof
CoRFu as a series of six nested loops, as illustrated in Fig. 8. Within each of theseloops issues can arise. The principle practical challenges encountered in the constructionand use of the
CoRFu framework (independent of the challenges in constructing itsindividual 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.E.1.
Incorporating evolution of the magnetic equilibrium
The
Trinity transport solver has the ability to evolve the magnetic equilibriuminternally. However, at present it is not able to treat the evolution of the magneticequilibrium implicitly, as it does the profiles of density, temperature and flow. Thereforeit 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 andmagnetic equilibrium solutions. Thus, provided
Trinity discovers a solution in whichboth 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 Fig. 8. In effect,
Trinity evolvesusing a fixed magnetic equilibrium until the profiles stop changing. At this point the newpressure profile is passed to
Chease which recalculates the equilibrium.
Trinity is thenre-run with new equilibrium and the cycle is repeated until the optimisation criterionconverges to within a specified tolerance (see Appendix E.2).E.2.
Determining whether a particular turbulence calculation had reached steady state
This is an extensive subject fraught with complication. Not only is it hard to automatewhat enters into human judgement in these cases (in many cases the best way ofdetermining 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 caseof slow zonal flow amplitude rise in ETG turbulence; Colyer et al. (2017)).The difficulty proceeds from having no a priori knowledge of the timescales of a givensituation: in particular, those of the linear growth, time to saturation, and the longesttimescales in the saturated state (usually related to zonal flows). A typical dilemma mightbe: is the heat flux no longer changing because(i) the simulation has reached a saturated nonlinear state, or2
E. G. Highcock, N. R. Mandell, M. Barnes and W. Dorland
Figure 8.
Schematic of
CoRFu operation showing the hierarchy of decision loops which mustbe walked through. ptimisation of fusion confinement O (10) turnover times ( a/v thi ). The turbulencecode driver then ran the calculation for two initial stages, and then proceeded to runeach calculation according to the decision tree presented in Fig. 9. This decision tree isextremely restrictive about allowing a simulation to be declared converged. This choicewas motivated by the observation that a bad value for the heat flux incurred such agreat 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.E.3.
Encountering failures of either Newton iteration, the turbulence code or theGS solver
It is of course possible for any component of
CoRFu to experience a failure duringthe very long operation typically required. This could be due to hardware failures, filecorruptions, and so on. However, the three most common failures that were encounteredwere • failure of the Newton iteration to converge due to overlarge timestep or a bad fluxvalue, • numerical instability from an overly large timestep in Gryfx , and • failure of Chease to converge on a solution, usually as a result of a non-monotonicpressure 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 inherentnoise in the turbulent fluxes. In particular, this was a problem if a large change in agradient (e.g. the temperature gradient κ T = R/L T ) produced only a small change ina transport coefficient (e.g. the heat flux Q ). Since we are inverting a Jacobian whichessentially depends on quantities such as ∆Q/∆κ T , if the true value of ∆Q is small, andthere is noise, of the same magnitude, which contrives to give ∆Q close to zero, there isan artificial stationary point created in the search domain which will lead the Newtoniteration to take a large jump incorrectly. In the future, it is planned to implement aLevenberg-Marquardt mixed search which will be more resistant to this; for the currentstudy it was necessary to be very stringent with the turbulence convergence conditions,and to use a large, fixed value of ∆κ T = 0 .
6. The reasoning behind this was that since Q ∝ κ T (Barnes et al. b ), a constant value of ∆κ T would continue to produce arespectable ∆Q as Q tended to 0, but would not produce an excessively large ∆Q farfrom the threshold.With this implemented, a remaining likely cause of a failure was a real stationary pointin the search domain, and the solution was to make the transport timestep smaller. Inthe next section we will discuss this in more detail with regard to efficiency. However,there was another consideration when resetting the timestep: what to do about the stateof the turbulence calculation.If the Newton iteration had deviated very far from the previous profiles it might be4 E. G. Highcock, N. R. Mandell, M. Barnes and W. Dorland
Figure 9.
Schematic of the
CoRFu turbulence code driver showing the logic used to determineif a turbulence calculation is saturated ( Q o is the heat flux from the previous iteration, Q p isthe heat flux from the previous stage of the turbulence calculation, and Q c is the heat flux fromthe current stage). reasonable to assume that since the turbulent state at the end of the failed iteration wasvery far from that at the beginning, it would be better to start the turbulence calculationfrom noise to avoid the turbulence displaying violent transients (particularly if the endof the failed iteration was in a strongly-driven regime). In fact, since (because of stifftransport) the steady-state profiles were never very far from marginal (with low growthrates), such a strategy would incur a large time penalty and run the risk of mistaking alow-growing turbulent case for a stable case.In the case of Gryfx , a timestep was chosen that was sufficiently small for almost ptimisation of fusion confinement
Gryfx simulationstage using the state of the previous simulation stage (since either it was continuing acalculation till convergence, or it was starting a new calculation with similar parametersto the previous); however this behaviour was changed if necessary (note that resetting theinitial condition close to threshold when the growth rate was low incurred a significanttime penalty). Otherwise, the following steps were taken to mitigate problems: • If the previous
Trinity
Newton iteration did not converge, values of fields andmoments were kept but time averages of fluxes were reset. • If the previous fluxes were very small ( < − ), all values were reset. • If the previous fluxes were either very large ( > e
10) or had
NaN values, all valueswere reset. Typically this was because a quirk of the initial condition pushed theturbulence amplitudes to values too large for the timestep and could be fixed by resetting,but if this occurred repeatedly the optimisation would be halted and require intervention. • If fewer than 3 turbulence stages had been run the moving averages would be resetto eliminate initial transients. • If greater than 3 turbulence stages had been run, the timescale of the moving averagewould gradually increase in successive stages (proportional to the total turbulence timeelapsed in these successive stages) to force eventual convergence of the heat flux. Thiswas a necessary to ensure that the transport calculation never became stuck, and had noimpact 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 anunconverged pressure profile,
CoRFu first attempted to run
Trinity for an extra periodof time using the previous equilibrium. Occasionally in preliminary studies (though notin the study reported here) the GS solver failed because the control parameters strayedinto a region of parameter space where a solution could not be found. In such cases ahighly unfavourable value was returned to prevent the algorithm pursuing that searchdirection. E.4.
Dealing with turbulent thresholds
To a user of the
CoRFu framework, it is an inconvenient reality that optimal solutionstypically 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. Thisis a problem because of the extra time required to calculate the turbulent transport. Ifthe previous transport iteration was below the threshold and the subsequent iterationis likewise, the turbulence calculation will end rapidly because any initial noise willdie 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 theprevious iteration was below threshold and the subsequent one above, (i.e the transportsolver crosses the zero-turbulence manifold from below) the turbulence calculation willnot converge until amplitudes have grown linearly to saturation amplitudes. Since weare close to threshold with small growth rates, this may take a very long time. Equallyunpleasant is when the zero-turbulence manifold is crossed from above: in this case thefree energy of the previous vigorous turbulence must be absorbed by only weakly dampedmodes.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. κ T ). Effectively we needed to make sure that, as far as possible, the lower part of thederivative 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 ofreducing the convergence rate of the Newton iteration and hence the size of the timestep).6 E. G. Highcock, N. R. Mandell, M. Barnes and W. Dorland
Fortunately, this requirement coincided with the requirement described in Appendix E.3for a large ∆Q . The second way to achieve this is to set a tight convergence criterionfor the Newton iteration to make sure that it rarely artificially (that is, at odds withthe 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 numberof turbulence calculations carried out, to strictly control the cost of the optimisationstudy. For many reasons, as given above, it is necessary to be overly cautious and runthe turbulence simulations until a strict set of criteria have been satisfied. Nonetheless, itwill always be the case that during every iteration of the transport calculation some fluxcalculations will converge before others. Therefore the flux driver is set up so that theoptimisation can run on an arbitrary number of GPUs, and that if the number of GPUsis less than the number of flux calculations carried out at a time, then at each stage of theflux calculation the GPUs will cycle through each of the flux calculations requested by thedriver. As the number of stages increases, more and more of the turbulence calculationswill 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 fourGPUs, this arrangement reduced the simulation cost by close to 50%.An equally important consideration with regard to cost is the choices made for iterationand time stepping parameters within
Trinity . When seeking a steady state solution, acertain amount of time will need to elapse while the profiles adjust to the new parametersprovided by the optimisation driver. A first assumption might be that increasing the sizeof the timestep will reduce the cost of evolving the profiles for that time (rememberingthat the
Trinity algorithm is implicit and therefore stable for large timesteps). However,the larger the timestep, the greater the difference between the profiles at the end andthe profiles at the beginning of the timestep, and thus the greater the chance that theNewton search will not be able to land upon the future profiles. If this occurs, the timestepmust be reduced and the step repeated, effectively doubling the cost of that timestep.Increasing the maximum number of steps in the Newton iteration would increase thechances of it reaching a solution; however, it would also increase the cost should theiteration still fail.By trial and error, the following choices were found to enable efficient and robustoperation. • The threshold for rejecting a Newton iteration was set high ( errtol =0.1) to avoidthe cost of repeating an iteration (this value was still low enough to avoid wacky iterationstriggered by a stationary point). • The threshold for saying that the Newton iteration was converged and thereforefewer than the allowed number of steps were needed was set low ( errflr =0.02) so thatthe iteration was always as accurate as possible, leading to an accurate time evolution. • The threshold for increasing the timestep was set very low ( flrfac =4 so that thethreshold was 0.02/4) so that the timestep was only increased if the iterations wereproceeding very smoothly. • The factor by which the timestep was reduced after a failed iteration was set high( deltadj =4) to avoid repeated failed iterations. For example, when crossing a turbulentthreshold it was often necessary to rapidly reduce the transport timestep: if deltadj wassmall, it would take a long time before the timestep was reduced sufficiently. • The maximum allowed number of steps in the Newton iteration was set to 4. This ptimisation of fusion confinement
Choosing resolution parameters for the turbulence code
The resolution required for the turbulence calculation, including perpendicular boxsize, perpendicular and parallel resolutions, and timestep, vary significantly with thevalues of the driving parameters and the magnetic geometry.In the present case, the GPUs that were available for the calculation were sufficientlypowerful to allow the increase of the
Gryfx resolution to a size that was more thanadequate for the majority of simulations (as determined by manual inspection of a sampleof turbulence simulations), without an unmanageable increase in cost.There are a posteriori ways of checking that resolutions were sufficient, such asexamining 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 donot change with increased resolution. However, the automation of such investigations,which are usually carried out by hand for individual turbulence calculations, has not yetbeen implemented in the
CoRFu framework, and would constitute significant additionalresearch.
REFERENCESAbel, I. G. & Cowley, S. C.
New Journal of Physics , arXiv: 1210.1417. Abel, I. G., Plunk, G. G., Wang, E., Barnes, M., Cowley, S. C., Dorland, W.& Schekochihin, A. A.
Reports on progress in physics. Physical Society(Great Britain) (11), 116201. Abramson, D., Peachey, T. & Lewis, A.
Proceedings of the 6th international conference on ComputationalScience (ICCS’06) , 720–727. Barnes, M., Abel, I. G., Dorland, W., Gorler, T., Hammett, G. W. & Jenko, F.
Physics ofPlasmas (5), 056109, arXiv: 0901.2868. Barnes, M., Parra, F. I., Highcock, E. G., Schekochihin, A., Cowley, S. C. & Roach,C. M. a Turbulent Transport in Tokamak Plasmas with Rotational Shear.
PhysicalReview Letters (17), 1–4.
Barnes, M., Parra, F. I. & Schekochihin, a. a. b Critically Balanced Ion TemperatureGradient Turbulence in Fusion Plasmas.
Physical Review Letters (11), 115003.
Beer, M. & Hammett, G.
VarennaProceedings pp. 1–11.
Beer, M. A., Cowley, S. C. & Hammett, G. W.
Physics of Plasmas (7), 2687–2700. Beer, M. A. & Hammett, G. W.
Physics of Plasmas (11), 4046. Bourdelle, C., Citrin, J., Baiocchi, B., Casati, A., Cottier, P., Garbet, X. &Imbeaux, F.
Plasma Physics and Controlled Fusion (1), 014036. Budny, R. V.
Nuclear Fusion (8), 85008. Burrell, K. H.
Physics of Plasmas (5), 1499–1518. Chang, C. S. & Hinton, F. L.
Physics of Fluids (1493), 3314. E. G. Highcock, N. R. Mandell, M. Barnes and W. Dorland
Chen, C., Wicks, R., Horbury, T. & Schekochihin, A.
The Astrophysical Journal Letters (2), L79.
Citrin, J., Breton, S., Felici, F., Imbeaux, F., Aniel, T., Artaud, J. F., Baiocchi, B.,Bourdelle, C., Camenen, Y. & Garcia, J. a Real-time capable first principlebased modelling of tokamak turbulent transport.
Nuclear Fusion (9), 92001, arXiv:1502.07402v1. Citrin, J., Garcia, J., G¨orler, T., Jenko, F., Mantica, P., Told, D., Bourdelle, C.,Hatch, D. R., Hogeweij, G. M. D., Johnson, T., Pueschel, M. J. & Schneider,M. b Electromagnetic stabilization of tokamak microturbulence in a high- β regime. Plasma Physics and Controlled Fusion (1), 014032. Citrin, J., Jenko, F., Mantica, P., Told, D., Bourdelle, C., Dumont, R., Garcia, J.,Haverkort, J., Hogeweij, G., Johnson, T. & Pueschel, M.
Nuclear Fusion (2), 023008. Colyer, G. J., Schekochihin, A. A., Parra, F. I., Roach, C. M., Barnes, M. A., Ghim,Y. C. & Dorland, W.
Plasma Physics and Controlled Fusion (5), arXiv: 1607.06752. Dimits, A. M., Bateman, G., Beer, M. A., Cohen, B. I., Dorland, W., Hammett,G. W., Kim, C., Kinsey, J. E., Kotschenreuther, M., Kritz, a. H., Lao, L. L.,Mandrekas, J., Nevins, W. M., Parker, S. E., Redd, a. J., Shumaker, D. E.,Sydora, R. & Weiland, J.
Physics of Plasmas (3), 969. Dorland, W. & Hammett, G. W.
Physics of Fluids B-Plasma Physics , 812–835. Dorland, W., Jenko, F., Kotschenreuther, M. & Rogers, B. N.
Physical review letters (26 Pt 1), 5579–5582. EJ Doyle, WA Houlberg, Kamada, Y., V Mukhovatov, TH Osborne, A Polevoi,Bateman, G., JW Connor, JG Cordey, T Fujita, X Garbet, TS Hahm, LDHorton, AE Hubbard, F Imbeaux, F Jenko, J. E. Kinsey, Y Kishimoto, J Li, TCLuce, Y Martin, M Ossipenko, V Parail, A Peeters, TL Rhodes, JE Rice, Roach,C. M., Rozhansky, V., Ryter, F., G Saibene, R Sartori, ACC Sips, JA Snipes, MSugihara, EJ Synakowski, H Takenaga, T Takizuka, K Thomsen, MR Wade, HRWilson, ITPA Transport Physics Topical Group, ITPA Confinement Database,Modelling Topical Group, ITPA Pedestal, Edge Topical Group, Physics), E.J. D. C. T., Database, W. A. H. C. C., Pedestal, Y. K. C., Physics), V. M. c.-C. T., Pedestal, T. H. O. c.-C., Database, A. P. c.-C. C., Bateman, G., Connor,J. W., (retired), J. G. C., Fujita, T., Garbet, X., Hahm, T. S., Horton, L. D.,Hubbard, A. E., Imbeaux, F., Jenko, F., Kinsey, J. E., Kishimoto, Y., Li, J., Luce,T. C., Martin, Y., Ossipenko, M., Parail, V., Peeters, A., Rhodes, T. L., Rice,J. E., Roach, C. M., Rozhansky, V., Ryter, F., Saibene, G., Sartori, R., Sips, A.C. C., Snipes, J. A., Sugihara, M., Synakowski, E. J., Takenaga, H., Takizuka,T., Thomsen, K., Wade, M. R., Wilson, H. R., Group, I. T. P. T., Database,I. C., Group, M. T., Pedestal, I., Group, E. T., Doyle, E. J., Houlberg, W. A.,Kamada, Y., Mukhovatov, V., Osborne, T. H., Polevoi, A. & Cordey, J. G.
Nuclear Fusion (6), S18. ESTECO, M.
Federici, G., Kemp, R., Ward, D., Bachmann, C., Franke, T., Gonzalez, S., Lowry,C., Gadomska, M., Harman, J., Meszaros, B., Morlock, C., Romanelli, F.& Wenninger, R.
FusionEngineering and Design (7-8), 882–889. Frieman, E. A. & Chen, L.
Physics of Fluids (3), 502. Galambos, J., Perkins, L., Haney, S. & Mandrekas, J.
Nuclear Fusion (5), 551–573. Glasstone, S. & Lovberg, R. H.
Hammett, G. W. & Perkins, F. W. ptimisation of fusion confinement application to the ion-temperature-gradient instability. Physical review letters (25),3019–3022. Highcock, E. G., Barnes, M., Parra, F. I., Schekochihin, A. A., Roach, C. M. &Cowley, S. C.
Physics of Plasmas (10), 102304. Highcock, E. G., Barnes, M., Schekochihin, A. A., Parra, F. I., Roach, C. M. &Cowley, S. C.
PhysicalReview Letters (21), 215003.
Highcock, E. G., Schekochihin, A. A., Cowley, S. C., Barnes, M., Parra, F. I., Roach,C. M. & Dorland, W.
PhysicalReview Letters (26), 265001.
Howes, G. G., TenBarge, J. M., Dorland, W., Quataert, E., Schekochihin, A. A.,Numata, R. & Tatsuno, T.
Physical review letters (3), 035004.
Jardin, S. C., Kessel, C. E., Mau, T. K., Miller, R. L., Najmabadi, F., Chan, V. S.,Chu, M. S., Lahaye, R., Lao, L. L., Petrie, T. W., Politzer, P., St. John, H. E.,Snyder, P., Staebler, G. M., Turnbull, A. D. & West, W. P.
Fusion Engineering and Design (1-4), 25–62. Kawamori, E.
Physical review letters (9), 095001.
Kinsey, J. E., Staebler, G. M., Candy, J., Waltz, R. E. & Budny, R. V.
Nuclear Fusion (8), 083001. Kotschenreuther, M., Rewoldt, G. & Tang, W. M.
Computer PhysicsCommunications (2-3), 128–140. Lewis, A.
Luce, T., Challis, C., Ide, S., Joffrin, E., Kamada, Y., Politzer, P., Schweinzer,J., a.C.C. Sips, Stober, J., Giruzzi, G., Kessel, C., Murakami, M., Na, Y.-S.,Park, J., a.R. Polevoi, Budny, R., Citrin, J., Garcia, J., Hayashi, N., Hobirk,J., Hudson, B., Imbeaux, F., Isayama, a., McDonald, D., Nakano, T., Oyama, N.,Parail, V., Petrie, T., Petty, C., Suzuki, T. & Wade, M.
Nuclear Fusion (1), 013015. L¨utjens, H., Bondeson, A. & Sauter, O.
Computer physics communications , 219–260. Maggi, C. F., Groebner, R. J., Oyama, N., Sartori, R., Horton, L. D., Sips, A. C.,Suttrop, W., Leonard, A., Luce, T. C., Wade, M. R., Kamada, Y., Urano, H.,Andrew, Y., Giroud, C., Joffrin, E. & De La Luna, E.
Nuclear Fusion (7), 535–551. Mandell, N. & Dorland, W.
Bulletin of the American Physical Society , p. Abstract CP8.039.
Mandell, N., Dorland, W. & Landreman, M.
Journal of Plasma Physics (1). Manousopoulos, P. & Michalopoulos, M.
European Journal of Operational Research (2),594–602.
Marinoni, A., Brunner, S., Camenen, Y., Coda, S., Graves, J. P., Lapillonne, X.,Pochelon, A., Sauter, O. & Villard, L.
Plasma Physics and Controlled Fusion (5), 055016. Meneghini, O., Snyder, P. B., Smith, S. P., Candy, J., Staebler, G. M., Belli, E. A.,Lao, L. L., Park, J. M., Green, D. L., Elwasif, W., Grierson, B. A. & Holland,C.
Physicsof Plasmas (4), 042507. Merle, A., Sauter, O. & Yu Medvedev, S. E. G. Highcock, N. R. Mandell, M. Barnes and W. Dorland negative triangularity using the EPED-CH model.
Plasma Physics and Controlled Fusion (10). Mukhovatov, V., Shimomura, Y., Polevoi, A., Shimada, M., Sugihara, M., Bateman,G., Cordey, J. G., Kardaun, O., Pereverzev, G., Voitsekhovitch, I., Weiland,J., Zolotukhin, O., Chudnovskiy, A., Kritz, A. H., Kukushkin, A., Onjun, T.,Pankin, A. & Perkins, F. W.
Nuclear Fusion (9), 942–948. Parail, V., Albanese, R., Ambrosino, R., Artaud, J. F., Besseghir, K., Cavinato, M.,Corrigan, G., Garcia, J., Garzotti, L., Gribov, Y., Imbeaux, F., Koechl, F.,Labate, C. V., Lister, J., Litaudon, X., Loarte, A., Maget, P., Mattei, M.,McDonald, D., Nardon, E., Saibene, G., Sartori, R. & Urban, J.
Nuclear Fusion (11), 113002. Parker, J. T., Highcock, E. G., Schekochihin, A. A. & Dellar, P. J.
Physics of Plasmas (7), 070703,arXiv: 1603.06968. Pochelon, A., Angelino, P., Behn, R., Brunner, S., Coda, S., Kirneva, N., Medvedev,S., Reimerdes, H., Rossel, J., Sauter, O., Villard, L., W´agner, D., Bottino,A., Camenen, Y., Canal, G. P., Chattopadhyay, P. K., Duval, B. P., Fasoli, A.,Goodman, T. P., Jolliet, S., Karpushov, A., Labit, B., Marinoni, A., Moret,J. M., Pitzschke, A., Porte, L., Rancic, M. & Udintsev, V. S.
Plasma andFusion Research (SPL.ISS.1), 1–8. Pochelon, A., Goodman, T., Henderson, M., Angioni, C., Behn, R., Coda, S.,Hofmann, F., Hogge, J.-P., Kirneva, N., Martynov, A., Moret, J.-M., Pietrzyk,Z. A., Porcelli, F., Reimerdes, H., Rommers, J., Rossi, E., Sauter, O., Tran,M. Q., Weisen, H., Alberti, S., Barry, S., Blanchard, P., Bosshard, P., Chavan,R., Duval, B. P., Esipchuck, Y. V., Fasel, D., Favre, A., Franke, S., Furno,I., Gorgerat, P., Isoz, P.-F., Joye, B., Lister, J., Llobet, X., Magnin, J.-C.,Mandrin, P., Manini, A., Marl´etaz, B., Marmillod, P., Martin, Y., Mayor, J.-M., Mlynar, J., Nieswand, C., Paris, P., Perez, A., Pitts, R., Razumova, K.,Refke, A., Scavino, E., Sushkov, A., Tonetti, G., Troyon, F., Toledo, W. V.& Vyas, P.
Nuclear Fusion (11Y), 1807–1818. Roach, C., Walters, M., Budny, R., Imbeaux, F., Fredian, T., Greenwald, M.,Stillerman, J., Alexander, D., Carlsson, J., Cary, J., Ryter, F., Stober, J.,Gohil, P., Greenfield, C., Murakami, M., Bracco, G., Esposito, B., Romanelli,M., Parail, V., Stubberfield, P., Voitsekhovitch, I., Brickley, C., Field, A.,Sakamoto, Y., Fujita, T., Fukuda, T., Hayashi, N., Hogeweij, G., Chudnovskiy,A., Kinerva, N., Kessel, C., Aniel, T., Hoang, G., Ongena, J., Doyle, E.,Houlberg, W. & Polevoi, A.
Nuclear Fusion (12), 125001. Rogers, B. N., Dorland, W. & Kotschenreuther, M.
Physical review letters (25),5336–9. Rosenbluth, M. & Hinton, F.
Physical review letters (4), 724. Schekochihin, A. A., Cowley, S. C., Dorland, W., Hammett, G. W., Howes, G. G.,Quataert, E. & Tatsuno, T.
The Astrophysical JournalSupplement Series (1), 310–377.
Schekochihin, A. A., Parker, J. T., Highcock, E. G., Dellar, P. J., Dorland, W. &Hammett, G. W.
Journal of Plasma Physics (02), 905820212, arXiv: 1508.05988. Simpson, T., Toropov, V., Balabanov, V. & Viana, F. (September), 1–22. ptimisation of fusion confinement Snyder, P. & Hammett, G.
Physics of Plasmas (7), 3199–3216. Sorbom, B. N., Ball, J., Palmer, T. R., Mangiarotti, F. J., Sierchio, J. M., Bonoli,P., Kasten, C., Sutherland, D. A., Barnard, H. S., Haakonsen, C. B., Goh, J.,Sung, C. & Whyte, D. G.
Fusion Engineeringand Design , 378–405, arXiv: 1409.3540.
Staebler, G. M. & John, H. E. S.
Nuclear Fusion (8), L6–L8. Staebler, G. M., Kinsey, J. E. & Waltz, R. E.
Physics of Plasmas (5). Stork, D., Agostini, P., Boutard, J. L., Buckthorpe, D., Diegele, E., Dudarev, S. L.,English, C., Federici, G., Gilbert, M. R., Gonzalez, S., Ibarra, A., Linsmeier,C., Puma, A. L., Marbach, G., Packer, L. W., Raj, B., Rieth, M., Tran, M. Q.,Ward, D. J. & Zinkle, S. J.
Fusion Engineeringand Design (7-8), 1586–1594. Sugama, H. & Horton, W.
Physics of Plasmas (7), 2560. Synakowski, E. J.
Plasma Physics and Controlled Fusion (5), 581–596. Tatsuno, T., Dorland, W., Schekochihin, A. A., Plunk, G. G., Barnes, M., Cowley,S. C. & Howes, G. G.
Physical Review Letters (1), 2–5, arXiv: 0811.2538.
Tatsuno, T., Plunk, G., Barnes, M., Dorland, W., Howes, G. & Numata, R.
Physics of Plasmas (12), 122305. Uckan, N.
ITER Documentation Series 10(Vienna: IAEA) . Wagner, F., Fussmann, G., Grave, T., Keilhacker, M., Kornherr, M., Lackner,K., McCormick, K., M¨uller, E. R., St¨abler, A., Becker, G., Bernhardi, K.,Ditte, U., Eberhagen, A., Gehre, O., Gernhardt, J., Gierke, G. v., Glock,E., Gruber, O., Haas, G., Hesse, M., Janeschitz, G., Karger, F., Kissel, S.,Kl¨uber, O., Lisitano, G., Mayer, H. M., Meisel, D., Mertens, V., Murmann, H.,Poschenrieder, W., Rapp, H., R¨ohr, H., Ryter, F., Schneider, F., Siller, G.,Smeulders, P., S¨oldner, F., Speth, E., Steuer, K. H., Szymanski, Z. & Vollmer,O.
Physical Review Letters (15), 1453–1456. Weisen, H., Alberti, S., Berry, S., Behn, R., Blanchard, P., Bosshard, P., B¨uhlmann,F., Chavan, R., Coda, S., Deschenaux, C., Dutch, M. J., Duval, B. P., Fasel, D.,Favre, A., Franke, S., Furno, I., Goodman, T., Henderson, M., Hofmann, F.,Hogge, J.-P., Isoz, P.-F., Joye, B., Lister, J. B., Llobet, X., Magnin, J.-C.,Mandrin, P., Marletaz, B., Marmillod, P., Martin, Y., Mayor, J.-M., Moret,J.-M., Nieswand, C., Paris, P., Perez, A., Pietrzyk, Z. a., Piffl, V., Pitts, R. a.,Pochelon, A., Razumova, K., Reimerdes, H., Refke, A., Rommers, J., Roy, I.,Sauter, O., Suttrop, W., Toledo, W. V., Tonetti, G., Tran, M. Q., Troyon,F., Vyas, P. & Ward, D. J.
Plasma Physics and Controlled Fusion (12B), B135–B144. Wenninger, R., Arbeiter, F., Aubert, J., Aho-Mantila, L., Albanese, R., Ambrosino,R., Angioni, C., Artaud, J.-F., Bernert, M., Fable, E., Fasoli, A., Federici, G.,Garcia, J., Giruzzi, G., Jenko, F., Maget, P., Mattei, M., Maviglia, F., Poli, E.,Ramogida, G., Reux, C., Schneider, M., Sieglin, B., Villone, F., Wischmeier, M.& Zohm, H.
NuclearFusion , 063003. White, A. E., Howard, N. T., Greenwald, M., Reinke, M. L., Sung, C., Baek, S.,Barnes, M., Candy, J., Dominguez, A., Ernst, D., Gao, C., Hubbard, a. E.,Hughes, J. W., Lin, Y., Mikkelsen, D., Parra, F. I., Porkolab, M., Rice, J. E.,Walk, J., Wukitch, S. J. & Team, A. C.-M. E. G. Highcock, N. R. Mandell, M. Barnes and W. Dorland at Alcator C-Mod and comparison with gyrokinetic simulations.
Physics of Plasmas (5), 056106. van Wyk, F., Highcock, E. G., Schekochihin, A. A., Roach, C. M., Field, A. R. &Dorland, W. Journal ofPlasma Physics (6), arXiv: 1607.08173. Xanthopoulos, P., Mynick, H. E., Helander, P., Turkin, Y., Plunk, G. G., Jenko, F.,G¨orler, T., Told, D., Bird, T. & Proll, J. H.
Physical Review Letters (15), 1–4.
Xiao, Y., Catto, P. J. & Dorland, W.
Physics of plasmas (5), 055910. Zohm, H., Angioni, C., Fable, E., Federici, G., Gantenbein, G., Hartmann, T.,Lackner, K., Poli, E., Porte, L., Sauter, O., Tardini, G., Ward, D. &Wischmeier, M.
Nuclear Fusion53