  We extend the positivity-preserving method of Zhang & Shu (2010, JCP, 229,
3091-3120) to simulate the advection of neutral particles in phase space using
curvilinear coordinates. The ability to utilize these coordinates is important
for non-equilibrium transport problems in general relativity and also in
science and engineering applications with specific geometries. The method
achieves high-order accuracy using Discontinuous Galerkin (DG) discretization
of phase space and strong stability-preserving, Runge-Kutta (SSP-RK) time
integration. Special care in taken to ensure that the method preserves strict
bounds for the phase space distribution function $f$; i.e., $f\in[0,1]$. The
combination of suitable CFL conditions and the use of the high-order limiter
proposed in Zhang & Shu (2010) is sufficient to ensure positivity of the
distribution function. However, to ensure that the distribution function
satisfies the upper bound, the discretization must, in addition, preserve the
divergence-free property of the phase space flow. Proofs that highlight the
necessary conditions are presented for general curvilinear coordinates, and the
details of these conditions are worked out for some commonly used coordinate
systems (i.e., spherical polar spatial coordinates in spherical symmetry and
cylindrical spatial coordinates in axial symmetry, both with spherical momentum
coordinates). Results from numerical experiments --- including one example in
spherical symmetry adopting the Schwarzschild metric --- demonstrate that the
method achieves high-order accuracy and that the distribution function
satisfies the maximum principle.
