  We present a weighted BFBT approximation (w-BFBT) to the inverse Schur
complement of a Stokes system with highly heterogeneous viscosity. When used as
part of a Schur complement-based Stokes preconditioner, we observe robust fast
convergence for Stokes problems with smooth but highly varying (up to 10 orders
of magnitude) viscosities, optimal algorithmic scalability with respect to mesh
refinement, and only a mild dependence on the polynomial order of high-order
finite element discretizations ($Q_k \times P_{k-1}^{disc}$, order $k \ge 2$).
For certain difficult problems, we demonstrate numerically that w-BFBT
significantly improves Stokes solver convergence over the widely used inverse
viscosity-weighted pressure mass matrix approximation of the Schur complement.
In addition, we derive theoretical eigenvalue bounds to prove spectral
equivalence of w-BFBT. Using detailed numerical experiments, we discuss
modifications to w-BFBT at Dirichlet boundaries that decrease the number of
iterations. The overall algorithmic performance of the Stokes solver is
governed by the efficacy of w-BFBT as a Schur complement approximation and, in
addition, by our parallel hybrid spectral-geometric-algebraic multigrid (HMG)
method, which we use to approximate the inverses of the viscous block and
variable-coefficient pressure Poisson operators within w-BFBT. Building on the
scalability of HMG, our Stokes solver achieves a parallel efficiency of 90%
while weak scaling over a more than 600-fold increase from 48 to all 30,000
cores of TACC's Lonestar 5 supercomputer.
