Towards optimal explicit time-stepping schemes for the gyrokinetic equations

The nonlinear gyrokinetic equations describe plasma turbulence in laboratory and astrophysical plasmas. To solve these equations, massively parallel codes have been developed and run on present-day supercomputers. This paper describes measures to improve the efficiency of such computations, thereby making them more realistic. Explicit Runge-Kutta schemes are considered to be well suited for time-stepping. Although the numerical algorithms are often highly optimized, performance can still be improved by a suitable choice of the time-stepping scheme, based on spectral analysis of the underlying operator. Here, an operator splitting technique is introduced to combine first-order Runge-Kutta-Chebychev schemes for the collision term with fourth-order schemes for the remaining terms. In the nonlinear regime, based on the observation of eigenvalue shifts due to the (generalized) $E\times B$ advection term, an accurate and robust estimate for the nonlinear timestep is developed. The presented techniques can reduce simulation times by factors of up to three in realistic cases. This substantial speedup encourages the use of similar timestep optimized explicit schemes not only for the gyrokinetic equation, but also for other applications with comparable properties.

💡 Research Summary

The paper addresses the computational bottleneck in large‑scale gyrokinetic simulations, which are essential for predicting turbulence in magnetized laboratory and astrophysical plasmas. Gyrokinetic equations evolve a five‑dimensional distribution function and contain three distinct operators: a linear operator L (parallel streaming, curvature and ∇B drifts, gradient drives), a nonlinear E×B advection operator N, and a linearized collision operator C. Traditional codes such as GENE, GKW, and others use a global explicit Runge‑Kutta (RK) integrator. Because explicit methods are limited by the most restrictive eigenvalue of the discretized operator, the timestep Δt is often forced to be very small, especially when fast Alfvén‑type waves or strong E×B advection dominate.

The authors propose two complementary strategies to enlarge the stable timestep without sacrificing accuracy. First, they apply operator splitting: the stiff collisional term C is integrated with a first‑order Runge‑Kutta‑Chebyshev (RKC) scheme, which possesses an extended real‑axis stability region, allowing explicit treatment even when collision frequencies are high. The remaining terms (L and N) are advanced with a fourth‑order RK scheme (RK4 or its variant RK4M), which provides a wide stability region along the imaginary axis where the high‑frequency wave modes lie. By matching each sub‑operator with the scheme that best fits its spectral characteristics, the overall method retains explicitness while dramatically increasing the allowable Δt.

Second, the paper develops a robust nonlinear timestep estimate based on the observation that the E×B advection term shifts eigenvalues purely along the imaginary axis. By freezing the advection velocity at the current timestep, measuring its maximum magnitude, and applying a CFL‑like condition Δt ≈ C·Δx/|vχ|max, the authors obtain a dynamic, physics‑based bound that is far less conservative than the traditional global bound. They also compute the exact linear eigenvalue spectrum using SLEPc to identify the most restrictive linear modes (fast kinetic shear Alfvén waves, electron streaming, curvature drifts) and combine these limits with the nonlinear estimate.

When implemented in the GENE code, the combined operator‑splitting and adaptive timestep strategy yields speed‑ups of up to a factor of three in realistic test cases, including core and edge tokamak simulations where collisionality varies widely. Accuracy is maintained because the overall error tolerance required for turbulent statistics is modest (relative error ≈10⁻³), making low‑order time integration acceptable. The additional computational cost of extra stages in the RKC scheme is offset by the larger Δt, and memory usage remains unchanged.

The authors argue that the methodology is not limited to gyrokinetics. Any system whose discretized operator exhibits a clear separation between real‑axis damping (e.g., diffusion, collisions) and imaginary‑axis oscillations (e.g., wave propagation) can benefit from a similar split‑scheme approach. Consequently, the work provides a practical pathway to more efficient, explicit time integration for a broad class of high‑dimensional plasma and fluid models on modern massively parallel supercomputers.