Scalable hybrid multigrid for staggered grid discretizations in geodynamics

The staggered grid finite difference is a robust discretization method for the high-resolution 3D  geodynamic simulations that involve heterogeneous material parameters. Achieving its scalability on parallel machines inevitably requires the application of multigrid solvers. In particular, a coupled velocity-pressure geometric multigrid preconditioner based on Galerkin coarsening scheme demonstrates very good results. However, this method relies on assembled matrices which have a significant memory imprint and prohibits achieving peak performance due to suboptimal use of the limited memory bandwidth. A geometric multigrid method, based on the re-discretization of linear operators on the coarser levels, converges generally slower, but can be implemented  in a completely matrix-free manner. It poses a valuable alternative to the Galerkin method, since an increased number of iterations can be compensated by a greater performance of the matrix-vector products computed on the fly without storing matrices in the memory.

Here we present a hybrid approach that allows optimal combination between various types of coarsening techniques for staggered grid discretizations. This work is performed within the framework of ChEESE-2p project (Centre of Excellence for Exascale in Solid Earth) and involves the flagship code LaMEM (Lithosphere and Mantle Evolution Model), which is based on Portable Extensible Toolkit for Scientific computation (PETSc), following an approach suggested for finite element discretizations by May et al. (2015). Here, we extend it to the staggered grid finite difference, discuss the optimal solver parameter selection, and document performance gains that can be achieved by using the matrix-free operators.

We typically start with a few levels of re-discretized matrix-free operators, followed by Galerkin geometric coarsening operating on assembled matrices. This approach ensures that most of the optimization and memory saving is already obtained at the top levels, whereas more robust Galerkin coarsening can be used at coarser levels without compromising the convergence. At the coarse grid level, we either utilize a parallel sparse direct solver or a black-box algebraic multigrid method.  The number of processors participating in a coarse grid solve can be optimally selected via PETSc sub-communicator framework (Telescope). 

D. A. May, J. Brown, L. Le Pourhiet, 2015. A scalable matrix-free multigrid preconditioner for finite element discretizations of heterogeneous Stokes flow, Comput. Methods Appl. Mech. Engrg., 290, 496–523.