Efficient Distributed Matrix-free Multigrid Methods on Locally Refined Meshes for FEM Computations

Which multigrid approach to choose depends on the way the mesh is generated and on the underlying finite-element space. If the mesh is generated by globally refining each cell of a coarse grid recursively, then it is a natural choice to apply geometric multigrid (abbreviated here as h-multigrid), which uses the levels of the resulting mesh hierarchy as multigrid levels. Alternatively, in the context of high-order finite elements, it is possible to create levels by reducing the polynomial order of the shape functions p of the elements, while keeping the mesh the same, as done by polynomial multigrid (abbreviated as p-multigrid). For a very fine, unstructured mesh with low-order elements, it is not as trivial to explicitly construct enough multigrid levels and one might need to fall back to non-nested multilevel algorithms [1, 22, 23] or algebraic multigrid (AMG; see the review by Stüben [97]). These basic multigrid strategies can be nested in hybrid multigrid solvers [32, 71, 84, 85, 89, 93, 98] and, in doing so, one can exploit the advantages of all of them regarding robustness. Most common are hp-multigrid, which combines h- and p-multigrid, and AMG as black-box coarse-grid solver of geometric or polynomial multigrid solvers.

All the above-mentioned multigrid variants are applicable to locally refined meshes. However, local refinement comes with additional options (local vs. global definition of the multigrid levels, resulting in local(-smoothing) and global(-coarsening) versions of geometric/polynomial multigrid) and with additional challenges, which are particularly connected with the presence of hanging-node constraints, as indicated in Figure 1.

Some authors of this study have been involved in previous contributions to the research field of multigrid methods. Their implementations are used and extended in this work.

In References [27, 28, 64, 66], matrix-free and parallel implementations of geometric local-smoothing algorithms from the deal.II finite-element library were investigated and compared with AMG, demonstrating the benefits of using matrix-free implementations for both CPU and GPU but also non-optimal scalability due to load imbalance for simple partitioning strategies. Local-smoothing algorithms have been a core module of deal.II for many years; however, they gained much maturity, flexibility, and performance due to the above-mentioned thorough studies. In Reference [66], their performance was compared to a state-of-the-art implementation of a Poisson solver from the literature [36]: A speedup of 30% for quadratic elements and a speedup of a factor of 8 for fourth-degree elements on comparable hardware underline the high node-level performance of the code.

Independently of local smoothing, an efficient hybrid multigrid solver for discontinuous Galerkin methods (DG) for globally refined meshes was developed and presented in Reference [32]. It is part of the deal.II-based incompressible Navier–Stokes solver ExaDG [10] and relies on auxiliary-space approximation [5], i.e., on the transfer into a continuous space, as well as on subsequent execution of p-multigrid, h-multigrid, and AMG. That solver was extended to leverage locally refined meshes, using a geometric/polynomial global-coarsening algorithm after the transfer into the continuous space, to simulate the flow through a lung geometry [60]. Its functionalities have been generalized and are the foundations of the new global-coarsening implementation in deal.II [7, 9], targeting both h- and p-multigrid, in addition to the established geometric local-smoothing implementation.

In this publication, we consider three well-known multigrid algorithms for locally refined meshes for continuous higher-order matrix-free FEM: geometric local smoothing, geometric global coarsening (both h-multigrid), and polynomial global coarsening (a variant of p-multigrid). We have implemented them into the same framework, which allows us to compare their implementation complexity and performance for a large variety of problem sizes. This has not been done in an extensive way in the literature, often using only one of them [28, 99]. Furthermore, we rely on matrix-free operator evaluation, which provides optimal, state-of-the-art implementations in terms of node-level performance on modern hardware [66] and hence exercise the methods in a challenging context in terms of communication costs and workload differences.

The algorithms presented in this publication have been integrated into the open source finite-element library deal.II [8] and are part of its 9.4 release [9]. All experimental results have been obtained with small benchmark programs leveraging on the infrastructure of deal.II. The programs are available on GitHub under https://github.com/peterrum/dealii-multigrid.

The results obtained in this publication for continuous FEM are transferable to the DG case, where one does not have to consider hanging-node constraints but fluxes between differently refined cells. In the case of auxiliary-space approximation [60], this difference only involves the finest level and the rest of the multigrid algorithm could be as described in this publication.

The remainder of this work is organized as follows. In Section 2, we give a short overview of multigrid variants applicable to locally refined meshes. Section 3 presents implementation details of our solver, and Section 4 discusses relevant performance metrics. Sections 5 and 6 demonstrate performance results for geometric multigrid and polynomial multigrid, and, in Section 7, the solver is applied to a challenging Stokes problem. Finally, Section 8 summarizes our conclusions and points to further research directions.

Geometric local-smoothing algorithms [19, 24, 28, 51, 53, 54, 70, 77, 90, 96] use the refinement hierarchy also for multigrid levels and perform smoothing refinement level by refinement level: Cells of less refined parts of the mesh are skipped (see the top part of Figure 1) so that hanging-node constraints do not need to be considered during smoothing. The authors of References [3, 20, 50, 77, 100, 102] investigate a version of local smoothing in which smoothing is also performed on a halo of a single coarse cell, which includes hanging-node constraints. We will not consider this form of local smoothing in the following.

The fact that domains on each level might not cover the whole computational domain results in multiple issues. Data need to be transferred in lines 2 and 13 in Algorithm 1 to and from all multigrid levels that are active, i.e., have cells that are not refined further. Moreover, internal interfaces (also known as “refinement edges”) can appear with the need for special treatment. For details, interested readers are referred to References [51, 64]. In the following, we summarize key aspects relevant for our investigations.

For the purpose of explanation, let us split the degrees of freedom (DoFs) associated with the cells on an arbitrary level l into the interior ones \(\boldsymbol {x}^{(l)}_S\) and the ones at the refinement edges \(\boldsymbol {x}^{(l)}_E\), leading to the following block structure in the associated matrix system \(\boldsymbol {A}^{(l)} \boldsymbol {x}^{(l)} = \boldsymbol {b}^{(l)}\):For presmoothing on level l, homogeneous Dirichlet boundary conditions are applied at the refinement edges (\(\boldsymbol {x}^{(l)}_{E} = \boldsymbol {0}\)), allowing us to skip the coupling matrices. However, when switching to a coarser or finer level, the coupling matrices need to be considered. The residual to be restricted becomes the following:Since \(\boldsymbol {b}^{(l)}_{E}\) has already been transferred to the coarser level by line 2 in Algorithm 1, one only has to restrict and add the result of term \(*\) to the coarser level. During postsmoothing, an inhomogeneous Dirichlet boundary condition is applied with boundary values prescribed by the coarser level. This can be achieved, e.g., by modifying the right-hand side (\(\boldsymbol {b}^{(l)}_{S} \leftarrow \boldsymbol {b}^{(l)}_{S} - \boldsymbol {A}^{(l)}_{SE} \boldsymbol {x}^{(l)}_{E}\)).

A natural choice to partition the multigrid levels for local-smoothing algorithms is to partition the active level and let cells on lower refinement levels inherit the rank of their children. A simple variant of it is the “first-child policy” (see Figure 2(a) or Reference [28]): It recursively assigns the parent cell the rank of its first child cell. Since in adaptive FEM codes the parents of locally owned and ghost cells are generally already available on processes due to treelike data-structure storage of adaptively refined meshes [15, 26], no additional data structures need to be constructed and saved, but the storage of an additional flag (multigrid rank of the cell) is enough, leading to low memory consumption. Furthermore, intergrid transfer operations are potentially cheap as data are mainly transferred locally. A disadvantage—besides having to consider the edge constraints—is the potential load imbalance on the levels, as discussed in Reference [28]. This load imbalance could be alleviated by partitioning each level for itself. Such an alternative partitioning algorithm would lead to similar problems as in the case of global-coarsening algorithms (discussed next) and, as a result, the needed data structures would become more complex. This would counter the claimed simplicity of the data structures of this method; hence, we will consider geometric local smoothing only with “first-child policy” in this publication.

Furthermore, the fact that the transfer to and from the multigrid levels involves all active levels prevents an early switch to a coarse-grid solver on levels finer than those that are indeed not locally refined anymore, i.e., \(\boldsymbol {A}^{(l)}_{SS} = \boldsymbol {A}^{(l)}\). However, the absence of hanging nodes allows the usage of smoothers that have been developed for uniformly refined meshes, e.g., patch smoothers [11, 55], which are superior for anisotropic meshes.

Since geometric local smoothing is the only local-smoothing approach we will consider here, we will call it simply—as common in the literature—local smoothing in the following.

Geometric global-coarsening algorithms [20, 21] coarsen all cells simultaneously, translating to meshes with hanging nodes also on coarser levels of the multigrid hierarchy (see the bottom part of Figure 1). The computational complexity—i.e., the total number of cells to be processed—is slightly higher than in the case of local smoothing and might be non-optimal for some extreme examples of meshes [19, 54].

The fact that all levels cover the whole computational domain has the advantage that no internal interfaces need to be considered and the transfer to/from the multigrid levels becomes a simple copy operation to/from the finest level. However, hanging nodes have to be considered during the application of the smoothers on the levels. While this typically does not require new algorithms for basic smoothers as the infrastructure to support adaptive meshes can be simply applied to coarser representations, the operator evaluation and the applicable smoothers might become more expensive per cell by the need to resolve hanging-node constraints [81]. However, global-coarsening approaches show—for comparable smoothers—a better convergence behavior, which improves with the number of smoothing iterations [13, 14].

As the work on the levels generally increases compared to local smoothing, it is a valid option to repartition each level separately (see Figure 2(b)). On the one hand, this implies a higher pressure on the transfer operators, since they need to transfer data between independent meshes,¹ requiring potentially complex internal data structures, which describe the connectivities, and involved setup routines.² On the other hand, it opens the possibility to control the load balance between processes and the minimal granularity of work per process (by removing processes on the coarse level in a controlled way, allowing us to switch to subcommunicators [99]) as well as to apply a coarse-grid solver on any level. For the same reason, the construction of full multigrid solvers (FMG), which visit the finest level only a few times, is easier, since any level can be used as a starting point for the recursion. Not selecting the actual coarsest grid but some more refined level allows us to relax the unfavorable complexity of FMG in terms of \(\mathcal {O}(L^2)\) by reducing the time spent on the coarse levels.

Polynomial global-coarsening algorithms [12, 17, 18, 25, 29, 33, 35, 37, 38, 40, 41, 43, 44, 46, 48, 49, 52, 67, 72, 73, 74, 75, 76, 83, 85, 87, 88, 91, 94, 95, 99] are based on keeping the mesh size h constant on all levels but reducing the polynomial degree p of shape functions, e.g., to \(p=1\). Hence, the multigrid levels in this case have the same mesh but different polynomial orders. There are various strategies to reduce the order of the polynomial degree [32, 42]: The most common is the bisection strategy, which repeatedly halves the degree \(p^{(c)}= \lfloor p^{(f)} / 2 \rfloor\). This strategy reduces the number of DoFs in the case of a globally refined mesh similarly to the geometric multigrid strategies and provides a compromise between the number of V-cycles necessary to reduce the residual norm below a desired threshold and the cost of a single V-cycle [32].

The statements made in Section 2.2 about geometric global coarsening are also valid for polynomial global coarsening. However, the number of unknowns is reduced uniformly on each subdomain in contrast to geometric global coarsening, which obviates repartitioning of the level grids. This leads to a transfer operation that mainly works on locally owned DoFs.

In the following, we call geometric global coarsening simply global coarsening and polynomial global coarsening polynomial coarsening. As a reference, polynomial local smoothing arises, e.g., in the context of p- or hp-refinement, where only cells with the finest degree are smoothened [79].

Constraints need to be considered—with slight differences—in the case of both local-smoothing and global-coarsening algorithms. First, we impose Dirichlet boundary conditions in a strong form and express them as constraints. Second, hanging-node constraints, which force the solution of the refined side to match the polynomial expansion on the coarse side, need to be considered to maintain \(H^1\) regularity of the tentative solution [92]. In a general way, these constraints can be expressed as \(x_i = \sum _j c_{ij} x_j + b_i\),

where \(x_i\) is a constrained DoF, \(x_j\) a constraining DoF, \(c_{ij}\) the coefficient relating the DoFs, and \(b_i\) a real value, which can be used to consider inhomogeneities. We do not eliminate constraints but use a condensation approach [16, 101].

Instead of assembling the sparse system matrix \(\boldsymbol {A}\) and performing matrix-vector multiplications of the form \(\boldsymbol {A} \boldsymbol {x}\), the matrix-free operator evaluation computes the underlying finite-element integrals to represent \(\mathcal {A}(\boldsymbol {x})\). This is motivated by the fact that many iterative solvers and smoothers do not need the matrix \(\boldsymbol {A}\) explicitly but only its action on a vector. On a high level, matrix-free operator evaluation can be derived in two steps: (1) loop switching \(\boldsymbol {v} = \boldsymbol {A} \boldsymbol {x} = (\sum _{e}(\mathcal {G}_e^{T} \circ \mathcal {C}_e^{T}) \boldsymbol {A}_e (\mathcal {C}_e \circ \mathcal {G}_e)) \boldsymbol {x} = \sum _{e} \mathcal {G}_e^{T} \circ \mathcal {C}_e^{T} (\boldsymbol {A}_e (\mathcal {C}_e \circ \mathcal {G}_e \boldsymbol {x}))= \sum _{e} \mathcal {G}_e^T \circ \mathcal {C}_e^{T} (\boldsymbol {A}_e \boldsymbol {x}_e)\), i.e., replacing the sequence “assembly–mat-vec” by a loop over cells with the three steps “gathering–application of the element stiffness matrix–scattering,” and (2) exploitation of the structure of the element stiffness matrix \(\boldsymbol {A}_e\). In this formula, the operator \(\mathcal {C}_e \circ \mathcal {G}_e\) represents the assembly of element-wise quantities into the global matrix and vectors, respectively [92], including the constraints \(\mathcal {C}_e\) as discussed in Section 3.1. The element stiffness matrix \(\boldsymbol {A}_e\) is constructed from three main ingredients in a finite element method, the shape functions and derivatives of the tentative solution tabulated at quadrature points, geometric factors and coefficients of the underlying differential operator at quadrature points, and the multiplication by the test functions and their derivatives as well as subsequent summation over quadrature points. The final structure of a matrix-free operator evaluation in the context of continuous finite elements in operator form is as follows:This structure is depicted in Figure 3. For each cell e, cell-relevant values are gathered with operator \(\mathcal {G}_e\), constraints are applied with \(\mathcal {C}_e\), and values, gradients, or Hessians associated to the vector \(\boldsymbol {x}\) are computed at the quadrature points with \(\mathcal {S}_e\). These quantities are processed by a quadrature-point operation \(\mathcal {Q}_e\); the result is integrated and summed into the result vector \(\boldsymbol {v}\) by applying \(\tilde{\mathcal {S}}_e^T\), \(\mathcal {C}_e^T\), and \(\mathcal {G}_e^T\). In this publication, we consider symmetric (self-adjoint) PDE operators with \(\tilde{\mathcal {S}}_e = \mathcal {S}_e\).

In the literature, specialized implementations for GPUs [4, 57, 64, 68, 69] and CPUs [4, 62, 63, 80, 82] for operations as expressed in Equation (3) have been presented. For tensor-product (quadrilateral and hexahedral) elements, a technique known as sum factorization [78, 86] is often applied, which allows us to replace full interpolations from the local solution values to the quadrature points by a sequence of one-dimensional (1D) steps. In the context of CPUs, it is an option to vectorize across elements [62, 63], i.e., perform \((\mathcal {S}^T \circ \mathcal {Q} \circ \mathcal {S})_e\) for multiple cells in different lanes of the same instructions. This necessitates the data to be laid out in a struct-of-arrays fashion. The reshuffling of the data from array-of-structs format to struct-of-arrays format and back can be done, e.g., by \(\mathcal {G}_e\), while looping through all elements [62]. Note that our implementation performs \(\mathcal {C}_e \circ \mathcal {G}_e \boldsymbol {x}\) for a single SIMD batch of cells at a time to keep intermediate results in caches and reduce global memory traffic. For the application of hanging-node constraints, we use the special-purpose algorithm presented in Reference [69, 81], which is based on updating the DoF map \(\mathcal {G}_e\) and applying in-place sum factorization for the interpolation during the application of edge and face constraints. Even though this algorithm is highly optimized and comes with only a small overhead for high-order finite elements (\(\lt 20\%\) per cell with hanging nodes), the additional steps are not for free particularly for linear elements. This might lead to load imbalances in a parallel setting when some processes have a disproportionately high number of cells with hanging nodes, see the quantitative analysis in Reference [81].

Despite relying on matrix-free algorithms overall, some of the multigrid ingredients (e.g., smoother and coarse-grid solver) need an explicit representation of the linear operator in the form of a matrix or a part of it (e.g., its diagonal). Based on the matrix-free algorithm (3) applied to the local unit vector \(e_i\), the local matrix can be computed column by column:The resulting element matrix can be assembled as usual, including the resolution of the constraints. Computing the diagonal is slightly more complex when attempting to immediately write into a vector for the diagonal without complete matrix assembly. In our code, we choose to compute the jth entry of the locally relevant diagonal contribution viai.e., we apply the local constraint matrix \(\boldsymbol {C}_e\) from the left to the ith column of the element matrix (computed via Equation (4)) and apply \(\boldsymbol {C}_e\) again from the right. This approach needs as many basis vector applications as there are shape functions per cell. The local result can be simply added to the global diagonal via \(d = \sum _e \mathcal {G}^T_e d_e\). For cells without constrained DoFs, Equation (5) simplifies to \(d_e(i)=\boldsymbol {A}_e(i,i)\).

As the derivation of good matrix-free smoothers is an active research topic, we use Chebyshev iterations around a point-Jacobi method for smoothing [2]. This setup needs an operator evaluation and its diagonal representation according to Section 3.2. It is run on the levels defined by either the global-coarsening or the local-smoothing multigrid algorithm. In the case of local smoothing, refinement edges can be treated as homogeneous Dirichlet boundaries under the assumption that the right-hand-side vector is suitably modified.

The algorithms described in Section 3.2 allow us to set up traditional coarse-grid solvers, e.g., a Jacobi solver, a Chebyshev solver, and direct solvers, but also AMG. In this publication, we mostly apply AMG as a coarse-grid solver, since we use it either on very coarse meshes (in this case, it falls back to a direct solver) or for problems discretized with linear elements, for which AMG solvers are competitive [66].

The prolongation operator \(\mathcal {P}^{(f,c)}\) prolongates the result \(\boldsymbol {x}\) from a coarse space c to a fine space f (between either different geometric grids or different polynomial degrees),According to the literature [14, 99], this can be done in three steps:with \(\mathcal {C}^{(c)}\) setting the values of constrained DoFs on the coarse mesh, particularly resolving the hanging-node constraints, \(\tilde{\mathcal {P}}^{(f,c)}\) performing the prolongation on the discontinuous space as if no hanging nodes were present, and the weighting operator \(\mathcal {W}^{(f)}\) zeroing out the DoFs constrained on the fine mesh.

To derive a matrix-free implementation, one can express Equation (6) for nested meshes as loops over all (coarse) cells (see also Figure 3):Here, \(\mathcal {C}^{(c)}_e \circ \mathcal {G}^{(c)}_e\) gathers the cell-relevant coarse DoFs and applies the constraints just as in the case of the matrix-free operator evaluation (3). The operator \(\mathcal {P}^{(f,c)}_e\) embeds the coarse-space solution into the fine space, whereas \(\mathcal {S}^{(f)}_e\) sums the result back to a global vector. Since multiple cells could add to the same global entry of the vector \(\boldsymbol {x}^{(f)}\) during the cell loop, the values to be added have to be weighted with the inverse of the valence of the corresponding DoF. This is done by \(\mathcal {W}^{(f)}_e\), which also ignores constrained DoFs (zero valence) to be consistent with Equation (6). Figure 4 shows, as an example, the values of \(\mathcal {W}^{(f)}\) for a simple mesh for a scalar Lagrange element of degree \(1\le p \le 3\).

We construct the element prolongation matrices aswhere \(\phi ^{(c)}\) are the shape functions on the coarse cell and \(\phi ^{(f)}\) on the fine cell. Instead of treating each “fine cell” on its own, we group direct children of the coarse cells together and define the fine shape functions on the appropriate subregions. As a consequence, the “finite element” on the fine mesh depends both on the actual finite-element type (like on the polynomial degree and continuity) and on the refinement type, as indicated in Figure 5. For local smoothing, there is only one refinement configuration, since all cells on a finer level are children of cells on the coarser level. In the case of polynomial global coarsening, the mesh stays the same and only the polynomial degree is uniformly changed for all cells. In the case of geometric global coarsening, cells on a finer level are either identical to the cells or direct children of cells on the coarser level, but the element and its polynomial degree stay the same. This necessitates two prolongation cases, an identity matrix for non-refined cells and the coarse-to-fine embedding. We define the set of all coarse–fine cell pairs connected via the same element prolongation matrix as category \(\mathcal {C}\).

Since \(\mathcal {P}^{(f,c)}_e=\mathcal {P}^{(f,c)}_{\mathcal {C}(e)}\), i.e., all cells of the same category \(\mathcal {C}{(e)}\) have the same element prolongation matrix,

and to be able to apply them for multiple elements in one go in a vectorization-over-elements fashion [62] as in the case of matrix-free loops (3), we loop over the cells type by type so that Equation (7) becomes

We choose the restriction operator as the transpose of the prolongation operator,

We conclude this subsection with discussing the appropriate data structures for transfer operators, beginning with the ones needed for both global geometric and polynomial coarsening. Since global-coarsening algorithms smoothen on the complete computational domain, data structures only need to be able to perform two-level transfers (8) independently between arbitrary fine (f) and coarse (c) grids. \(\mathcal {C}^{(c)}_e \circ \mathcal {G}^{(c)}_e\) is identical to \(\mathcal {C}_e \circ \mathcal {G}_e\) in the matrix-free loop (3) so that specialized algorithms and data structures [81] can be applied and reused. \(\mathcal {S}^{(f)}_e\) needs the indices of DoFs for the given element e to be able to scatter the values, and \(\mathcal {W}^{(f)}_e\) stores the weights of DoFs (or of geometric entities—see also Figure 4) for the given element e. \(\mathcal {P}^{(f,c)}_{\mathcal {C}(e)}\) (and \(\mathcal {R}^{(c,f)}_{\mathcal {C}(e)}\)) need to be available for each category. We regard them as general operators and choose the most efficient way of evaluation based on the element types: simple dense-matrix representation vs. some substructure with sum factorization based on smaller 1D prolongation matrices. The category of each cell has to be known for each cell.

Alongside these process-local data structures, one needs access to all constraining DoFs on the coarse level required during \(\mathcal {C}^{(c)}_e \circ \mathcal {G}^{(c)}_e\) and to the DoFs of all child cells on the fine level during the process-local part of \(\mathcal {S}^{(f)}_e\), which is concluded by a data exchange. If external vectors do not allow access to the required DoFs, then we copy data to/from internal temporal global vectors with appropriate ghosting [62]. We choose the coarse-side identification of cells due to its implementation simplicity and structured data access at the price of more ghost transfer.³ For setting up the communication pattern, we use consensus-based sparse dynamic algorithms [6, 47].

For the sake of separation of concerns, one might create three classes to implement a global-coarsening transfer operator as we have done in deal.II. The relation of these classes is shown in Figure 6: The multigrid transfer class (MGTransferGlobalCoarsening) delegates the actual transfer tasks to the right two-level implementation (MGTwoLevelTransfer), which performs communications needed as well as evaluates (8) for each category and cell by using category-specific information from the third class (MGTransferSchemes).

The discussed data structures are applicable also to local smoothing. The fact that only a single category is available in this case allows us to simplify many code paths. However, one needs data structures that match DoF indices on the active level and on all multigrid levels that share a part of the mesh with the active level in addition to the two-level data structures. For further details on implementation aspects of transfer operators for local smoothing, see References [28, 64] and the documentation of MGTransferMatrixFree class in deal.II.

In the experimental Sections 5 and 6, we consider two types of static 3D meshes, as has been also done in Reference [28]. They are obtained by refining a coarse mesh consisting of a single cell defined by \([-1,1]^3\) according to one of the following two solution criteria:

These two meshes are relevant in practice, since similar meshes occur in simulations of flows with far fields and of multi-phase flows with bubbles [59] or any kind of interfaces. All the refinement procedures are completed by a closure after each step, ensuring one-irregularity in the sense that two leaf cells may only differ by one level if they share a vertex. Figure 7 shows the considered meshes and provides numbers regarding the cell count for \(3\le L \le 13\). All meshes are partitioned along space-filling curves [15, 26] with the option to assign cells weights.

Tables 1 and 2 give—as examples—evaluated numbers for geometrical metrics of the two considered meshes for a single process and for 192 processes with cells constrained by hanging nodes weighted with the factor of 2 for partitioning, compared to the rest of the cells. For a single process, only workload and memory consumption are shown. In the tables, bold indicates the most beneficial behavior among the listed variants.

Starting with the octant case, one can see that the serial workload in the case of local smoothing and global coarsening is similar, with local smoothing having consistently less work. The workload of each level is depicted in Figure 8. The behavior of the memory consumption is similar to the one of the workload: Global coarsening has a slightly higher memory consumption, since it explicitly needs to store the coarser meshes as well; however, the second finest mesh has already approximately one eighth of the size of the coarsest triangulation so that the overhead is small. In the parallel case, the memory consumption differences are higher: This is related to the fact that—in the case of global coarsening—overlapping ghosted forests of trees need to be saved. In contrast, the workload in the parallel case is much lower in the case of global coarsening: While global coarsening is able to reach efficiencies higher than \(90\%\), the efficiency is only approximately 50–60% in the case of local smoothing. The high value in the global-coarsening case was to be expected, since we repartition each level during construction. The low value in the case of local smoothing can be seen in Figure 9, where the minimum, maximum, and average workload are shown for each level. The workload on the finest level is optimally distributed between processes with cells; however, there are some processes without any cells, i.e., any work, on that level. On the second level, most processes can reduce the number of cells nearly optimally by a factor of 8, but processes idle on the finest level start to participate in the smoothing process with a much higher number of cells, similarly to what other processes process on the finest level. This pattern of discrepancy between the minimum and the maximum number of cells on the lower levels continues and, as a consequence, the maximum workload per level is higher, increasing the overall parallel workload. Figure 9 might also give the impression of a load imbalance in the case of global coarsening, since the minimum workload on the finest level is the half of the maximum value. This is related to the way we partition—penalizing cells that have hanging nodes with a weight of 2—and to the fact that a lot of cells with hanging nodes are clustered locally. The actual resulting load imbalance is small, as shown in Figure 12. The difference between local smoothing and global coarsening in the parallel workload comes at a price. While the vertical communication efficiency is—by construction—high in the local-smoothing case, this is not true in the case of global coarsening: Twenty percent or less is not uncommon, requiring the permutation of data during transfer. The horizontal communication costs are similar in the case of local smoothing and global coarsening. Assuming that (1) pre- and postsmoothing are the most time-consuming steps, (2) the parallel workload is the relevant metric, and (3) the transfer between levels is not dominant in this case, the values indicate that global coarsening might be twice as fast as local smoothing. However, if the transfer is the bottleneck, then the picture might look differently so that a conclusive statement only based on geometrical metrics is not possible in this case.

The observations for the octant case are also made, in a more pronounced form, in the shell case. Here, the vertical communication is more favorable in the case of global coarsening, and the workload efficiency is significantly worse in the local-smoothing case so that one can expect a noticeable speedup when using global coarsening.

In summary, one can state the following: In the serial case, global coarsening has to process at least as many cells as local smoothing so that one can expect that the latter has—with the assumption that the number of iterations is the same—an advantage regarding throughput, particularly since no hanging-node constraints have to be applied. With an increasing number of processes, the workload might not be well distributed in the case of local smoothing, leading to a drop in parallel efficiency. In contrast, this is—by construction—not an issue in the global-coarsening case, since the work is explicitly redistributed between the levels. The price is that one might need to send around a lot of data during restriction and prolongation. The number of levels in the case of local smoothing and global coarsening is the same, leading to potentially same scaling limits \(\mathcal {O}(L)\). However, with appropriate partitioning of the levels, one could decrease the number of participating processes on the coarser levels and switch to a coarse-grid solver at an earlier stage in the global-coarsening case, leading to a better scaling limit due to lower latency. In the following, we show experimentally that the above statements can be verified and take a more detailed look at the influence of the mutually contradicting requirements “good workload balance” and “cheap transfer” on the parallel efficiency.

Making definite general conclusions is difficult as they depend on the number of processes as well as on the type of the coarse mesh and of the refinement. While the two meshes considered here are prototypical for many problems we have encountered in practice, it is clear that the statements made in this publication cannot hold in all cases, since it is easy to construct examples that favor one multigrid variant over the other. However, the two meshes clearly demonstrate the dominating aspects at various problem sizes; the actual crossover point, however, might be problem- and hardware-specific.

Figure 10 gives an overview of the time shares during the solution process in a serial shell simulation for local smoothing and global coarsening. Without going into details of the actual numbers, one can see that most of the time is spent in the multigrid preconditioner in the case of both local smoothing and global coarsening (72%/76%). It is followed by the other operations in the outer conjugate-gradient solver (21%/22%). The least time is spent for transferring data between preconditioner and solver (7%/2%). The higher time of the transfer in the local-smoothing case is not surprising, since the transfer involves all multigrid levels sharing cells with the active level. Since the overwhelming share of solution time is taken by the multigrid preconditioner, all detailed analysis in the remainder of this work concentrates on the multigrid V-cycle.

One should note that spending 72%/76% of the solution time within the multigrid preconditioner is already low, particularly taking into account that the outer conjugate-gradient solver also performs its operator evaluations (1 matrix-vector multiplication per iteration) in an efficient matrix-free fashion. The low costs are related to the usage of single-precision floating-point numbers and to the low number of pre-/postsmoothing steps, which results in a total of six to seven matrix-vector multiplications per level and iteration.

Tables 4 and 5 show the number of iterations and the time to solution for the octant, and the shell test cases run serially with local smoothing and global coarsening as well as with the polynomial degrees \(p=1\) and \(p=4\).

It is well visible that local smoothing has to perform at least as many iterations as global coarsening, with the difference in iterations limited to 1 in the examples considered. This difference is due to the additional smoothening applied to certain cells in the course of global coarsening. Note that global coarsening also benefits from the simple setup of a smooth solution with artificial refinement, see Table 7 for a somewhat more realistic test case. This comes with a higher serial workload for global coarsening, which is also visible in the times of a single V-cycle (not shown). Nevertheless, the lower number of iterations leads to a smaller time to solution in the case of global coarsening in some instances. Generally, the global-coarsening V-cycle is relatively more expensive in the case of the shell simulation with linear elements: This is not surprising due to the higher number of cells with hanging-node constraints (see also Figure 7) in the shell case and the higher overhead of linear elements for application of hanging-node constraints, as analyzed in Reference [81].

Figure 11 shows the distribution of times spent on each multigrid level and in each multigrid stage for the octant case with \(L=8\) and \(p=4\). While the overall runtimes are similar for local smoothing and global coarsening, distinct (and expected) differences are visible for the individual ingredients: The restriction and prolongation steps take about the same time in both cases. The presmoothing, residual, and postsmoothing steps are slightly more expensive in the global-coarsening case, which is related to the observation that the evaluation of the level operator \(\boldsymbol {A}^{(l)}\) (2/1/3 times) is the dominating factor. For local smoothing, the computation of \(\boldsymbol {A}_{ES}^{(l)}\boldsymbol {x}_S^{(l)}\) is realized as a side product of the application of \(\boldsymbol {A}^{(l)}\), and hence limiting its impact on the residual step. The only visible additional cost of local smoothing is \(\boldsymbol {A}_{SE}^{(l)}\boldsymbol {x}_E^{(l)}\) for the modification of the right-hand-side vector for postsmoothing to incorporate inhomogenous Dirichlet boundary conditions (edge step). However, its evaluation is less expensive than the application of \(\boldsymbol {A}^{(l)}\), since only DoFs in proximity to the interface are updated.

We report our findings of simulations with 192 processes on four compute nodes. Tables 4 and 5 confirm the significance of a good workload balance for the time to solution, as indicated by Tables 1 and 2. The number of processes does not influence the number of iterations due to the chosen smoother. Speedups—compared to local smoothing—are reached: up to 2.3/1.7 for the octant (\(p=1\)/\(p=4\)) and up to 2.1/1.6 for the shell case if the different iteration numbers are considered. Normalized per solver iteration, the advantage of global coarsening in these four cases is 2.0, 1.3, 1.7, and 1.6, respectively. Figure 12 illustrates this behavior by showing the distribution of the minimum/maximum/average times spent on each multigrid level and each multigrid stage for the octant case with \(L=8\). In the case of global coarsening, it is well visible that the load is equally distributed and the time spent on the levels is significantly reduced level by level for the higher levels. For the finest ones, local smoothing shows a completely different picture. On the finest level, there are processes with hardly any work, but nevertheless the average work is close to the maximum value, indicating that the load is well balanced among the processes with work. However, on the second finest level, a significant workload imbalance—by a factor of 2.8—is visible also among the processes with work. This leads to the situation that the maximum time spent on the second finest level is just slightly less than the one on the finest level, contradicting our expectation of a geometric series and leading to the observed increase in the total runtime.

Figures 13 and 14 show results of scaling experiments starting with 1 compute node (48 processes) up to 3,072 nodes (147,456 processes). Besides the times of a single V-cycle, we plot the normalized throughput (DoFs per process and time per iteration) against the time per iteration. The plot can be read in the following way: Far from the scaling limit (right top corner of the plot), simulations take longer but are more efficient (parallel efficiency of one); at the scaling limit (left bottom corner of the plots), simulations reach shorter times at the cost of lower efficiencies. Here, a horizontal behavior from the right to the left corresponds to ideal strong scaling. If different lines coincide, then we observe ideal weak scaling. Based on this plot, one can judge the resource utilization and parallel efficiency of a simulation if a certain time to solution should be obtained.

The throughput for the polynomial degree \(p=4\) is higher by a factor of 3 than for \(p=1\). This is expected and related to the used matrix-free algorithms and their node-level performance, which improves with the polynomial order [66]. Just as in the moderately parallel case (see Section 5.3), we can observe better timings in the case of global coarsening for a large range of configurations (maximum speedup: octant 1.9/1.4 for \(p=1\)/\(p=4\), shell: 2.4/2.4). The number of iterations of local smoothing is 4 for both cases and all refinement numbers. For global coarsening, it is 4 for the sphere case and decreases from 4 to 2 with increasing number of refinements in the case of quadrant so that the actual speedups reported above are even higher. The high speedup numbers of global coarsening in the shell simulation case are particularly related to its high workload and vertical efficiency, as shown in Table 5.

The normalized plots give additional insights. Apart from the obvious observation that the minimal time to solution increases with increasing number of refinements (left bottom corner of the plots in Figures 13 and 14) and global coarsening starts with higher throughputs, one can also see that the decrease in parallel efficiency is more moderate in the case of global coarsening. For the example of the octant case with \(p=1\)/\(L=10\), one can increase the number of processes by a factor of 16 and still obtain a throughput per process that is higher than the one in the case of single-node computations of local smoothing, which normally shows a kink in efficiency at early stages (particularly visible in the shell case, which matches the findings made in Reference [28]). A further observation is that the lines for global coarsening overlap far from the scaling limit, i.e., the throughput is independent of the number of processes and the number of refinements. This is not the case for local smoothing, where the throughput deteriorates with the number of levels, indicating load-balance problems. The simulations with \(p=4\) show similar trends, but the lines are not as smooth, possibly due to the decreased granularity for higher orders.

With Equation (9), one can approximate a latency \(\Delta t_{\text{Latency}}^{\texttt {V-Cycle}}=3.6\) ms for \(L=9\) and \(\Delta t_{\text{Latency}}^{\texttt {vmult}}=50~\mu\)s. This value is visualized with vertical lines in the diagrams. One can see that the approximation indeed matches well for \(p=1\), indicating that the solvers are limited by the latency in these ranges. Due to the similarity between the numbers of communication steps in the case of local smoothing and global coarsening, the minimal times reached are similar. For \(p=4\), the discrepancy between the approximated limit and the measurements is higher; the minimal times reached by global coarsening are smaller, indicating that the latency limit has not been reached yet and load imbalances still have a noticeable effect. For less refinements, e.g., \(L=7\), better correspondence can be observed. Apart from these observations, we emphasize that, in practice, one would run at higher efficiencies (left-top corner in the normalized plot), for which both \(p=1\) and \(p=4\) are dominated by load-imbalance effects.

Section 5.3 indicates that local smoothing with first-child policy might suffer from deteriorated reduction rates of the maximum number of cells on each level; in particular, there might be processes without any cells, i.e., any work, increasing the critical path, although the vertical efficiency is optimal. In this section, we consider the first-child policy as an alternative for partitioning of the levels for global coarsening.

Figure 15 presents the timings of large-scale octant simulations for (1) local smoothing, (2) global coarsening with default partitioning, and (3) global coarsening with first-child policy for \(p=1\)/\(p=4\). The timings of the latter show similar trends as global coarsening with default partitioning and are lower than the ones of local smoothing. To explain this counterintuitive observation, Figure 15 provides the maximum number of cells on each multigrid level of the three approaches for \(L=11\) on 256 nodes. Since global coarsening with first-child policy does not perform any repartitioning, better reduction trends as in the local-smoothing case on all levels cannot be expected; however, one can observe that, for the local section of the refinement tree, the number of cells on the finest levels is reduced nearly as well as in the case of the default partitioner, i.e., the parallel workload and the parallel efficiency are not much worse; nonetheless, with higher vertical efficiency. This is not possible for the lower levels, but the behavior of the finest levels dominates the overall trends, since they take the largest time of the computation. One should note that local smoothing has access to the same cells but simply skips them during smoothing of a given level, missing the opportunity to reduce the problem size locally. In our experiments, it is not clear whether an optimal reduction of cells is crucial for all configurations: For \(p=1\), global coarsening with default partitioning is faster for most configurations (up to 20%) and, for \(p=4\), global coarsening with first-child policy is faster (up to 10%). We could trace this difference back to the different costs of the transfer, where, for \(p=4\), a reduced data transfer (both within the compute node and across the network) and, for \(p=1\), a better load balance is beneficial.

In the shell case (Figure 16), we observe that both global-coarsening partitioning strategies result in very similar reduction rates, which can be traced back to the fact that both strategies lead to comparable partitionings of the levels (see also Table 5, which shows a very high vertical efficiency of the default partitioning) and—as a consequence—to very similar solution times (\(\pm\)10%). Local smoothing only reduces the maximum number of cells per level optimally once all locally refined cells have been processed and the levels that had been constructed via global refinement have been reached. This case stresses the issue of load imbalances related to reduction rates significantly differing between processes.