Page 7 - Feasibility study of microwave electron heating on the C-2 field-reversed configuration device

012509-7 Fulton et al.
Phys. Plasmas 23, 012509 (2016)
The difference between the magnetic geometries of the tokamak and FRC makes it necessary to reformulate the Poisson solver. In the gyrokinetic formulation, electric and magnetic fields are solved on the plane perpendicular to magnetic field lines. In tokamaks, the field is primarily toroidal with a small poloidal component, while in FRCs, the field is entirely poloidal. GTC has several methods for solving the Poisson equation,63 and for the FRC, we solve the gyrokinetic Poisson equation with the Pad e approxima- tion. The normalized form of the gyrokinetic Poisson equa- tion and its representation in FRC geometry may be expressed as
this scheme requires only minor modification of the parti- cle boundary conditions, which are already periodic in the full torus case. Particles passing out through 0 enter from 2p=n and vice versa. The mode can now be well resolved without sacrificing computational efficiency. The total number of particles, particles per wavelength, number of grid points per wavelength and, correspondingly, com- putational cost and mode resolution are all independent of the mode wavelength.
In future non-linear simulations, when inclusion of mul- tiple n modes is desirable, the partial torus domain may still be used to sample a single toroidal mode number plus its har- monics: 0, n, 2n, 3n,. The full torus domain will be applied in instances where all n modes are kept.
B. Gyroaveraging
For turbulence and transport, the plasma phenomena of interest often evolve on a time-scale longer than the ion gyro-period. The salient feature of gyrokinetic formula- tions is averaging over the gyro-phase angle, which reduces the phase space dimensionality from six dimen- sions to five dimensions and allows a much coarser simula- tion time step. In order to retain realistic finite Larmour radius effects, we must perform accurate gyro-averaging on particles.
In the FRC implementation, the gyro-average of a func- tion is split into radial gyro-averaging and toroidal gyro- averaging. In the toroidal direction, spatial sampling is expensive, because domain decomposition requires message passing between toroidal neighbors. Fortunately, when simu- lating single-n modes, mode variation in the toroidal direc- tion is simple and easy to approximate. By contrast, spatial sampling in the radial direction does not require expensive communication, but all radial wave numbers are retained making approximation of the radial mode variation impracti- cal. These qualitative differences lead us to treat the two dimensions separately.
Radial gyro-averaging is represented by gyro-particles sampling different locations of the gyro-ring. In the case where the radial and toroidal wavenumbers are similar in magnitude, kr kf, radial two-point averaging may be used to include gyroaveraged effects.64 In radially local simula- tions, where kr kf, radial gyro-averaging is omitted, since its effects are negligible. In both cases, toroidal gyro- averaging is represented as the multiplication of a Bessel function of the first kind J0 ðkf qc Þ, where the arguments are kf 1⁄4 n=R and qc 1⁄4 mv?=jqjB. Here, n is the toroidal mode number, R is the major radius measured from the machine axis to the particle position, m is the particle mass, q is the particle charge, v? is the particle velocity perpendicular to the magnetic field, and B is the magnetic field strength at the particle position. The value of the Bessel function is calcu- lated using an intrinsic Fortran function for each individual particle, and the multiplication of this factor is self- consistently applied to each particle in the accumulation of the charge density on the grid and in the time-advancement of particle positions. Timing tests indicate that use of the Fortran intrinsic Bessel function does not significantly affect
r2U1⁄4ð1 r2Þðdn dn!Þ; ??ie
(19) 1⁄4 g @w2þg @f2 U; (20)
ww @2 ff @2
where dni is the perturbed ion density, dne is the perturbed electron density, U is electrostatic potential, and gnanb rna rnb is the metric tensor element.63
Our initial focus is on the linear properties of FRC insta- bilities. In this work, we present simulations of single-n modes, where n, the toroidal mode number, is an integer related to the mode’s toroidal wave number by ktor 1⁄4 n=R0. Characterizing linear instabilities in the FRC requires param- eter scans consisting of many such single-n simulations. To make these scans feasible, we implement optimizations in the Poisson solver. The first optimization is the reduction of the toroidal simulation domain to a partial-torus which elimi- nates grid size dependence on the toroidal mode number as detailed in Section III A. The gyro-kinetic treatment of par- ticles is detailed in Section III B. In Section III C, the second optimization, the semi-spectral operator used to solve the gyrokinetic Poisson equation, is outlined.
A. Partial torus domain
To accurately resolve shorter toroidal wavelength modes, a finer simulation grid is necessary, incurring higher computational costs. We make use of the toroidal periodicity of single-n modes to attain a computational cost that is inde- pendent of the toroidal mode number, n, in the linear simula- tion of a single n mode.
The computational cost of a simulation is proportional to the number of grid points used. In the field solver, increased number of grid points means that the matrix rep- resenting the Laplacian operator becomes larger. For fixed particles per grid cell, the total number of particles in the simulation is also proportional to the number of grid points. In a full torus, the total number of grid points is a product of grid points per wavelength and wavelengths per torus. While grid points per wavelength determines how well the mode in question is resolved, wavelengths per torus adds computational cost without benefit. The perio- dicity of a single-n mode in the toroidal direction allows reduction of the domain size from a full torus, 1⁄20; 2p , to a partial torus, 1⁄20; 2p=n . This partial torus domain corre- sponds to one wavelength of the mode. Implementation of