Toroidal Code Research Articles

Implementation and performance evaluation of a communication-avoiding GMRES method for stencil-based code on GPU cluster

In this study, a communication-avoiding generalized minimum residual method (CA-GMRES) is implemented on a hybrid CPU–GPU cluster, targeted for the performance acceleration of iterative linear system solver in the gyrokinetic toroidal five-dimensional Eulerian code (GT5D). In the GT5D, its sparse matrix-vector multiplication operation (SpMV) is performed as a 17-point stencil-based computation. The specialized part for the GT5D is only in the SpMV, and the other parts are usable also for other application program codes. In addition to the CA-GMRES, we implement and evaluate a modified variant of CA-GMRES (M-CA-GMRES) proposed in the previous study Idomura et al. (in: Proceedings of the 8th workshop on latest advances in scalable algorithms for large-scale systems (ScalA ’17), 2017. https://doi.org/10.1145/3148226.3148234) to reduce the amount of floating-point calculations. This study demonstrates that beneficial features of the CA-GMRES are in its minimum number of collective communications and its highly efficient calculations based on dense matrix–matrix operations. The performance evaluation is conducted on the Reedbush-L GPU cluster, which contains four NVIDIA Tesla P100 (Pascal GP100) GPUs per compute node. The evaluation results show that the M-CA-GMRES or CA-GMRES for the GT5D is advantageous over the GMRES or the generalized conjugate residual method (GCR) on GPU clusters, especially when the problem size (vector length) is large so that the cost of the SpMV is less dominant. The M-CA-GMRES is 1.09 ×, 1.22 × and 1.50 × faster than the CA-GMRES, GCR and GMRES, respectively, when 64 GPUs are used.

Verification of gyrokinetic particle simulation of current-driven instability in fusion plasmas. IV. Drift-tearing mode

The drift-tearing instability due to diamagnetic drift effects is verified using the Gyrokinetic Toroidal Code (GTC). First, the classical (2,1) resistive tearing mode is verified in a cylindrical geometry with a fluid model. The dependence of the growth rate of the resistive tearing mode on the beta value of the plasma is obtained and is found to qualitatively agree with the theoretical prediction. A drift-tearing mode is subsequently generated when the equilibrium pressure gradient is significant. In this mode, diamagnetic drift effects result in a reduced growth rate and a real frequency equal to the electron diamagnetic frequency. The scaling relation between the diamagnetic frequency and the growth rate of the drift-tearing mode has been calculated. This relation shows good agreement with the theoretical prediction for a relatively small resistivity; however, an obvious deviation arises when the resistivity is large.

Global gyrokinetic simulation of turbulence driven by kinetic ballooning mode

Nonlinear evolution of the kinetic ballooning mode (KBM) is investigated by extending the global toroidal gyrokinetic simulation code (GKNET) to an electromagnetic regime. It is found that the saturation process of KBM, which is unstable at high normalized pressure β, is significantly different from the ion temperature gradient (ITG) mode, which is unstable at low β. The KBMs get saturated by producing zonal flows and zonal magnetic fields. The production of zonal flow is weak in the initial saturation phase of KBM, which is in contrast to the ITG mode which produces strong zonal flows in the initial saturation phase. However, strong zonal flows are produced in the subsequent evolution of KBM, and a quasisteady state of KBM turbulence is established. In addition to the zonal flows, some low toroidal number modes, which are linearly stable against the KBM, dominate the KBM turbulence. The strong zonal magnetic field is also produced by the KBM. These zonal modes regulate the KBM turbulence.

Kinetic particle simulations in a global toroidal geometry

The gyrokinetic toroidal code has been upgraded for global simulations by coupling the core and scrape-off layer regions across the separatrix with field-aligned particle-grid interpolations. A fully kinetic particle pusher for high frequency waves (ion cyclotron frequency and beyond) and a guiding center pusher for low frequency waves have been implemented using cylindrical coordinates in a global toroidal geometry. The two integrators correctly capture the particle orbits and agree well with each other, conserving energy and canonical angular momentum. As a verification and application of this new capability, ion guiding center simulations have been carried out to study ion orbit losses at the edge of the DIII-D tokamak for single null magnetic separatrix discharges. The ion loss conditions are examined as a function of the pitch angle for cases without and with a radial electric field. The simulations show good agreement with past theoretical results and with the experimentally observed feature in which high energy ions flow out along the ion drift orbits and then hit the divertor plates. A measure of the ion direct orbit loss fraction shows that the loss fraction increases with the ion energy for DIII-D in the initial velocity space. Finally, as a further verification of the capability of the new code, self-consistent simulations of zonal flows in the core region of the DIII-D tokamak were carried out. All DIII-D simulations were performed in the absence of turbulence.

Role of thermal conduction on resistive tearing mode in Tokamaks

We investigate the role of thermal conductivity on the evolution and saturation of the m/n = 2/1 resistive tearing mode (m and n are, respectively, the poloidal and toroidal Fourier mode numbers), by using a 3D toroidal MHD code (CLT). It is found that thermal conductivity has great influence on pressure and current profiles, and further affects the dynamic evolution of the tearing mode. Our simulation results indicate that the linear growth rate of the tearing mode increases, but the saturation level of magnetic islands decreases with increase of thermal conductivity. With a small thermal conductivity, a flattened distribution of the thermal pressure inside magnetic islands leads to a large pressure gradient at the edge that drives a ballooning mode to be unstable. We further found that the radial electric field at the magnetic island boundary and the vortex-like flow inside the magnetic island lead to a poloidally asymmetric transport barrier and reduce the thermal conductivity at the magnetic island boundary, which contributes to formation of an internal transport barrier.

Core-crash sawtooth associated with m/n = 2/1 double tearing mode in Tokamak

The core-crash sawtooth caused by m/n = 2/1 double tearing mode (DTM) in Tokamak is investigated using a three-dimensional toroidal MHD code (CLT). It is found that mode–mode coupling effect plays a crucial role during the core-crash sawtooth. The n = 0 component and n > 1 components generated by mode–mode coupling of the n = 1 modes take up more than 40% of the total kinetic energy and become unneglectable during the core-crash sawtooth. The n = 0 component is crucial for the saturation of DTM. Without the n = 0 component, the DTM will not saturate. The high n components are responsible for the explosive growth of the DTM that finally results in the fast crash of the core plasma temperature.

Loss of bootstrap current in vicinity of magnetic islands

Profiles of the ion density and bootstrap current in the vicinity of magnetic islands are investigated based on the first principles gyro-kinetic particle simulation via the gyro-kinetic toroidal code. The physics on the recovery of the ion density gradient inside the islands in various collision regimes is discussed. Simulation results show that for small magnetic islands, the ion density gradient can survive inside the island due to the combination effect of both the banana-orbit width of trapped ions and the drift-orbit displacement of passing ions. It is suggested that the recovery of the pressure gradient inside small islands may play a more important role in the reduction of driving force of the ion bootstrap current in the evolution of the neoclassical tearing mode, rather than the so-called finite banana-orbit effects.

Global gyrokinetic simulation of microturbulence with kinetic electrons in the presence of magnetic island in tokamak

An electrostatic model has been formulated and implemented in the gyrokinetic toroidal code to study the nonlinear ion temperature gradient (ITG) turbulence in the presence of an n = 1, m = 2 magnetic island. The ions are described by the gyrokinetic equation while the electrons are treated with the drift-kinetic equation. In our simulation, an n = 1, m = 2 electrostatic mode is formed with the same vortex structure of the magnetic island. When the magnetic island flattening effect is turned on, the island vortex mode is well preserved and couples to the n = 0, m = 0 geodesic acoustic mode. Simulation shows that the magnetic island can suppress the ITG turbulence at the island O-point and strengthen it near the X-point. We show that the vortex mode can generate a substantial helical shear flow around the island. We also find that the turbulence and transport are suppressed inside the island and enhanced at the island X-point.

Dynamic evolution of resistive kink mode with electron diamagnetic drift in tokamaks

The dynamic evolution of the m/n = 1/1 resistive kink mode with electron diamagnetic drift is investigated by employing a three-dimensional toroidal Hall-MHD code CLT. It is found that the dependence of the linear growth rate of the resistive kink mode on the electron diamagnetic drift velocity is associated with thermal conductivity. For a large thermal conductivity, the linear growth rate of the resistive kink mode monotonically decreases with increasing electron diamagnetic drift velocity. But, for a small thermal conductivity, the linear growth rate increases first and then decreases with increasing electron diamagnetic drift velocity, which is different from previous studies. The saturation level of the resistive kink mode decreases with increasing electron diamagnetic drift velocity.

Stabilization of tearing modes by modulated electron cyclotron current drive

The influence of modulated-ECCD on m/n=2/1 resistive tearing mode is investigated by a three-dimensional toroidal and non-reduced MHD code CLT. It is found that, after applying a modulated-ECCD, tearing mode instabilities are suppressed and magnetic islands are gradually reduced to a low level, then the width of the magnetic islands exhibit periodic oscillation with the time scale of ECCD modulation frequency. The minimum width of magnetic islands decreases with the decrease of ECCD modulation frequency and increases with the increase of the buildup time of the driven current.

Gyrokinetic simulations of toroidal Alfvén eigenmodes excited by energetic ions and external antennas on the Joint European Torus

The gyrokinetic toroidal code (GTC) has been used to study toroidal Alfvén eigenmodes (TAEs) in high-performance plasmas. Experiments performed at the Joint European Torus (JET), where TAEs were driven by energetic particles arising from neutral beams, ion cyclotron resonant heating, and resonantly excited by dedicated external antennas, have been simulated. Modes driven by populations of energetic particles are observed, matching the TAE frequency seen with magnetic probes in JET experiments. A synthetic antenna, composed of one toroidal and two neighboring poloidal harmonics has been used to probe the modes’ damping rates and quantify mechanisms for this damping in GTC simulations. This method was also applied to frequency and damping rate measurements of stable TAEs made by the Alfvén eigenmode active diagnostic in these discharges.

Destabilization of resistive plasma resistive wall mode by anisotropic thermal transport

Both an analytic model and the toroidal MARS-F code [Liu et al., Phys. Plasmas 7, 3681 (2000)] have been employed to study the effect of the anisotropic thermal transport on the stability of the resistive wall mode (RWM) in resistive tokamak plasmas. The results show that thermal transport destabilizes the resistive plasma RWM, by effectively eliminating the Glasser-Greene-Johnson favorable average curvature stabilization associated with the resistive layer [Glasser et al., Phys. Fluids 7, 875 (1975)]. Modification of the mode eigenfunction is also observed in MARS-F computations, due to the presence of the anisotropic thermal transport.

Gyrokinetic particle simulation of electrostatic microturbulence with impurity ions

Impurity is an important factor that can affect significantly turbulent transport in tokamaks. In order to study the impurity physics, we implement a new impurity module in the gyrokinetic particle simulation code GTC (Gyrokinetic Toroidal Code). With an improved numerical scheme, we expand the validity of gyrokinetic Poisson equation in the GTC to the short wavelength region, for both non-zonal and zonal parts of the perturbed Poisson equation. Verifications of this new scheme are carried out on the linear instability and zonal flow response. The linear simulation of the ion temperature gradient (ITG) instability including the impurity ions shows that the new Poisson solver can obtain the correct linear growth rate and frequency at the thermal ion gyro-radius scale. The residual zonal flow with impurities obtained via the new zonal flow solver is consistent with the numerical and analytical predictions in the large aspect-ratio limit. The nonlinear simulation of the ITG turbulence shows that the turbulent transport is significantly reduced by the impurity ions through decreasing the linear growth rate of the instability.

Development of Finite Element Field Solver in Gyrokinetic Toroidal Code

Development of Finite Element Field Solver in Gyrokinetic Toroidal Code

Multiple branches of resistive wall mode instability in a resistive plasma

The resistive wall mode (RWM) instability is numerically investigated for a toroidal resistive plasma, with results compared to those of an analytic cylindrical model. The full toroidal code MARS-F [Liu et al., Phys. Plasmas 7, 3681 (2000)] is applied for a computational study. The results indicate that there are two branches of unstable RWMs, when the toroidal favorable average curvature effect (the GGJ effect) is taken into account in the resistive layer. In addition, the GGJ physics not only directly affects the mode growth rate, but also indirectly modifies the mode stability by changing the continuum damping through modifying the mode frequency in the plasma frame. Furthermore, the plasma resistivity can either stabilize or destabilize the RWM, depending on the regime of key plasma parameters (e.g., the plasma rotation). Similarly, the plasma rotation can stabilize or destabilize the RWM, depending on the plasma resistivity. These numerical results from MARS-F are qualitatively confirmed by an analytic theory model which includes the GGJ effect.

A conservative scheme of drift kinetic electrons for gyrokinetic simulation of kinetic-MHD processes in toroidal plasmas

A conservative scheme of drift kinetic electrons for gyrokinetic simulations of kinetic-magnetohydrodynamic processes in toroidal plasmas has been formulated and verified. Both vector potential and electron perturbed distribution function are decomposed into adiabatic part with analytic solution and non-adiabatic part solved numerically. The adiabatic parallel electric field is solved directly from the electron adiabatic response, resulting in a high degree of accuracy. The consistency between electrostatic potential and parallel vector potential is enforced by using the electron continuity equation. Since particles are only used to calculate the non-adiabatic response, which is used to calculate the non-adiabatic vector potential through Ohm's law, the conservative scheme minimizes the electron particle noise and mitigates the cancellation problem. Linear dispersion relations of the kinetic Alfvén wave and the collisionless tearing mode in cylindrical geometry have been verified in gyrokinetic toroidal code simulations, which show that the perpendicular grid size can be larger than the electron collisionless skin depth when the mode wavelength is longer than the electron skin depth.

Nonlinear co-existence of beta-induced Alfvén eigenmodes and beta-induced Alfvén-acoustic eigenmodes

The nonlinear co-existence of β–induced Alfvén eigenmode (BAE) and β–induced Alfvén-acoustic eigenmode (BAAE) is found in simulations using the gyrokinetic toroidal code, which provides a new mechanism responsible for BAAE excitation in tokamaks. Here, the normalized pressure β is the ratio between plasma thermal pressure and magnetic pressure. The nonlinear simulation results show that the BAAE branch emerges after the BAE branch is saturated. The mode structure's evolution shows that existence of BAAE will change the original BAE mode structure. The perturbed distribution functions in the velocity phase space show that a new resonance region manifesting the wave-particle resonance in the BAAE branch appears during the nonlinear co-existence stage.

Gyrokinetic particle simulations of the effects of compressional magnetic perturbations on drift-Alfvenic instabilities in tokamaks

The compressional component of magnetic perturbation δB∥ can play an important role in drift-Alfvenic instabilities in tokamaks, especially as the plasma β increases (β is the ratio of kinetic pressure to magnetic pressure). In this work, we have formulated a gyrokinetic particle simulation model incorporating δB∥, and verified the model in kinetic Alfven wave simulations using the Gyrokinetic Toroidal Code in slab geometry. Simulations of drift-Alfvenic instabilities in tokamak geometry shows that the kinetic ballooning mode (KBM) growth rate decreases more than 20% when δB∥ is neglected for βe=0.02, and that δB∥ has stabilizing effects on the ion temperature gradient instability, but negligible effects on the collisionless trapped electron mode. The KBM growth rate decreases about 15% when equilibrium current is neglected.

Effects of electron cyclotron current drive on magnetic islands in tokamak plasmas

The effects of the electron cyclotron current drive on magnetic islands in tokamak plasmas are studied using gyrokinetic simulations. By investigating the effects of different characteristics of the driven current, such as current density distribution and deposition location, the factors which can determine the suppression effect on the resistive tearing modes have been explored. It is found that an electron cyclotron wave (ECW) driven current with a larger peak value and more focused deposition region has a better stabilization effect. When the ECW-driven current is closer to the rational surface, it has a better stabilizing effect. These gyrokinetic toroidal code (GTC) linear simulations in the electron fluid limit of the tearing modes in the cylindrical geometry agree well with the magnetohydrodynamic codes. Furthermore, the optimal timing control of the current deposition on resistive tearing modes is demonstrated.

Modern gyrokinetic particle-in-cell simulation of fusion plasmas on top supercomputers

The gyrokinetic toroidal code at Princeton (GTC-P) is a highly scalable and portable particle-in-cell (PIC) code. It solves the 5-D Vlasov–Poisson equation featuring efficient utilization of modern parallel computer architectures at the petascale and beyond. Motivated by the goal of developing a modern code capable of dealing with the physics challenge of increasing problem size with sufficient resolution, new thread-level optimizations have been introduced as well as a key additional domain decomposition. GTC-P’s multiple levels of parallelism, including internode 2-D domain decomposition and particle decomposition, as well as intranode shared memory partition and vectorization, have enabled pushing the scalability of the PIC method to extreme computational scales. In this article, we describe the methods developed to build a highly parallelized PIC code across a broad range of supercomputer designs. This particularly includes implementations on heterogeneous systems using NVIDIA GPU accelerators and Intel Xeon Phi (MIC) coprocessors and performance comparisons with state-of-the-art homogeneous HPC systems such as Blue Gene/Q. New discovery science capabilities in the magnetic fusion energy application domain are enabled, including investigations of ion–temperature–gradient driven turbulence simulations with unprecedented spatial resolution and long temporal duration. Performance studies with realistic fusion experimental parameters are carried out on multiple supercomputing systems spanning a wide range of cache capacities, cache-sharing configurations, memory bandwidth, interconnects, and network topologies. These performance comparisons using a realistic discovery-science-capable domain application code provide valuable insights on optimization techniques across one of the broadest sets of current high-end computing platforms worldwide.