We present a GPU-accelerated, arbitrary-order, nearly quadrature-free, Runge-Kutta (RK) discontinuous Galerkin (DG) approach to interface capturing for atomizing multiphase flows via the conservative level set (CLS) method [16, 17] An arbitrary-order DG numerical method is utilized for both advection and reinitialization, further developing the ideas of Czajkowski and Desjardins  by implementing a quadrature-free approach allowing for arbitrary polynomial degree, and treating the normal function in a DG sense. For effective use of processing power, the method is executed with the dual narrow band overset mesh approach of the refined level set grid method . Computation is performed in parallel on either CPU or GPU architectures to make the method feasible at high order. Finally, by using normalized Legendre polynomial basis functions, we are able to pre-compute volume and surface integrals analytically. The resulting sparse integral arrays are stored in the compressed row storage format to take full advantage of parallelism on the GPU, where performance relies heavily on well-managed memory operations. The accuracy, consistency, and convergence of the resulting method is demonstrated using the method of manufactured solutions (MMS). Using MMS, we demonstrate k + 1 order spatial convergence for kth order normalized Legendre polynomial basis functions on both advection and reinitialization. MMS is also used to demonstrate the benefits of GPU hardware, where advection is found to provide a speedup factor >45x comparing a 2.0GHz Intel Xeon E5-2620 in serial against a NVIDIA Tesla K20c, with speedup increasing with polynomial degree. Arbitrarily high convergence rates combined with speedup factors that increase with polynomial degree motivate the development and use of a GPU accelerated, arbitrary-order DG method.