We investigate the system of partial differential equations used in resistive magnetohydrodynamic modeling of fusion plasmas. This system couples the Euler and Maxwell equations for evolution of a charged fluid in an electromagnetic field, hence the magnetic field in the resulting PDE system must evolve on a divergence-free constraint manifold. As traditional numerical solution approaches often violate these constraints, we investigate a reformulation of the resistive MHD system to allow for accurate evolution of the continuum-level equations, while simultaneously ensuring that the solution satisfies the solenoidal constraint.