Numerical calculation of neoclassical distribution functions in an axisymmetric torus

We 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 complete Landau collision operators. Solubility conditions on the DKEs determine the relevant non-Maxwellian pieces of $f$ (called $f_{NM}$). We work in a 4D phase space in which $\psi$ defines a flux surface, $\theta$ is the poloidal angle, $v$ is the total velocity, and $\lambda$ is the pitch angle parameter. We expand $f_{NM}$ in finite elements in both $v$ and $\lambda$. The Rosenbluth potentials, $\Phi$ and $\Psi$, which define the collision operator, are expanded in Legendre series in $\cos \chi$, where $\chi$ is the pitch angle, Fourier series in $\cos \theta$, and finite elements in $v$. At each $\psi$, we solve a block tridiagonal system for $f_{NMi}$ (independent of $f_e$), then solve another block tridiagonal system for $f_{NMe}$ (dependent on $f_i$). We demonstrate that such a formulation can be accurately and efficiently solved. Results will be benchmarked against other codes (e.g., NEO) and could be used as a kinetic closure for an MHD code (e.g., M3D-C1). Results will also include the lowest-order collisionality correction and the use of generalized Grad-Shafranov equilibria.