  Fast and accurate integration of geodesics in Kerr spacetimes is an important
tool in modeling the orbits of stars and the transport of radiation in the
vicinities of black holes. Most existing integration algorithms employ
Boyer-Lindquist coordinates, which have coordinate singularities at the event
horizon and along the poles. Handling the singularities requires special
numerical treatment in these regions, often slows down the calculations, and
may lead to inaccurate geodesics. We present here a new general-purpose
geodesic integrator, GRay2, that overcomes these issues by employing the
Cartesian form of Kerr-Schild coordinates. By performing particular
mathematical manipulations of the geodesic equations and several optimizations,
we develop an implementation of the Cartesian Kerr-Schild coordinates that
outperforms calculations that use the seemingly simpler equations in
Boyer-Lindquist coordinates. We also employ the OpenCL framework, which allows
GRay2 to run on multi-core CPUs as well as on a wide range of GPU hardware
accelerators, making the algorithm more versatile. We report numerous
convergence tests and benchmark results for GRay2 for both time-like (particle)
and null (photon) geodesics.
