  We present augmented Lagrangian Schur complement preconditioners and robust
multigrid methods for incompressible Stokes problems with extreme viscosity
variations. Such Stokes systems arise, for instance, upon linearization of
nonlinear viscous flow problems, and they can have severely inhomogeneous and
anisotropic coefficients. Using an augmented Lagrangian formulation for the
incompressibility constraint makes the Schur complement easier to approximate,
but results in a nearly singular (1,1)-block in the Stokes system. We present
eigenvalue estimates for the quality of the Schur complement approximation. To
cope with the near-singularity of the (1,1)-block, we extend a multigrid scheme
with a discretization-dependent smoother and transfer operators from
triangular/tetrahedral to the quadrilateral/hexahedral finite element
discretizations $[\mathbb{Q}_k]^d\times \mathbb{P}_{k-1}^{\text{disc}}$, $k\geq
2$, $d=2,3$. Using numerical examples with scalar and with anisotropic
fourth-order tensor viscosity arising from linearization of a viscoplastic
constitutive relation, we confirm the robustness of the multigrid scheme and
the overall efficiency of the solver. We present scalability results using up
to 28,672 parallel tasks for problems with up to 1.6 billion unknowns and a
viscosity contrast up to ten orders of magnitude.
