Gyrokinetic Simulation

A focus of the last grant period was the completion of the core GYRO [Candy 03a] solver. GYRO went through various developmental incarnations before converging on the current algorithm set, which has remained unchanged since early 2003. While we had very good results with a fully-explicit solver in the case of adiabatic electron dynamics (as reported during the last grant period), electromagnetic operation became robust only after a switch to an implicit-explicit Runge-Kutta (RK) scheme [Pareschi 02]. The latter class of RK schemes have recently generated significant attention within the applied mathematics community [Kennedy 03]. Now, GYRO operates with all of the physics originally planned: trapped and passing electrons with pitch angle collisions, finite- fluctuations, real flux-surface shape, linear and nonlinear ExB rotation, parallel flow shear, and fully gyrokinetic impurity dynamics. It operates at finite using a WKB-like formulation to treat the small corrections due to radial variation of , , , etc, not present in the flux-tube approximation. Simulations can be carried out at timesteps which are limited by the nonlinear physics and not the electron Courant time. Over the past grant period, we have made substantial progress on a number of fundamental gyrokinetic fronts. We began with a systematic study of the mechanisms which lead to the emergence of sub-gyroBohm scaling (or, the "breaking" of gyroBohm scaling) [Waltz 02] within the context of the adiabatic electron model - the model for which the lion's share of worldwide gyrokinetic simulations have been performed. This landmark study examined the effects of temperature gradient scale length on -scaling, discussed the dynamics and prevention of artificial profile relaxation, and noted the production of transport in linearly-stable regions of the plasma (nonlocal transport). It is the earliest and most exhaustive work on gyrokinetic transport scaling.

Figure 4: Effect of profile shear, S, on transport scaling, as
riginally published in [Waltz 02].

Figure 5: Transition from Bohm-scaled transport at the center to gyroBohm transport at the right simulation
boundary. These results are taken from [Waltz 02].

Next, using more aggressive and numerically expensive electromagnetic simulations, we have been able to reproduce, within experimental uncertainty, the transport observed in DIII-D L-mode plasmas [Candy 03b] as summarized in Fig. 6. In particular, we noted the strong stabilizing effect of electron collisions as well as equilibrium sheared ExB rotation, and were able to numerically recover the Bohm-like scaling separating discharges at and .

Figure 6: From [Candy 03b]. Frame (a) shows the impact of selected physical effects on total energy transport (at ). A dramatic stabilization due to collisions and equilibrium sheared ExB rotation is apparent. On this plot, equal transport values at and 400 indicate gyroBohm scaling, while a ratio of 1:1.6 between 250 and 400 indicate Bohm scaling. Frame (b) shows a rerun of the full-physics case in (a) after reducing the ion temperature gradient by 10%. There is remarkable agreement over the entire radially-simulated domain.

Next, we attempted to resolve a long-standing theoretical controversy in connection with the role of an point in transport barrier formation. Global simulations show that transport due to ITG microturbulence is smooth across an point [Candy 04b], by virtue of the appearance of nonresonant ITG modes in the weak-shear region. The results are supplemented by a detailed multiple-scale perturbative analysis of the ITG-Schrödinger equation in the limit.

Figure 7: Results from [Candy 04b] showing the appearance of a gap mode
(rightmost panel) in the absence of rational surfaces.

We also resolved a second controversy in connection with the local limit of global simulations. Theoretical concerns about the validity of the flux-tube approximation motivated us to undertake a global GYRO -scan of a popular test case over the range [Candy 04a]. We confirmed that the local limit is properly recovered - to with about 5% accuracy - at selected interior radii of the global radial domain. Further, the comparison showed close intercode agreement between GS2 (Dorland), PG3EQ (Dimits) and GYRO.

Figure 8: Comparison of local (GS2, PG3EQ) and global (GYRO, GTC) results for the Cyclone base case. This figure is derived from results published in [Candy 04a].

GYRO is widely ported, runs on local GA and PPPL clusters, the IBM-SP Power3 at NERSC, IBM-SP Power4 at ORNL, and has recently been ported to the ORNL Cray-X1. Performance on the latter cutting-edge platform has been noteworthy [Dunigan 03].

Our primary goal, in the most general sense, is to focus on the application of GYRO to transport simulations for both existing and future devices. Having recently completed studies of L-mode DIII-D discharges [Candy 03b] where we have reproduced the Bohm-like scaling of and , we will attempt to explain the transition to gyroBohm seen in H-mode discharges. Preliminary work shows that H-mode discharges are very close to threshold (that is, to zero ion transport), highlighting the critical need for full physics (collisional, electromagnetic, real geometry) simulations. In parallel, we are preparing for global ITER transport simulations including the dynamics of kinetic electrons and realistic shaping. We are refining the fixed-flux operational mode (in particular, the feedback loop) of GYRO for use in these transport studies. Very-long-time simulations (up to 10 times longer than at fixed-gradient) are needed for fixed-flux operation in order to slowly correct the ion and electron driving gradients to match experimental power flows. GYRO is uniquely suited to these studies, since it (1) uses an Eulerian scheme free from statistical noise and the growing weight problem inherent to the PIC algorithm, and (2) can dynamically change gradients without the need for a code restart. A secondary goal is to continue a broad range of more theoretical gyrokinetic studies, including but not limited to: gyrokinetic impurity transport, transport at high-, effects of weak and reversed shear, supression of turbulence by equilibrium sheared ExB flows, spectral analysis of ITG-to-zonal-flow energy transfer in saturated ITG turbulence, transition from core to edge turbulence. GYRO calculates not only the ion thermal transport, but in fact both the particle and energy flows from all species, as well as momentum and energy exchange. All these are required for experimentally-relevant comparisons. We have also recently incorporated neoclassical driving terms to GYRO and a conservative ion-ion collision operator in the hope of exploring the self-consistent coupling between neoclassical and turbulent transport. There is also a pressing need to support and guide a growing community of users who have already begun to apply GYRO to studies of experimental results from existing tokamaks. Some of these users are also contributing tools and advice which facilitiate simulations based on experimental profile data. These users further motivate code enhancements, bug fixes and new minor version releases. Finally, straightforward modifications to the core solver are planned. These include the addition of compressional magnetic perturbations into the semi-implicit advance, for which no new numerical algorithms are required. Also, in cooperation with Mark Fahey from ORNL, optimizations are planned for the current run-time bottlenecks, including the collision step and the sparse field solve. The former issue may require algorithmic modification in addition to general scalar and vector code optimizations. For the latter, we envision adding an extra layer of parallelism to the field-solve, as it (unlike the advance of the kinetic equation) is only partially distributed. This would be active only for very large processor grids. GYRO is a premiere X1-application code [Dunigan 03] and ORNL enthusiastically endorses such code optimization.