{'fraction_non_alphanumeric': 0.027080136746614983, 'fraction_numerical': 0.011359659748580043, 'mean_word_length': 4.865172087148721, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Direct multiscale coupling of a transport code to gyrokinetic turbulence codes

Direct coupling between a transport solver and local, nonlinear gyrokinetic calculations using the multiscale gyrokinetic code TRINITY [M. Barnes, Ph.D. thesis, arxiv:0901.2868] is described. The coupling of the microscopic and macroscopic physics is done within the framework of multiscale gyrokinetic theory, of which we present the assumptions and key results. An assumption of scale separation in space and time allows for the simulation of turbulence in small regions of the space-time grid, which are embedded in a coarse grid on which the transport equations are implicitly evolved. This leads to a reduction in computational expense of several orders of magnitude, making first-principles simulations of the full fusion device volume over the confinement time feasible on current computing resources. Numerical results from TRINITY simulations are presented and compared with experimental data from JET and ASDEX Upgrade plasmas.


I. INTRODUCTION
A fundamental challenge of fusion science is to maximize fusion power, which is determined primarily by macroscopic profiles of density and temperature. These profiles, which vary spatially on the system scale and evolve on the energy confinement time scale, drive turbulence at micro-scales in space and time. In the absence of MHD instability, this micro-turbulence is the dominant source of heat flux observed in standard tokamaks, which sets rigid constraints on the macroscopic profiles. 1 Consequently, it is of critical importance to understand the self-consistent interaction between the macroscopic profiles and the micro-turbulence. This is a challenging task due to both the wide range of scales involved and the high dimensionality of the system. The electron turbulence space scale, which is comparable to the electron Larmor radius, will be on the order of 0.1 millimeters for a device such as ITER, 2 which has a minor radius of about two meters. Similarly, the electron turbulence time scale is approximately one microsecond, much smaller than the expected energy confinement time of several seconds. Additionally, the instabilities driving the turbulence are kinetic in nature, requiring treatment of the velocity space in addition to the configuration space. Including all of these dynamics directly in a single simulation is not feasible on current computing resources.
As long as all relevant dynamics occur on time scales long compared to particle gyro-motion, it is possible to average over the gyro-orbits and eliminate the gyro-angle as a phase space variable. Furthermore, for magnetic confinement devices with sufficiently small ρ * (ratio of ion Larmor radius to plasma minor radius), there is expected to be a separation between micro-and macro-scales in space and time. 3,4 The δf -gyrokinetic model 5,6,7 exploits these scale separations, simplifying the system considerably. Given a set of fixed macroscopic profiles, it allows for the calculation of turbulent fluxes in a 5-D phase space. Such calculations have been performed numeri-cally in gyrokinetic codes for more than two decades, providing much insight into the nature of kinetic instabilities and micro-turbulence. First-principles δf -gyrokinetic simulations have steadily advanced in their sophistication and physical fidelity and have become routine in recent years.
However, these codes only model the effect of macroscopic profiles on micro-turbulence: They do not provide quantitative information on how this turbulence subsequently affects the evolution of the macroscopic profiles. Accurately calculating fluxes from experimental profiles is also problematic because of the acute sensitivity of the fluxes to small changes in the input profile gradients (which would arise due to uncertainty in experimental measurements). 8 In an attempt to address these issues, full-f gyrokinetic codes have been developed, which do not explicitly assume the scale separations listed above. Subsequently, they are able to take into account the twoway interaction between turbulence and equilibrium thermodynamic profiles. However, well-resolved full-f simulations would be extremely expensive numerically because of the wide range of scales described above and because of the necessity of calculating the distribution function and fields to very high order. It has also recently been argued that this approach as currently formulated leads to an unphysical source of toroidal angular momentum. 9 An alternative approach commonly used to study the interaction of the micro-and macro-physics is to solve fluid transport equations with a reduced model for the turbulent fluxes. Typically, the transport is modeled as a diffusive process, with turbulent diffusivities coming from a wide range of models, including empirical fits to experiment or simulation and theory-based estimates. 1,10,11,12,13,14 These reduced models have provided a basic qualitative understanding of the multiscale interaction and are capable in some cases of giving good quantitative agreement with first-principles, nonlinear gyrokinetic calculations. 15,16 However, such models do not permit detailed validation studies because they do not produce fluctuation spectra or other data that can be experimentally challenged. Without careful validation, even first-principles reduced models (with no adjustable parameters) cannot be fully trusted for predictions of performance in new operational regimes, and reduced models with fit parameters adjusted for current experiments have even less credibility. There are known cases for which even first-principles reduced models currently available are inadequate. 17,18 Consequently, it is desirable to couple not only reduced flux models, but also direct numerical simulations of the turbulence, to transport solvers. This has been done in Ref. 19 for a Hasegawa-Wakatani 2D fluid model of the turbulence using an implicit relaxation technique and allowing for nonlocality. Here, we describe coupling to local, nonlinear, 5D gyrokinetic turbulence calculations using a Newton method (similar to that described in Ref. 20), which accelerates convergence by more than an order of magnitude for typical parameters. This coupling is achieved using the multiscale gyrokinetic code, TRINITY, 21 which can use nonlinear fluxes from the continuum gyrokinetic codes GS2 22 or GENE. 23,24 Our approach is similar to that employed in TGYRO, 25 with the key distinction that TRINITY evolves the macroscopic profiles in time, whereas TGYRO assumes a steady-state and solves the volume-integrated transport equations for profile gradients.
It should be noted that some meso-scale phenomena in space and time are not formally considered in the standard multiscale gyrokinetic model. Our simulations ignore low-order magnetic islands and so are directly applicable only during MHD-quiescent periods of plasmas. Rapid cold/heat pulse propagation (such as following a sawtooth or ELM crash) is possible in the TRINITY code because of the presence of stiff critical gradients in ITG and TEM turbulence, though the transport time step would of course have to be reduced in TRINITY during such a transient event to be able to follow its propagation. It should be noted that flux-tube simulations include the contribution to the heat flux of avalanches on all scales up to the radial size of the flux tubes. If even longer wavelength avalanches were important, then the heat flux would increase as the flux-tube simulation domain was made larger, so convergence studies can be used to test this. PIC and continuum flux-tube simulations generally find that the flux converges with sufficiently large simulation size, of order the sizes we are using here. Previous gyrokinetic studies have found that some modest non-local turbulence spreading may occur over distances of a few radial eddy sizes, 26 but in the core region of the large tokamaks we are studying here (and even more so at reactor scales) this should usually be a small effect. It should be acknowledged that the separation of scales assumed for the core plasma in this paper may break down in the edge region of the plasma because gradient scale lengths and eddy sizes may not be very different near the edge, so non-local effects may be important there.
This paper is organized as follows: in the next section we state the fundamental assumptions of the multiscale model and present the closed system of equations that results from gyrokinetic expansion of the Maxwell-Boltzmann system. In Sec. III, we describe the numerical scheme used in TRINITY to simulate the multiscale gyrokinetic system of equations. We also give estimates for the space and time domain savings provided by the multiscale scheme. Sec. IV contains results from TRINITY simulations of L-mode and H-mode discharges from JET and ASDEX Upgrade. We show that the numerical data from TRINITY is in good quantitative agreement with reconstructions of experimental data. Finally, we conclude in Sec. V with a summary and a discussion of possible future directions for research.

II. THEORETICAL FRAMEWORK
In this section, we state the fundamental assumptions present in our multiscale model and present the resulting closed system of equations that must be solved. These equations, which are a rigorous asymptotic limit of the full Maxwell-Boltzmann system, have been derived in detail in Refs. 27 and 28. We include a brief overview of the key results here for completeness.
As the starting point for our analysis, we begin with the coupled system consisting of Maxwell's equations and the driven Fokker-Planck equation: where f s = f s (x, v, t) represents the distribution of particles of species s in position (x) and velocity (v) space, C represents the effect of two-particle Coulomb interactions, and S s is a source term accounting for the external injection of particles, momentum, and energy. This system of equations describes all of the important dynamics in fusion plasmas and is consequently intractable, both analytically and numerically. Since we are interested in studying the interaction of the plasma micro-turbulence with the macroscopic profiles, we simplify the system by adopting a variant of the standard δf -gyrokinetic ordering. 6 This ordering imposes constraints on the relative amplitudes and space-time scales of the macro-and micro-physics. We first decompose the distribution function into macroscopic and microscopic quantities by ensemble averaging: f = f + δf , with the angled brackets denoting an ensemble average. Defining the smallness parameter to be ≡ ρ/L, where ρ is the Larmor radius and L is a macroscopic scale length, we order each of the terms in the Maxwell-Boltzmann equations. The assumptions employed in ordering the terms are as follows: (1) The fluctuations are assumed to be low amplitude compared to macroscopic quantities, such that δf ∼ f . This is in good agreement with core measurements from a number of modern fusion experiments, which find density and temperature fluctuations 1% of the macroscopic densities and temperatures; 29,30,31 (2) The microturbulence is assumed to be spatially anisotropic with macro-scale variations along and micro-scale (i.e. Larmor radius) variations across the equilibrium magnetic field. Experimental measurements of turbulence parallel and perpendicular correlation lengths support this hypothesis; 29,32,33 (3) All frequencies of interest are assumed to be well below the ion cyclotron frequency, and the evolution of the macroscopic profiles is taken to be much slower than the turbulent fluctuations (∂ f /∂t ∼ ∂δf /∂t). Again, this transport ordering is in agreement with experimental evidence; 32 (4) We order ∂δf /∂t ∼ C[δf ], with the collision frequency, ν, defined such that ν ∂ ln δf /∂t. Consequently, δf is allowed characteristic scales in velocity space, δv, of size √ v th δv v th . This collision frequency ordering is satisfied even in very collisionless plasmas such as anticipated in ITER; (5) Macroscopic flows are assumed to be comparable to the ion thermal speed, with microscopic flows (i.e. E × B velocity) taken to be much smaller. For simplicity, we consider in this paper the case in which the Mach number of the macroscopic flow is taken to be small as a subsidiary ordering; (6) The external particle, momentum, and energy sources are assumed to affect the system evolution on the confinement time scale, consistent with experiment.
Expanding the distribution function and fields in and applying the above ordering assumptions to the Maxwell-Boltzmann system results in a hierarchy of equations that is ultimately closed by ensemble and flux surface averaging and taking moments of the evolution equation for the lowest order (macroscopic) distribution function, f 0 = f 0 . One finds that f 0 is a gyroangle-independent, shifted Maxwellian whose evolution is governed by the following transport equations: where R is the major radius, φ is the physical toroidal angle, ψ is the flux label, the overline denotes a flux surface average, and V = dV /dψ, with V being the volume enclosed by the flux surface. The evolved quantities n s , p s , and L are the ensemble-averaged density, pressure, and species-summed toroidal angular momentum, respectively. In the low Mach limit we are considering, the density and pressure are constant on flux surfaces. Terms denoted by S represent external sources, with subscripts indicating the relevant injected quantity. The collisional energy exchange frequency, ν ε su is given in Ref. 34. The terms Γ, π, Q, and H are fluxes and heating generally consisting of classical, neoclassical, and turbulent contributions. They are given by · ∇b contains the magnetic drifts, e is the species charge, ρ is the gyroradius vector, and D/Dt = ∂/∂t + u · ∇, with u = Rω(ψ)φ the equilibrium flow velocity.
The components of the first order distribution function, f 1 , and the fluctuating potentials are obtained by solving the neoclassical and gyrokinetic equations, coupled to the low-frequency Maxwell's equations. These equations are all self-consistently obtained as part of the multiscale gyrokinetic expansion. The neoclassical equation governing f 1 is The gyrokinetic equation determining the evolution of δf 1 is where is the safety factor, . R denotes a gyroaverage at fixed guiding center position, R, and the subscript on the D/Dt operator indicates that it is to be evaluated at R. The low-frequency Maxwell's equations are given by The final elements needed to close this system are equations for the evolution of the magnetic geometry. In particular, the magnetic field can be specified in an axisymmetric system if given the poloidal flux, ψ, and the toroidal flux function, I(ψ). The toroidal flux function is evolved by taking the toroidal component of Faraday's Law and flux surface averaging: with the B · ∂A 0 /∂t term obtainable from the neoclassical equation (Eq. 9). Eqs. 9 and 14 are then coupled to the Grad-Shafranov equation, which uses the updated macroscopic pressure from Eq. 4 to obtain ψ and close the system:

III. NUMERICAL METHOD
We describe in this section the numerical model we have developed for solving the system described by Eqs. 2-15. This model is implemented in the multiscale gyrokinetic transport solver TRINITY. 21 Currently, there are a few additional assumptions employed in TRINITY to simplify the multiscale system presented in the previous section. In what follows, we provide a numerical prescription for solving the full system, pointing out the places where the TRINITY model has been simplified.
A simple sketch of our multiscale numerical model is given in Fig. 1. A direct numerical simulation would require us to have a fine space-time mesh over the full device volume and over at least a confinement time. We exploit scale separation present in the system to drastically reduce the domain over which a fine mesh is required. Our assumption of time scale separation between the turbulence and the equilibrium allows us to fix equilibrium quantities while we evolve the turbulence to saturation. Additionally, it allows us to use steady-state, time-averaged fluxes in our transport equations. Consequently, we need only resolve turbulence time scales for short periods of time, between which we can take large time steps characteristic of the confinement time. As the separation of scales gets wider, the simulation domain savings from this approach grows: The cost of simulating small ρ * devices is no greater than that for moderate ρ * devices. The time domain savings for a device like ITER is a factor of hundreds.
Similarly, spatial scale separation allows us to assume that macroscopic quantities (and their associated gradient scale lengths) are constant across the radial domain in which we simulate turbulent dynamics. As long as the turbulence simulation domain is wide enough in each dimension, the turbulence at the ends of the domain is uncorrelated. Statistically periodic boundary conditions then apply. The result of this local approximation is a flux tube simulation domain for the turbulence (Fig. 1), which can be used to periodically map out a flux surface. Comparisons between local and global gyrokinetic simulations have shown that the local approximation is valid for small ρ * , 35 as it must be for the gyro-Bohm scaling suggested by high confinement experiments 32 to hold. Once again, the spatial domain savings increases with the scale separation. On small ρ * devices, the simulation volume is reduced by approximately a factor of one hundred.

A. Discretization of the transport equations
The transport equations (Eqs. 2-4) are stiff, nonlinear partial differential equations. In order to take the large time steps required by our multiscale scheme, we must treat them implicitly. We allow for a general, single-step time discretization, but we primarily use first-order backwards differences for steady-state systems and secondorder backwards differences for time-dependent systems. An adaptive time step is employed, allowing for accurate time evolution with large time steps. The nonlinear terms are treated implicitly by linearizing them using a standard, multi-iteration Newton's method similar to that given in Ref. 20. For instance, the normalized heat fluxes at the m + 1 time level, which are nonlinear functions of the macroscopic profiles, are expanded as follows: where p denotes the iteration index within each time step, Q ≡ (Q/pv th )(a/ρ) 2 , and all quantities are understood to be evaluated at time level m + 1. The vector y contains the profiles of the fundamental macroscopic, timedependent quantities in the simulation. This consists of the two free flux functions from the Grad-Shafranov equation, ψ and I(ψ), as well as the species density and pressure and the species-summed toroidal angular momentum. We are not currently evolving the magnetic equilibrium in TRINITY, so that ψ and I(ψ) are fixed in time.
Discretizing Eq. 16 in space, we obtaiñ where the subscript denotes the spatial index. In the local approximation, the fluxes depend only on the local values of macroscopic quantities and their gradients. The above expression thus reduces tõ In TRINITY, we make the further simplifying assumption that the fluxes depend more strongly on profile gradients than the local values themselves. We then neglect ∂Q/∂y, giving us the final expressioñ If this assumption is not satisfied, it affects only the rate of convergence of the solution, not its accuracy. In order to numerically calculate ∂Q j /∂y j , we must employ finite differences. This requires us to computẽ Q j at multiple values of y j . This is equivalent to calculating the fluxes for both the nominal profiles and for additional profiles corresponding to each gradient that must be perturbed. The total number of flux tube calculations required during each transport time step, N , is given by where n r is the number of radial grid points in the transport solver and n p is the number of macroscopic profiles being evolved.
Once the transport equations are linearized, it is straightforward to implicitly evolve them. The expense of the implicit evolution is negligible when compared with the cost of the nonlinear turbulence calculations so that there is no need to use an approximate Jacobian. Because inversion of the Jacobian is essentially free, a dense matrix, arising from high order spatial derivatives, is tractable.

B. Schematic
To begin the numerical calculation, the initial state of the plasma must be specified. In particular, enough information must be given to calculate the fluxes and heating from Eqs. 5-8. This requires local information about the magnetic equilibrium, as well as values for the macroscopic density, flow, and temperature and their gradients at each of the flux surfaces comprising the radial grid for the transport solver. TRINITY is currently capable of both an analytic and numerical specification of these quantities, with experimental values taken from the publicly accessible ITER profile database. 36 Once these quantities are obtained at each of the radial grid points, TRINITY calls a solver for the fluxes. For the ion neoclassical heat flux, TRINITY currently uses the simplified analytic model given in Ref. 37. All other neoclassical and classical fluxes, which are typically small compared to turbulent fluxes, 38,39 are neglected. For a more accurate treatment of neoclassical effects, one could interface to a code such as that given in Ref. 40, which solves Eq. 9 directly. This will be the subject of future work.
For the turbulent fluxes, there are interfaces within TRINITY to two widely-used, nonlinear gyrokinetic codes, GS2 22 and GENE. 23,24 Additionally, there are options to use the IFS-PPPL model 1 and other simpler analytic models. Because we are using the local approximation, the flux calculations at each radius are independent for a given transport time step. Consequently, each flux tube calculation can be run in parallel, with the only communication occurring when the fluxes are gathered to advance the transport equations. Since each flux tube calculation parallelizes with high efficiency to thousands of processors, our scheme with tens of flux tubes can easily scale to hundreds of thousands of processors.
Once the steady-state turbulent fluxes are computed, they are time-and flux tube-averaged and passed back to the transport solver. The discretized transport equations are then solved to obtain updated densities, pressures, and the species-summed toroidal angular momentum. Boundary conditions are required to obtain unique solutions. At the outermost radius in the simulation, we fix the values of the thermodynamic profiles, typically to values taken from experiment. As the other boundary condition, we take the product of V with the fluxes to vanish at the magnetic axis. Within each transport time step, the equations are iterated until the relative error upon successive iterations is less than a user-specified tolerance. In practice, we find that two iterations is sufficient to obtain accurate results.
Currently, we only consider static magnetic equilibria in TRINITY. This eliminates several terms in Eqs. 2-4, and foregoes the necessity of solving Eqs. 14 and 15. This approximation is strictly valid in the limit of β = 8πp/B 2 m e /m i , where the magnetic geometry evolves on a resistive time which is much longer than the energy confinement time. However, if one were to evolve the magnetic equilibrium, then the next step would be to use the parallel current obtained from Eq. 9 in Faraday's Law (Eq. 14) to calculate I(ψ). When combined with the updated pressure gradient, this allows us to solve for ψ in Eq. 15. With the updated thermodynamic profiles and magnetic equilibrium, we complete the feedback loop by solving for updated fluxes. This process is repeated as many times as is necessary to evolve the macroscopic profiles beyond the time of interest. For steady-state simulations, approximately ten to fifteen transport time steps are typically required in TRINITY for convergence.

IV. SIMULATION RESULTS
We illustrate the utility of our multiscale model in this section by presenting numerical results from TRINITY simulations and comparing these results with experimental data. We consider an L-mode discharge from JET (shot #19649) and H-mode discharges from JET (shot #42982) and ASDEX Upgrade (shot #13151). For simplicity, we consider time slices taken from approximately steady-state periods of each discharge. Initial thermo- dynamic profiles and external sources used in TRINITY are taken during these steady-state time slices from the TRANSP 41 or ASTRA 42 reconstructions provided in the ITER profile database. Some key experimental parameters for these shots are given in Table I. In all gyrokinetic simulations, we consider electrostatic turbulence with gyrokinetic ions and electrons.
First, we consider JET shot #19649. 43 This was a standard JET L-mode pulse with 9.2 M W of neutral beam heating. For the TRINITY simulation, nonlinear fluxes from GS2 were used in the transport solver. We used a Miller local equilibrium model 44 for the magnetic geometry in these GS2 simulations, with the necessary parameters taken from the ITER profile database. In physical space, we used 16 grid points along the equilibrium magnetic field and a 40 × 25 grid in the perpendicular plane, with perpendicular box widths at the outboard midplane of 64 ρ i . In the dealiased Fourier space, this corresponds to covering |k θ ρ i | = 0, 0.1, 0.2, ...0.8 and |k r ρ i | = 0, 0.1, 0.2, ..., 1.3. The velocity space grid consisted of 10 velocities and 32 pitch angles, giving approximately 10 7 mesh points for each two-species GS2 simulation. These simulations employed a hyperdiffusion operator whose magnitude is scaled by the shearing rate of the turbulence, which provides a sub-grid model of the cascade of fluctuations to smaller scales. 45,46 Previous tests of these kinds of sub-grid models have found that they can reduce the needed resolution, but more detailed convergence studies in the future could be useful. Note also that because of the nonlinearities and stiffness of the transport, changes in the turbulence level of a few tens of percent at fixed temperature gradient have little effect on the final self-consistent temperature profile.
Each flux tube simulation was run for 10 4 time steps, corresponding to average physical simulation times at each radius of approximately 400-1400 L Ti /v th,i . Electron density and ion and electron pressures were evolved at 8 radial locations, giving a total of 32 flux tube simulations per transport time step. The total number of mesh points required for each transport time step was thus approximately 3 × 10 8 . With a total of 15 transport time steps taken, the simulation lasted approximately 4 hours on 5760 CRAY XT4 processors.
A comparison of the steady-state profiles calculated in TRINITY with the TRANSP-reconstructed experimental profiles is given in Fig. 2. We find good agreement for all profiles across the simulation domain (r/a = 0.053 − 0.8), with an RMS relative error amongst all profiles ( (σ 2 n + σ 2 Ti + σ 2 Te )/3) of 12% (Table II). The total and incremental stored energies differ from TRANSP values by 9% and 12%, respectively. A comparison of the gradient scale lengths is given in Fig. 3. For r/a < 0.6, the ion temperature gradient scale lengths from TRINITY and TRANSP match almost perfectly. However, the electron density and temperature gradient scale lengths do not agree as well, despite reasonable profile agreement. This illustrates a difficulty in using experimental profile measurements in standalone gyrokinetic turbulence calculations: small changes in experimental profiles can lead to significant changes in experimental gradient scale lengths, which drastically affect the turbulent fluxes. An example of this is seen in Fig. 4. Here, we compare the volume integrated source terms to the flux surface integrated fluxes, which should be equal in steady state. We see that this balance is satisfied for the self-consistent profiles obtained in TRINITY, but not with the profiles taken from TRANSP. As an example of the sort of diagnostic information available in our multiscale gyrokinetic calculations, we show the steady-state fluctuation amplitudes for density, temperature, and electrostatic potential in Fig. 5. All fluctuation amplitudes are small compared to equilibrium quantities, giving us a check on our δf f assumption. Interestingly, we see peaked electron temperature fluctuations at the magnetic axis, while the electrostatic potential increases with radius.
Next, we consider JET shot #42982, 47,48 which achieved a record fusion energy yield of 22 M J with 21.5 M W of neutral beam heating. The plasma in this shot was a 50-50 deuterium-tritium mixture, operating in ELMy H-mode. Again, GS2 was used to calculate the turbulent fluxes, with a Miller local equilibrium model for the magnetic geometry. The perpendicular spatial grid used was the same as for the L-mode discharge, but the resolution in the remaining dimensions was increased to 24 parallel grid points, 12 velocities, and 40 pitch angles. We again employed a hyperdiffusion operator to prevent a cascade to sub-grid scales. Each flux tube simulation was run 10 4 time steps, with average simulation times ranging from 250-1000 L Ti /v th,i . This time the electron density and ion and electron pressures were evaluated at 10 radial locations, resulting in 40 flux tube simulations per transport time step. A total of 20 transport time steps were taken, and the full simulation ran just under ten hours on 8640 CRAY XT4 processors.
As with the L-mode case, all profiles show relatively good quantitative agreement with their TRANSP counterparts (Fig. 6), with an RMS relative error averaged over all profiles of 12%. In this case, there appears to be a systematic over-prediction of the ion and electron heat fluxes over the outer half of the minor radius, with the profiles pinned by the critical gradient (Fig. 7). A possible explanation for this discrepancy is the fact that the radial flow shear, which has been found to have a significant effect in this discharge, 49 was not included in our GS2 simulations. This capability does exist in GS2 50,51 and will be included in future TRINITY studies with evolving toroidal angular momentum. Again, fluctuation amplitudes are on the order of a percent of equilibrium amplitudes, with electron temperature fluctuations peaking on-axis and electrostatic potential fluctuations increasing with radius (Fig. 8).
In order to study the effect of the edge temperature on the profiles, we repeated this simulation with electron and ion edge temperatures increased by 20%. The results are shown in Figs. 9-11. We see that the 20% increase in edge temperature leads to an increase of approximately 14% at the magnetic axis, as expected from the stiff profiles indicated in Fig. 7. The gradient scale lengths are similar to the base case across most of the minor radius, with the only significant discrepancies occurring near the edge where the stiffness of the profiles is less pronounced (Fig. 10). Fluctuation levels also do not differ signifi- cantly from those obtained for the base case (Fig. 11).
Finally, we consider ASDEX Upgrade shot #13151, 52 which was an ELMy H-mode discharge with 5 M W of neutral beam heating. Here, we used GENE to calculate the turbulent fluxes, with a numerical magnetic equilibrium generated by TRACER. 53 We used 16 parallel grid points and a 64 × 48 grid in the perpendicular spatial plane, with perpendicular box widths of approximately 64 ρ i at the outboard midplane. The velocity space was sampled with 32 parallel velocities and 8 magnetic moments, resulting in a total of approximately 2.5 × 10 7 mesh points for each GENE simulation. A linearized Landau-Boltzmann operator was used to model the effect of collisions. 54 The number of time steps taken in each flux tube simulation varied from about 3-4×10 4 depending on radial location, corresponding to simulation times of approximately 400-1000 L Ti /v th,i . Electron density and ion and electron pressures were again evolved at 8 radial locations, resulting in 32 flux tube calculations per transport time step. Within each transport time step, the total mesh points required was thus ≈ 8 × 10 8 . Due to the increased spatial resolution in the perpendicular domain and the use of a physical collision operator, the GENE simulations took somewhat longer than the earlier GS2 simulations. In total, the simulation took just under 24 hours for 16 transport time steps on 16384 CRAY XT4 processors.
The electron density and temperature agree relatively well across the minor radius, but the ion temperature is under-predicted near the magnetic axis (Fig. 12). This is reflected in the profile gradient scale lengths shown in Fig. 13. The RMS relative error averaged over profiles is nonetheless only 12%. A lack of flow shear in the calculation of the turbulent fluxes is again a possible explanation for the discrepancy in ion temperature near-axis. Another possibility is the fact that no MHD model was used for the q profile computed by ASTRA, which resulted in q < 1 inside r/a ≈ 0.4, with q(0) < 0.5. Consequently, there is significant uncertainty in modeling the system in this region.

V. CONCLUSIONS
We have presented in this paper a complete theoretical and numerical model for the interaction of microand macro-physics in axisymmetric fusion devices. This model arises from a rigorous asymptotic expansion of the full Maxwell-Boltzmann system. The ordering assumptions used in the expansion were given in Sec. II, along with the resulting closed system of equations for the evolution of the macroscopic thermodynamic profiles and magnetic geometry (Eqs. 2-15). These profiles depend on fluxes and heating (Eqs. 5-8) arising from classical, neoclassical, and turbulent dynamics, requiring solution of the neoclassical and gyrokinetic equations (Eqs. 9 and 10). In order to evolve the magnetic geometry and close the system, the toroidal component of Ampere's Law (Eq. 14) and the Grad-Shafranov equation (Eq. 15) also have to be solved.
In Sec. III, we described the numerical scheme used in TRINITY to solve the multiscale gyrokinetic system of equations. The key idea in the approach is to use scale separations in space and time to embed a fine mesh for turbulent dynamics in a coarse grid for the transport solver. A local, nonlinear gyrokinetic code is used to calculate the turbulent fluxes, which are then passed to the transport solver. The macroscopic thermodynamic profiles are then evolved in a manner consistent with the small ρ * limit of the orderings used to derive the nonlinear gyrokinetic equation. Thus, TRINITY can be used to simulate time-dependent experimental phenomena. A Newton method is used to devise an implicit time stepping algorithm, allowing for large time steps characteristic of the confinement time. With updated pressure profiles, one could then couple to solvers for the toroidal component of Faraday's Law and the Grad-Shafranov equation to evolve the magnetic geometry.
Simulation results from TRINITY are provided in Sec. IV, with comparisons to JET and ASDEX Upgrade plasmas. Relatively good agreement is found for all evolved profiles (density and electron/ion pressures), with an RMS relative error averaged over all profiles of 12%. Fluctuation levels for density, temperature, and electrostatic potential are on the order of a percent across the minor radius, in general agreement with experimental evidence.
While currently capable of faithfully simulating a range of interesting experimental conditions, there are several useful additions that could be made to the numerical model implemented in TRINITY. Coupling to a neoclassical code which solves Eq. 9 would provide more accurate neoclassical fluxes, which can be important in certain experimental regimes (for instance, when strong shear flow partially suppresses turbulent flux levels). Additionally, it would allow for the calculation of the parallel current necessary to evolve the safety factor and the poloidal flux. Coupling to a Grad-Shafranov solver would then allow for a study of the effects associated with evolving magnetic geometry. Additionally, coupling to a linear MHD code could be useful in monitoring macroscopic profiles to ensure that MHD stability boundaries are not crossed.
In conclusion, we emphasize that the multiscale gyrokinetic model presented here is not meant to be a comprehensive model for all physics present in a tokamak discharge. Instead, it is meant to be used as a tool for studying the self-consistent interaction between microturbulence and macroscopic profiles. The ultimate goal of this approach is to obtain both a better qualitative and quantitative understanding of this interaction in order to enhance our ability to suppress turbulence and improve confinement in fusion experiments.