Implementation techniques for multigrid solvers for high-order Discontinuous Galerkin methods

Matrix-free geometric multigrid solvers for elliptic PDEs that have been discretised with Higher-order Discontinuous Galerkin (DG) methods are ideally suited to exploit state-of-the-art computer architectures. Higher polynomial degrees offer exponential convergence, while the workload fits to vector units, is straightforward to parallelise, and exhibits high arithmetic intensity. Yet, DG methods such as the interior penalty DG discreisation do not magically guarantee high performance: they require non-local memory access due to coupling between neighbouring cells and break down into compute steps of widely varying costs and compute character. We address these limitations by developing efficient execution strategies for $hp$-multigrid. Separating cell- and facet-operations by introducing auxiliary facet variables localizes data access, reduces the need for frequent synchronization, and enables overlap of computation and communication. Loop fusion results in a single-touch scheme which reads (cell) data only once per smoothing step. We interpret the resulting execution strategies in the context of a task formalism, which exposes additional concurreny. The target audience of this paper are practitioners in Scientific Computing who are not necessarily experts on multigrid or familiar with sophisticated discretisation techniques. By discussing implementation techniques for a powerful solver algorithm we aim to make it accessible to the wider community.