  Upon application of a sufficiently strong electric field, electrons break
away from thermal equilibrium and approach relativistic speeds. These highly
energetic runaway electrons (MeV) play a significant role in tokamak disruption
physics, and therefore their accurate understanding is essential to develop
reliable mitigation strategies. For this purpose, we have developed a fully
implicit solver for the 0D-2P (i.e., including two momenta coordinates)
relativistic nonlinear Fokker-Planck equation (rFP). As in earlier implicit rFP
studies (NORSE, CQL3D), electron-ion interactions are modeled using the Lorentz
operator, and synchrotron damping using the Abraham-Lorentz-Dirac reaction
term. However, our implementation improves on these earlier studies by 1)
ensuring exact conservation properties for electron collisions, 2) strictly
preserving positivity, and 3) being scalable algorithmically and in parallel.
Key to our proposed approach is an efficient multigrid preconditioner for the
linearized rFP equation, a multigrid elliptic solver for the Braams-Karney
potentials [Braams and Karney, Phys. Rev. Lett. 59, 16 (1987)], and a novel
adaptive technique to determine the associated boundary values. We verify the
accuracy and efficiency of the proposed scheme with numerical results ranging
from small electric-field electrical conductivity measurements to the accurate
reproduction of runaway tail dynamics when strong electric fields are applied.
