  Fully realizing the potential of multigrid solvers often requires custom
algorithms for a given application model, discretizations and even regimes of
interest, despite considerable effort from the applied math community to
develop fully algebraic multigrid (AMG) methods for almost 40 years. Classic
geometric multigrid (GMG) has been effectively applied to challenging,
non-elliptic problems in engineering and scientifically relevant codes, but
application specific algorithms are generally required that do not lend
themselves to deployment in numerical libraries. However, tools in libraries
that support discretizations, distributed mesh management and high performance
computing (HPC) can be used to develop such solvers.
  This report develops a magnetohydrodynamics (MHD) code in PETSc (Portable
Extensible Toolkit for Scientific computing) with a fully integrated GMG solver
that is designed to demonstrate the potential of our approach to providing fast
and robust solvers for production applications. These applications must,
however, be able to provide, in addition to the Jacobian matrix and residual of
a pure AMG solver, a hierarchy of meshes and knowledge of the application's
equations and discretization. An example of a 2D, two field reduced resistive
MHD model, using existing tools in PETSc that is verified with a ``tilt"
instability problem that is well documented in the literature is presented and
is an example in the PETSc repository (\path{src/ts/tutorials/ex48.c}).
Preliminary CPU-only performance data demonstrates that the solver can be
robust and scalable for the model problem that is pushed into a regime with
highly localized current sheets, which generates strong, localized
non-linearity, that is a challenge for iterative solvers.
