Bringing global gyrokinetic turbulence simulations to the transport timescale using a multiscale approach
Jeffrey B. Parker Lynda L. LoDestro Daniel Told Gabriele Merlo Lee F. Ricketson Alejandro Campos Frank Jenko Jeffrey A. F. Hittinger






parker68@llnl.gov
Lawrence Livermore National Laboratory, Livermore, CA 94550, USA

Lawrence Livermore National Laboratory, Livermore, CA 94550, USA

Max-Planck-Institut für Plasmaphysik, Boltzmannstraße 2, 85748 Garching, Germany

University of California, Los Angeles, California 90095, USA

Lawrence Livermore National Laboratory, Livermore, CA 94550, USA

Lawrence Livermore National Laboratory, Livermore, CA 94550, USA


Max-Planck-Institut für Plasmaphysik, Boltzmannstraße 2, 85748 Garching, Germany
University of California, Los Angeles, California 90095, USA

Lawrence Livermore National Laboratory, Livermore, CA 94550, USA




The vast separation dividing the characteristic times of energy confinement and turbulence in the core of toroidal plasmas makes first-principles prediction on long timescales extremely challenging.  Here we report the demonstration of a multiple-timescale method that enables coupling global gyrokinetic simulations with a transport solver to calculate the evolution of the self-consistent temperature profile.  This method, which exhibits resiliency to the intrinsic fluctuations arising in turbulence simulations, holds potential for integrating nonlocal gyrokinetic turbulence simulations into predictive, whole-device models.


Bringing global gyrokinetic turbulence simulations to the transport timescale using a multiscale approach
    Jeffrey A. F. Hittinger
    
=========================================================================================================







§ INTRODUCTION

First-principles simulations of tokamaks will play an increasingly important role in the mission to understand and control magnetized fusion plasmas.  Validated whole-device models aid in interpreting experimental results and help guide the planning of experiments to come.  Predictive modeling is particularly important for reactor-scale burning-plasma machines, with their enormous cost, enormous stored energy, and constrained number of shots.  While reduced modeling will be used for mapping out parameter space and initial shot planning in these devices, high-fidelity simulations will be required before finalizing shots, to ensure as well as possible that each shot will evolve as planned, and in particular that dangerous operating limits will not be approached.

One area of importance is high-fidelity prediction of heat and particle transport, which tend to be dominated by turbulent motions.  In this area, the workhorse framework is the gyrokinetic theory <cit.>.  Gyrokinetic simulations of turbulence for tokamak core plasma have reached a state of maturity <cit.>, enabling a recent increase in emphasis on validating simulations to experiment <cit.>.  These validation efforts have had some success; for example, when experimentally measured quantities such as the temperature and density profiles are provided as simulation inputs, numerical predictions of heat transport are often in quantitative agreement with experimentally inferred levels.  These results suggest that turbulence simulations are nearing a point where quantitative predictive capability is becoming plausible.

For such predictive value to fully materialize, gyrokinetic simulations must be incorporated into modeling on long experimental timescales.  Recently, some full-f gyrokinetic simulations have been performed for an energy confinement time <cit.>.  Because simulating to the energy confinement timescale is computationally expensive, however, gyrokinetic simulations often last for shorter durations and do not take into account the self-consistent evolution of plasma parameters occurring on the longer timescale.  Due to the difference in timescales between turbulence and transport, a direct approach to self-consistent simulation from turbulence to energy-confinement timescales is exceedingly challenging.  An estimate of the required computational resources, assuming gyro-Bohm scaling of the turbulent flux <cit.>, the use of magnetic-field-aligned coordinates, and a fixed radial resolution relative to the gyroradius, yields a computational-cost scaling as ∼^-3, where  is the ratio of the ion thermal gyroradius to the minor radius a of the tokamak.  This cost renders such direct simulations of large reactors infeasible at the present time.  To illustrate, proceeding from transport-timescale gyrokinetic simulations of relatively small tokamaks (≈ 1/100) to ITER (≈ 1/1000)  <cit.> increases the computational cost of an already expensive scheme by three orders of magnitude.

A multiple-timescale approach could provide significantly more advantageous scaling to larger tokamaks such as ITER, by leapfrogging the timescale gap and reducing the computational-cost scaling from ^-3 to ^-1.  The ^-2 difference in cost scaling results from simulating turbulence phenomena and transport phenomena only on their natural timescales and avoiding the simulation of all intermediate timescales.  A multiple-timecale scale approach has previously been demonstrated with local, flux-tube gyrokinetic simulations <cit.>.  However, local simulations miss significant nonlocal effects such as turbulence spreading, transport avalanches, and transport in the core of high-performance discharges or internal transport barriers <cit.>.  A fully predictive capability requires global simulations.  Some early work focused on global gyrokinetic simulations with profile adjustments for matching experimental power balance <cit.>.  In this letter, we focus on global gyrokinetic simulations, presenting a unique method to bridge the timescale separation.



§ THEORETICAL FRAMEWORK

Beginning from the plasma Fokker–Planck equation and assuming a separation of timescales, one can derive coupled equations describing the fast microturbulence and the slowly evolving background Maxwellian.  Sugama and Horton first carried out this procedure, and more recently, Abel et al. have elucidated it in greater detail <cit.>.  We provide a brief overview here.

The formal orderings assume _turb/ ∼^2, 1/_turb∼, /̊L ∼∼, and f/f_0 ∼, where _turb and  are the characteristic timescales of turbulence and transport, respectively,  is the ion cyclotron frequency, $̊ is the ion gyroradius,Lis a characteristic scale length of the background density, temperature, or magnetic field, andf_0andfare the background and fluctuating distribution functions.  The smallness ofis well established in the core of a tokamak.  In the pedestal and edge region, where background profiles can vary on scales comparable to the gyroradius and fluctuations may be large, these orderings may break down.  An asymptotic expansion inyields the gyrokinetic equation and evolution equations for the density, angular momentum, and pressure of each species.  The latter determine the density, rotation velocity, and temperature of the background Maxwellianf_0, which is in local thermal equilibrium on each flux surface.

On the microturbulence timescale,f_0remains approximately constant.  The microturbulence and turbulent fluxes are assumed to be determined by the macroscopic profiles.  While the turbulence rapidly fluctuates, appropriate spatial and temporal averages over these quantities are assumed to evolve only on the slow timescale.





 §.§ Multiple-timescale numerical method

Although gyrokinetic codes for simulating turbulence on the fast timescale are well established, solving the set of coupled gyrokinetic and transport equations is still an open problem.  Even though the gyrokinetic equation is five-dimensional (5D) and the transport equations are 1D, solving the transport equations self-consistently is quite challenging.  The difficulty is that the turbulent flux is strongly nonlinearly dependent on the profiles.  One is forced to confront this nonlinearity when using implicit time-stepping, which is highly desirable to enable long time steps.

We use the LoDestro method <cit.> for coupling global gyrokinetic simulations to a transport equation.  We describe the method, which comprises several elements, as applied to the ion pressure equation:
	
    3/2 V' pt + x V'  = V' S,

wherepis the ion pressure,xis a radial flux coordinate,V' = dV/dxis the differential volume of a flux surface,= Q̌ ·∇xis the averaged radial projection of the heat flux,·denotes a flux surface and intermediate time average,Q̌ = ∫dv̌  12 mv^2 v̌ fis the heat flux, andSrepresents other local terms and sources.  For simplicity in this work, some of the terms in the full equation, such as turbulent heating and evolution of the magnetic equilibrium, have been neglected.  The heat flux includes both turbulent and neoclassical physics.

In the LoDestro method, the first element is to represent the turbulent component of, denoted[p]to emphasize it is a functional of the pressure profilep, as the sum of diffusive and convective contributions.  Introducing a subscriptmto represent the time index and a subscriptlto represent the index of iteration within a timestep, we write
	
    _m,l→ -D_m,l-1 (∂_x p_m,l) + c_m,l-1 p_m,l.

This key step of lagging the effective transport coefficientsD,c(evaluating them at the previous iteratel-1) provides a practical scheme to iterate and converge to the nonlinear solution.  There is some freedom in determining how to split the flux betweenDandc.  For instance, we write
	
    D_m,l-1   = -þ_l-1[p_m,l-1] / ∂_x p_m,l-1, 
    
    		c_m,l-1   = (1-þ_l-1) [p_m,l-1] / p_m, l-1,

where one may choose the flux-split parameterþin a variety of ways.  Convective contributions can be important if the flux is nondiffusive.  Here, we chooseþto be 1 or nearly 1, unless fluxes run up a gradient, in which case we chooseþlocally to be zero.  The iteration scheme becomes
	
    3/2 V' p_m,l - p_m-1/ t + ∂_x [    V' ( -D_m,l-1∂_x p_m,l
        + c_m,l-1 p_m,l) ] = V' S_m,

wheretis the time step, and for notational convenience only the turbulent component of the flux has been written.  The scheme is fully time-dependent, although the first demonstration presented here seeks a steady-state.  Equation iterationequation is linear in the new iteratep_m,land hence straightforward to solve.  Despite the simplicity in representing the turbulent flux, if the iteration inlconverges, then a correct answer is obtained.  This can be seen by substituting fluxrepresentation into fluxsplit, in which case a converged answer means the representation for the flux on the right-hand side is equal to the actual turbulent flux on the left-hand side.

The second element modifies the basic idea in fluxrepresentation in order to stabilize the iteration.  It was recognized in Ref. <cit.> that a diffusion coefficient's depending strongly on the gradient of the profile causes the iteration of Eqs. (<ref>–<ref>) to be numerically unstable.  The iteration can be stabilized by a relaxation of the diffusion coefficient, or similarly, of the turbulent flux.  Equation fluxrepresentation is replaced by
	
    D_m,l-1   = -þ_l-1_m,l-1/ ∂_x p_m,l-1, 
    
    		c_m,l-1   = (1-þ_l-1) _m,l-1/ p_m, l-1,

where
	
    _m,l-1   =  [p_m, l-1] + (1-) _m,l-2, 
    p_m,l-1   =  p_m, l-1 + (1-)  p_m, l-2.

In addition to stabilizing the iteration, a secondary purpose of the relaxation parameteris to provide an effective averaging over previous iterates, reducing the statistical variance from the turbulence simulation.

The third element of the method attempts to minimize the amount of computational time spent in the turbulence simulation, which dominates the computational cost.  In our procedure, we run the turbulence simulation for a few eddy times and calculate an average flux_l-1over that duration (suppressing the time-index subscriptmhere for simplicity).  The relaxed flux_l-1is computed via relaxedflux, which is then used to find a new profile iteratep_lusing Eqs. (<ref>–<ref>) and a newp_lusing relaxedprofile.  Then, the process repeats for the next iteration.  The turbulence simulation runs with the new background profilep_l(background Maxwellianf_0), while using the turbulent stateffrom the end of the previous iteration as the initial condition.

This procedure efficiently takes advantage of thefrepresentation in gyrokinetic simulations and avoids the need to start an entirely new simulation for each iteration.  The time-saving action of restarting a simulation with a newf_0and the oldfcan result in transient bursts associated with the fast relaxation to a self-consistentf, but if the variation inf_0per iteration is small enough, the bursts have minor impact.  Computational efficiency is also promoted by averaging the flux over just a few turbulence eddy times per iteration.  This strategy may result in the turbulent flux lagging a couple iterations behind the profile iterates, but should not cause significant delay in convergence.

We now compare the LoDestro method with a Newton-based method, which was used by earlier works with local gyrokinetic simulations <cit.>.  Those works developed an iteration scheme by expanding the flux as_l = [p_l] ≈[p_l-1] + / p |_p_l-1 ·(p_l - p_l-1).  Such a Newton-based approach, when applied to global simulations, becomes much more challenging because the flux depends in principle on the profile over the entire radial domain, rather than just its local value and gradient as in the local limit.  The Jacobian/ pis therefore significantly higher dimensional; the  number of nonzero elements is proportional toN_r^2in the global problem rather than toN_ras in the local problem, whereN_ris the number of radial grid points.  Moreover, each Jacobian-vector product is an expensive prospect because it requires an additional full turbulence simulation to implement a finite-difference estimate of the derivative.  In addition to being costly, computation of a derivative by finite difference is challenging as it amplifies error inarising from turbulent fluctuations.

The LoDestro method, which employs a variant of fixed-point, or Picard, iteration, avoids these pitfalls.  Not needing to project differences offers a large gain in efficiency, as the lower sensitivity to error allows the use of as little as a couple eddy times of gyrokinetic simulation per iteration.  This gain compensates for the primary disadvantage of fixed-point iteration, that an unaccelerated scheme theoretically achieves linear rather than superlinear asymptotic convergence. 
	


 §.§ Simulation Setup

A newly developed 1D transport solver, , implements the LoDestro method.  , a Python code, is open source <cit.>.  Here,  is coupled to the global gyrokinetic code GENE <cit.> to simulate ion-temperature-gradient-driven turbulence.  

The GENE simulations are collisionless, electrostatic, and use deuterium ions and adiabatic electrons.  Circular magnetic geometry is used <cit.>, with safety factorq = 0.868 + 2.2 ^̊2, where=̊ r/aand the radiusrplays the role of the flux coordinate.  The toroidal rotation velocity is set to zero.  The density profile is prescribed and fixed in time.  The electron temperature is held equal to the ion temperature.  GENE is run in the “gradient-driven” mode, in which a Krook-type source prevents the average profile from varying too much within an iteration as the gyrokinetic equation is integrated in time <cit.>.

 solves the ion pressure equation in the domain0 ≤≤̊0.9, and GENE solves the gyrokinetic equation in the domain0.1 ≤≤̊0.9.   spatially discretizes with second-order finite differences and uses Neumann boundary conditions at the magnetic axis and Dirichlet conditions (fixed temperature) at the outer radius=̊ 0.9.

Two simulations are described here, with key parameters given in Table <ref>.  In each simulation,  seeks a steady-state (t = ∞) solution to Eq. transportequation.  The tokamak in Simulation 1 is roughly half the size (≈1/149) as that in Simulation 2 (≈1/292).  In Simulation 1,120 ×16 ×16 ×78 ×48grid points are used for the radial, binormal, coordinate along the magnetic field, parallel velocity, and magnetic moment dimensions, respectively, in GENE, and 27 radial grid points are used in .  In Simulation 2,240 ×16 ×16 ×78 ×48grid points are used in GENE and 53 radial grid points are used in .  In lieu of a full neoclassical model, a small uniform diffusivity is included in the heat flux.  Further details on the simulations can be found in the Supplemental Material.



The setup for 's simulation is as follows.  In each iteration, GENE runs for 50R_0 / v_T_i, whereR_0is the major radius andv_T_iis the ion thermal velocity.  The simulation is run for 50 iterations.  The relaxation parameter is= 0.3.  Prior to commencing iterations coupled with , GENE is run standalone for 600R_0 / v_T_ito seed the initial condition.





§ RESULTS

The results from Simulation 1 are shown in Figure <ref>, which depicts the applied heat source, the temperature profile, the normalized temperature gradient, and the total heat flow (turbulent+added diffusivity)H = V'.  The initial condition (dashed curve) and a few early iterates (grey curves) are plotted for both temperature and heat flow.  After roughly 10–20 iterations, a noise floor is reached, where fluctuations inprevent any further convergence.  More iterations do not substantially help because the relaxation in Eqs. (<ref>–<ref>) has a finite memory.

Focusing on iterates 31–50, we take the mean as the solution (red curve) and the standard deviation (light blue shaded region) as a rough measure of uncertainty.  The maximum coefficient of variation across radius is 7.1% for the temperature.  For the heat flow, the maximum is 11.3% for≥̊0.28; greater relative fluctuations occur at smaller$̊ because the heat flow is near zero.  There is a small bump in temperature at ≈̊0.22, at which the heat flux flows up the temperature gradient.  The bump is likely due to the strong localization of the heat source in this region.  The heat flux is still outward, an effect probably due to turbulence spreading from the adjacent region.
	
In steady state, the heat flow H can be determined exactly by integrating Eq. (<ref>) over radius; H(r) must equal ∫_0^r dr  V'(r) S(r).  The average of the numerically simulated heat flow over the final 20 iterations agrees well with this exact result (black dots).  The heat flows in both Simulation 1 and 2 are insensitive to the value of the ad-hoc diffusivity, as the heat flux resulting from this diffusivity is miniscule compared to the turbulent flux.
		
Figure <ref> presents similar results for Simulation 2.  The multiscale coupling has also found the solution here, confirmed by the consistency of the numerical heat flow with the integral of the source, as shown in Figure <ref>(d).  The maximum coefficient of variation across radius is 5.7% for the temperature and 8% for the heat flow (≥̊0.35).  It is interesting to note that the self-consistent heat flow is quite different from the initial heat flow, even though the final temperature profile does not appear that different from the initial temperature profile.  This behavior reflects the well-known nature of stiff transport, where small profile variations can lead to large changes in the flux.

Figure <ref> shows results of simulations, using the parameters of Simulation 1, but starting from three different initial temperature profiles.  The steady-state temperature profile, normalized temperature gradient, and heat flow obtained in the three cases are nearly identical.

As a test of sensitivity to numerical parameters, we varied the flux-split parameter þ and redid Simulation 1 using constant þ = 0.5 or þ=0.8.  After iterating to convergence, we obtained the same steady-state solution as before.

	

	

	




§ DISCUSSION

The energy confinement time can be estimated from _E ≈ W_i / P_in, where W_i is the stored energy in the ions in the steady-state solution.  We find _E ≈ 26 ms for Simulation 1 and _E ≈ 55 ms for Simulation 2.  The ratio of confinement times does not follow gyro-Bohm scaling / _turb∼^-2 exactly, as is known for small tokamaks <cit.>.


The total time used for gyrokinetic simulation is a useful measure of computational cost and efficiency of the multiscale method.  In both cases, this expenditure is 3100 R_0 / v_T_i, which corresponds to a physical time of 23 ms for both Simulations 1 and 2.  Hence, our calculation has successfully bridged the timescale separation.  Additionally, opportunity exists to optimize performance through parameter tuning.

Fluctuations in the turbulent flux create a tradeoff between computational expediency and precision.  An example of the fluctuations can be seen in Fig. <ref>(a) (inset); the largest bursts can be three times as large as the mean.  These fluctuations, inherent to the equations regardless of the method of simulation, set a noise floor and limit the achievable precision.  The variance of the fluctuations can be reduced through greater averaging of the flux, which in practice is accomplished by increasing the simulation time per iteration or decreasing the relaxation parameter , either of which increases the computational cost.  It is worth noting that uncertainties in experimental measurements and limitations in physics models may set bounds on the precision desirable in practice, as additional precision requires greater computational resources while providing little added value beyond a certain point.


In summary, we have demonstrated multiple-timescale coupling of global gyrokinetic simulations to a transport solver.  Using the LoDestro method, which prescribes an algorithm for solving the nonlinear implicit timestep of the transport equation, we have solved for the steady-state temperature profile.  This method holds potential to play an important role in the integration of nonlocal gyrokinetic turbulence simulations into predictive, whole-device models.  Extensions of this work that allow for evolution of ion pressure, electron pressure, and density simultaneously are underway.



Useful discussions with Ian Abel are acknowledged.  This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract No. DE-AC52-07NA27344, in part through the SciDAC Partnership for Multiscale Gyrokinetic Turbulence.  This research was supported by the Exascale Computing Project (17-SC-20-SC), a collaborative effort of the U.S. Department of Energy Office of Science and the National Nuclear Security Administration.  This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.


apsrev4-1
 


