GORILLA: Guiding-center ORbit Integration with Local Linearization Approach

1 Institut für Theoretische Physik Computational Physics, Technische Universität Graz, Petersgasse 16, 8010 Graz, Austria 2 Institute of Plasma Physics, National Science Center, “Kharkov Institute of Physics and Technology,” Akademicheskaya str. 1, 61108 Kharkov, Ukraine 3 Department of Applied Physics and Plasma Physics, V. N. Karazin Kharkov National University, Svobody sq. 4, 61022 Kharkov, Ukraine ¶ Corresponding author DOI: 10.21105/joss.04890


Introduction
Extremely hot plasmas with a temperature of the order of hundred million degrees Celsius are needed to produce energy from nuclear fusion. Under these conditions, hydrogen isotopes are fused, and energy is released. The energy release from 1 kg of fusion fuel corresponds approximately to that of 10000 tons of coal. This high energy density has made fusion the subject of worldwide research projects. The confinement of such hot plasmas, however, poses major physical and technological problems for researchers. In particular, complex numerical methods are necessary to understand the physics of such plasmas in complicated toroidal magnetic fields.
An important kinetic approach for simulating the collective behavior of a plasma uses direct modeling of particle orbits. A well-known approximation for computing the motion of electrically charged particles in slowly varying electromagnetic fields is to reduce the dynamical equations by separating the relatively fast circular motion around a point called the guiding-center, and primarily treat the relatively slow drift motion of this point. This drift motion is described by the guiding-center equations in non-canonical Hamiltonian form; see, e.g., (Boozer, 1980), (Littlejohn, 1983) and (Cary & Brizard, 2009).
Here, we provide GORILLA: an efficient code for the purpose of solving the guiding-center equations. This code is a numerical implementation of the novel, quasi-geometric integration method in 3D geometry described by M. Eder et al. (2020). A related method for axisymmetric geometry has been realized in the code K2D (Kasilov et al., 2016). Apart from this quasigeometric approach there exists a variety of codes for tracing guiding-center drift motion including ORBIT (White & Chance, 1984), DCOM (Wakasa et al., 2001), MOCA (Tribaldos, 2001), FORTEC-3D (Satake et al., 2005), VENUS (Isaev et al., 2006), NEO-MC (Allmaier et al., 2008), XGC (Ku et al., 2009), ANTS (Drevlak, 2009), ASCOT (Kurki-Suonio et al., 2009), VENUS-LEVIS (Pfefferlé et al., 2014), SIMPLE  etc; see also (Beidler et al., 2011). Most of these codes solve the guiding-center equations primarily in the plasma core of toroidal fusion devices. Due to its formulation in general curvilinear coordinates, GORILLA is not limited by the field topology. That means that the computation domain of GORILLA covers both the closed field line region (i.e. the plasma core) and the open field line region (i.e. the scrape-off layer).

Summary
GORILLA is a Fortran code that computes guiding-center orbits for charged particles of given mass, charge and energy in toroidal fusion devices with three-dimensional field geometry. Conventional methods for integrating the guiding-center equations use high order interpolation of the electromagnetic field in space. In GORILLA, a linear interpolation of a specific representation of quantities in the guiding-center system employing a spatial mesh is used for the discretization of the electromagnetic field. This leads to locally linear equations of motion with piecewise constant coefficients. As shown by M. Eder et al. (2020), this local linearization approach retains the Hamiltonian structure of the guiding-center equations. For practical purposes this means that the total energy, the magnetic moment and the phase-space volume are conserved. Furthermore, the approach reduces computational effort and sensitivity to noise in the electromagnetic field. In GORILLA, guiding-center orbits are computed without taking into account collisions between particles. Two examples of guiding-center orbits obtained with GORILLA can be seen in Figure 1 where the magnetic field of a real-world fusion device is used, specifically the tokamak "ASDEX Upgrade".

Statement of need
GORILLA is designed to be used by researchers in scientific plasma physics simulations in the field of magnetic confinement fusion. In such complex simulations, a simple interface for the efficient integration of the guiding-center equations is needed. Specifically, the initial position in five-dimensional phase-space is provided in terms of guiding-center coordinates (i.e. guiding-center position, parallel and perpendicular velocity) and the main interest is to retrieve the phase-space position after a prescribed time step. Such a pure "orbit time step routine" acting as an interface with a plasma physics simulation is provided. The integration process itself, however, can also be of great interest. Therefore, a program allowing the detailed analysis of guiding-center orbits, the time evolution of their respective invariants of motion and Poincaré plots is also available. A unique feature of GORILLA is the direct way in which fluxes through cell boundaries can be directly obtained while conserving physical invariants. These properties result from the geometric approach of linear approximation in tetrahedal cells. Thus the computation of moments of the distribution function such as current densities from Monte-Carlo sampling can be much more efficient and stable than in usual orbit tracers. Furthermore, the mesh is agnostic to whether an orbit is inside or outside the magnetic separatrix, and does not rely on flux coordinates. Therefore GORILLA is a flexible tool that is useful in a variety of applications for non-axisymmetric devices, in particular stellarators and tokamaks with 3D perturbations.
GORILLA has already been used by M. Eder et al. (2020) for the application of collisionless guiding-center orbits in an axisymmetric tokamak and a realistic three-dimensional stellarator configuration. There, the code demonstrated stable long-term orbit dynamics conserving invariants. In the same publication, GORILLA was further applied to the Monte Carlo evaluation of transport coefficients. There, the computational efficiency of GORILLA was shown to be an order of magnitude higher than with a standard fourth order Runge-Kutta integrator. Currently, GORILLA is part of the "EUROfusion Theory, Simulation, Validation and Verification Task for Impurity Sources, Transport, and Screening", where it is tested for the kinetic modelling of the impurity ion component. First results of this research project and, in addition, GORILLA's application to the computation of fusion alpha particle losses in a realistic stellarator configuration have been reported by M. . The source code for GORILLA has been archived on Zenodo with the linked DOI: (Michael Eder et al., 2021)