The research introduces a GPU-accelerated, matrix-free multigrid-preconditioned conjugate-gradient solver tailored for SBP-SAT finite-difference discretizations of elliptic PDEs. By designing SBP-preserving multigrid operators and specialized matrix-free GPU kernels that efficiently handle boundary stencils, the method achieves near-constant iteration counts and significant runtime improvements, delivering up to 5x speedups over matrix-explicit implementations and outperforming PETSc and AmgX solvers on large problems while using far less memory