Numerical calculation of neoclassical distribution functions and current profiles in low collisionality, axisymmetric plasmas

A new code, the Neoclassical Ion-Electron Solver (NIES), has been written to solve for stationary, axisymmetric distribution functions (f) in the conventional banana regime for both ions and electrons using a set of drift-kinetic equations (DKEs) with linearized Fokker-Planck-Landau collision operators. Solvability conditions on the DKEs determine the relevant non-adiabatic pieces of f (called h). We work in a 4D phase space in which ψ defines a flux surface, θ is the poloidal angle, v is the magnitude of the velocity referenced to the mean flow velocity, and λ is the dimensionless magnetic moment parameter. We expand h in finite elements in both v and λ. The Rosenbluth potentials, Φ and Ψ, which define the integral part of the collision operator, are expanded in Legendre series in cos χ, where χ is the pitch angle, Fourier series in cos θ, and finite elements in v. At each ψ, we solve a block tridiagonal system for h _i (independent of f _e), then solve another block tridiagonal system for h e (dependent on f _i). We demonstrate that such a formulation can be accurately and efficiently solved. NIES is coupled to the MHD equilibrium code JSOLVER [J. DeLucia, J. Comput. Phys. 37, 183-204 (1980)] allowing us to work with realistic magnetic geometries. The bootstrap current is calculated as a simple moment of the distribution function. Results are benchmarked against the Sauter analytic formulas and can be used as a kinetic closure for an MHD code (e.g., M 3 D - C 1 [S. C. Jardin, Comput. Sci. Discovery 5, 014002 (2012)]).