<p>Scalable preconditioners for saddle point problems are essential to the solution of problems in geodynamics and beyond. Recent years have produced a wealth of research into efficient solvers for finite element methods. &#160;These solvers are also effective, however, for orthogonal-grid finite volume discretizations of saddle point problems, also know as "staggered grid" or "marker and cell (MAC)" methods. Perhaps, ironically, due to the highly-structured nature of these discretizations, the use of advanced solvers is stymied due to the lack of a uniform topological abstraction, which is required for most scalable solvers, such as geometric multigrid. &#160;We present new software to allow experimentation with and composition of these advanced solvers. &#160;We focus on variable-viscosity Stokes problems with discontinuous coefficient jumps. &#160;In particular, we attempt to demonstrate how the important know robust preconditioners may be employed, and how new variants may be experimented with. &#160;Important solvers are compositions of block factorizations and multigrid cycles. &#160;We demonstrate as many of these as possible, including triangular block preconditioners with nested multigrid solves, and monolithic multigrid solves with cellwise (Vanka) or field-based (Distributed Gauss-Seidel, Braess-Sarazin) smoothers. &#160;Implementations are provided as part of the PETSc library, using the new DMStag component, and examples from the StagBL library are also shown where appropriate. &#160;These tools are intended to help break down the barrier between cutting-edge research into advanced solvers (which is only becoming more complex, as multi-phase problems are further explored) and practical usage in geophysical research and production codes.</p>