Solving advection equations with reduction multigrids on GPUs

Methods for solving hyperbolic systems typically depend on unknown ordering (e.g., Gauss-Seidel, or sweep/wavefront/marching methods) to achieve good convergence. For many discretisations, mesh types or decompositions these methods do not scale well in parallel. In this work we demonstrate that the combination of AIRG (a reduction multigrid which uses GMRES polynomials) and PMISR DDC (a CF splitting algorithm which gives diagonally dominant submatrices) can be used to solve linear advection equations in parallel on GPUs with good weak scaling. We find that GMRES polynomials are well suited to GPUs when applied matrix-free, either as smoothers (at low order) or as an approximate coarse grid solver (at high order). To improve the parallel performance we automatically truncate the multigrid hierarchy given the quality of the polynomials as coarse grid solvers. Solving time-independent advection equations in 2D on structured grids, we find 66-101% weak scaling efficiency in the solve and 47-63% in the setup with AIRG, across the majority of Lumi-G, a pre-exascale GPU machine.