Scaling to the stars – a linearly scaling elliptic solver for p-multigrid

Abstract

High-order methods gain increased attention in computational fluid dynamics. However, due to the time step restrictions arising from the semi-implicit time stepping for the incompressible case, the potential advantage of these methods depends critically on efficient elliptic solvers. Due to the operation counts of operators scaling with with the polynomial degree p times the number of degrees of freedom nDOF, the runtime of the best available multigrid solvers scales with O(p · nDOF). This scaling with p significantly lowers the applicability of high-order methods to high orders. While the operators for residual evaluation can be linearized when using static condensation, Schwarz-type smoothers require their inverses on fixed subdomains. No explicit inverse is known in the condensed case and matrix-matrix multiplications scale with p · nDOF. This paper derives a matrix-free explicit inverse for the static condensed operator in a cuboidal subdomain. It scales with p 3 per element, i.e. nDOF globally, and allows for a linearly scaling additive Schwarz smoother, yielding a p-multigrid cycle with an operation count of O(nDOF). The resulting solver uses fewer than four iterations for all polynomial degrees to reduce the residual by ten orders and has a runtime scaling linearly with nDOF for polynomial degrees at least up to 48. Furthermore the runtime is less than one microsecond per unknown over wide parameter ranges when using one core of a CPU, leading to time-stepping for the incompressible Navier-Stokes equations using as much time for explicitly treated convection terms as for the elliptic solvers. * Corresponding author: Immo.Huismann@tu-dresden.de arXiv:1808.03595v1 [math.NA] 10 Aug 2018 to h-multigrid. In both cases, smoothers are required and overlapping Schwarz smoothers are a good option [18,29,41]. These decompose the grid into overlapping blocks on which the exact inverse is applied, producing exceptional smoothing rates. But these methods not only require the super-linearly scaling operator for residual evaluation, but a super-linearly scaling solver as well. This super-linear scaling prevents the respective multigrid techniques from achieving their full potential and prohibits the usage of large polynomial degrees.The goal of this paper is to devise a linearly scaling multigrid cycle for three-dimensional structured Cartesian grids, where residual evaluation, smoothing, restriction, and prolongation all scale linearly with the number degrees of freedom, without a further factor of p present in other methods. To this end the residual evaluation derived in [22] for the static condensed case is combined with the p-multigrid block smoother technique proposed in [18,17], leaving only the inverse on a 2 3 element block as non-matrix-free, super-linearly scaling operator. This operator is embedded in the full system and then factorized and rearranged, attaining a matrix-free inverse that is then factorized to linear complexity. The inverse leads to a linearly scaling multigrid cycle, which is confirmed in run...

Linearizing the hybridizable discontinuous Galerkin method: A linearly scaling operator

“…However, the pre-and postprocessing occur, by definition, once per solution process, whereas the iteration process is reiterated tens if not hundreds of times. For the continuous method, the iteration, not pre-and postprocessing are the most costly component for polynomial degrees up to p > 48 [24], and these lie far outside the range of currently employed ones. Therefore, the solver can be described as scaling linearly for all relevant polynomial degrees.…”

“…These two operations scale with O p 4 n e and are only evaluated once during the solution process, such that they not impede the linear scaling. They are not expected to dominate the runtime behaviour until at least p > 48 [24], allowing to disregard the issue here.…”

“…Assuming a constant iteration count, only tensor-product preconditioners working on the faces remain a possibility. In general, a multigrid approach leads to very efficient solvers, as done for the continuous discretization in [14,24,17]. There, Schwarz-type preconditioners using tensor-products were employed to generate an efficient smoother with linear scaling.…”

“…For a solver to scale linearly with the number of unknowns when increasing the polynomial degree, a linearly scaling operator is mandatory. While these linearly scaling operators are available for the continuous discretization [22,24,17], they have so far only been postulated for the discontinuous one [29], and, to the knowledge of the authors, no implementation thereof exists.…”

“…In previous work [22,24], the present authors have considered a continuous spectralelement discretization, developing a linearly scaling operator which, in conjunction with multigrid techniques, allows for a linearly scaling elliptic solver. The resulting runtime scales with O(n DOF ) independent of the polynomial degree.…”