Recent progress on the development of scalable differential variational inequality multigrid-based solvers for the phase-field approach to mesoscale materials modeling is described. We have developed an active-set method for variational inequalities in PETSc, leveraging experience by the optimization community in TAO. A geometric multigrid solver in PETSc is used to solve the resulting linear systems. We present strong and weak scaling results for 2D coupled Allen-Cahn/Cahn-Hilliard systems.