We present a custom geometric multigrid preconditioned conjugate gradient method that applies summation-by-parts(SBP)-preserving interpolations and a custom matrix-free GPU kernel that achieves up to 5x speedup compared to solvers from PETSc and AmgX