We extend the positivity-preserving method of Zhang & Shu ( ) 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 ( ) 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 .