A two-scale approach for efficient on-the-fly operator assembly in massively parallel high performance multigrid codes

Large scale matrix-free finite element implementations save memory and are often significantly faster than implementations using classical sparse matrix techniques. They are especially well suited for massively parallel geometric multigrid solvers in combination with hierarchical hybrid grids on polyhedral domains. In the case of constant coefficients, the number of different stencil entries depends only on the coarse grid and does not increase with the number of refinement levels. However, for non-polyhedral domains the situation changes. Then even for the Laplace operator, the element mapping leads to fine grid stencils that can vary from grid point to grid point. Traditional matrix-free techniques that are based on an element-wise assembly then result in a considerably increase in computational cost. To compensate for this shortcoming, we introduce a new two-scale approach that uses a surrogate operator. It exploits a piecewise polynomial approximation of the entries of the stencil of the fine grid operator with respect to the coarse mesh size. The low-cost evaluation of these surrogate polynomials results in an efficient stencil assembly on-the-fly for non-polyhedral domains. We discuss and illustrate numerically two-scale a priori bounds. The accuracy of the approximate solution can be further improved if combined with a double discretization technique. A careful performance analysis in combination with a hardware-aware code optimization based on the Execution-Cache-Memory model yields a significant speed up. Weak and strong scaling results illustrate the potential of this new two-scale approach within large scale PDE simulations.