PaperlayerView exact revision on arXiv

Source-linked AI summary

Adaptive Multilevel Discontinuous Galerkin Methods on GPUs

Michał Wichrowski

arXiv:2609.04534v1math.NAcs.NA

TL;DR

Adaptive refinement makes uniform GPU face evaluation and local multigrid treatment difficult because refinement edges require special handling. Shadow cells convert those interfaces into matching faces, and experiments on an NVIDIA A100 report efficient execution with constant iteration count at four local refinement levels.

  • Problem

    Adaptive meshes require separate edge kernels, constrained transfers, and non-matching face treatments that introduce irregular addressing and divergent execution paths.

    Source: §1 · p. 2

  • Method

    Shadow cells represent coarse polynomials on virtual fine children, making refinement-edge faces matching; primal–dual pairing and frozen shadows integrate them into Krylov iteration and local multigrid.

    Source: §1 · p. 2, §7 · p. 22

  • Results

    On an NVIDIA A100, the face exchange and blocked layout reach 81% of peak device memory bandwidth; at four local refinement levels, throughput is 71.5% of uniform-mesh throughput with constant iteration count.

    Source: §7 · p. 22

  • Takeaways & Limitations

    Shadow evaluation reproduces standard non-matching face integrals and quadrature exactly, preserving the standard discrete problem and its approximation properties.

    Source: §7 · p. 22

    The implementation and measurements assume 2:1 balance; deeper level jumps are proposed as an extension but have not been tested.

    Source: §3.2 · p. 8

Abstract

from arXiv · show

I present a matrix-free symmetric interior penalty discontinuous Galerkin method for adaptively refined Cartesian meshes on GPUs. At non-matching interfaces, auxiliary \emph{shadow cells} represent the adjacent coarse polynomial on the fine level. This approach converts non-matching interfaces into matching faces, permitting uniform face evaluation throughout the mesh. I prove that this construction is equivalent to the standard non-matching formulation. Furthermore, this representation yields a local geometric multigrid method in which frozen shadows provide inter-level boundary data and carry residual contributions to coarser levels. A primal--dual pairing eliminates shadow assembly from the Krylov iteration. Numerical experiments with cubic elements ($p=3$) on an NVIDIA A100 show stable multigrid convergence under increasing refinement depth and efficient GPU execution.

1 Introduction

The paper makes adaptive DG computation more uniform by replacing non-matching interfaces with matching-face exchanges and shadow cells. It combines this representation with a local multigrid method and primal–dual pairing, and reports efficient execution under adaptive refinement.

1 Introduction

  • Adaptive meshes complicate uniform GPU execution: existing local smoothing and non-matching face treatments use separate edge kernels, constrained transfers, subface quadrature, and irregular addressing.Local smoothing follows the refinement hierarchy rather than globally refined meshes, but its level operators must account for refinement-edge coupling.Source: §1 · pp. 1–2
  • Structured face exchanges enumerate faces arithmetically and process each interface once, using direct neighbor loads and race-free two-sided writes without atomics or coloring.Indirection remains only on the skeleton between macro-blocks.Source: §1 · p. 2
  • Under 2:1 balance, shadow cells store prolonged coarse polynomials as virtual fine children, making every face matching and reproducing standard non-matching integrals and quadrature exactly.Shadow contributions return to the parent by transposed prolongation, without introducing an approximation.Source: §1 · p. 2
  • A primal–dual pairing lets conjugate gradients use split active and shadow data while eliminating explicit shadow assembly from every Krylov operator application.Convergence is monitored through the duality pairing rather than the standard Euclidean residual norm.Source: §1 · pp. 2–3
  • In local geometric multigrid, frozen shadows supply inter-level boundary data and carry outward-flux residual contributions to coarser levels.Shadow values are refreshed once per level visit; subsequent smoothing sweeps and residual evaluations use the uniform face loop.Source: §1 · p. 3
  • At four levels of local refinement, the operator retains 71.5% of uniform-mesh throughput while iteration count stays constant; the A100 face kernel reaches 81% of peak memory bandwidth.The approach uses surface-order shadow storage and one shadow refresh and assembly per level visit, not per operator application.Source: §1 · p. 3
  • The solver is implemented in Triton, whose dense tile model matches the blocked cell-wise layout and supports compiler-managed GPU memory and tensor-core mapping.The same kernel source can target NVIDIA and AMD GPUs.Source: §1 · p. 3
  • The blocked cell-wise storage layout is treated as a performance parameter, decoupling persistent data arrangement from mesh topology for DG elements.The paper benchmarks blocked layouts for the interior penalty operator.Source: §1 · p. 3

2 Discontinuous Galerkin Discretization

The paper formulates a matrix-free SIPG discretization for Poisson problems on adaptively refined Cartesian meshes, then separates volume and face work to accommodate shadow cells. Its Cartesian setting and blocked storage layout are chosen with data movement in mind.

2 Discontinuous Galerkin Discretization

  • The model is Poisson on a bounded 2D or 3D domain with homogeneous Dirichlet conditions, using SIPG on axis-aligned Cartesian cells and tensor-product polynomials Q_p.The method and kernels are formulated for general degree p, while implementation measurements use cubic elements, p = 3.Source: §2 · p. 3
  • The matrix-free operator recomputes volume and face integrals by sum-factorization, keeping the two kernel families separate under adaptive refinement.Volume kernels run over cells with degrees of freedom, while face sweeps include auxiliary shadow cells so the face loop can remain uniform.Source: §2 · p. 4
  • The Cartesian-grid restriction avoids per-point geometry traffic that can erode bandwidth performance on deformed meshes; coupling complex-geometry layers to these kernels is outside the stated scope.On Cartesian cells, kernels stream no geometry data and face integrals factorize into precomputed one-dimensional matrices.Source: §2 · p. 4
  • Blocked storage arranges field data as [B, nx, ny, nz, ne], tuning the balance between coalesced volumetric access and compact face-trace extraction.The optimal block size depends on the mix of volume and face work, which varies with polynomial degree.Source: §2 · pp. 4–5
  • On adaptive meshes, a primal–dual formulation removes shadow assembly from the Krylov loop while the active-cell residual needs no assembly.The residual on active cells requires no assembly on a uniform DG mesh either.Source: §2 · p. 5

3 Evaluation of Face Integrals

The method organizes face evaluation around structured, axis-wise sweeps and changes the basis so face data can be read from two adjacent layers. Shadow cells turn refinement-edge interfaces into matching faces, preserving the standard SIPG contribution while enabling uniform face evaluation.

3 Evaluation of Face Integrals

  • Structured, axis-wise face sweeps use direct neighbor loads and race-free two-sided writes, evaluating each face once without atomics or coloring.Sequential axis sweeps resolve accumulation onto shared layers; only the macro-block skeleton retains indirect neighbor lookup.Source: §3 · p. 5

3.1 Tensor-Product Evaluation and the Hermite Basis

  • Precomputed one-dimensional face-coupling matrices reduce runtime face integration to dense contractions, with four matrices shared across faces and axes.The complete face-integration data comprises 4(p + 1)^2 numbers per level, and no face quadrature is performed at runtime.Source: §3 · pp. 5–6

3.1.1 The Hermite-Type Basis and Two-Layer Face Access

  • The Hermite-type basis makes traces and normal derivatives available in two face-adjacent layers, reducing face traffic by a factor of (p + 1)/2.The change of basis leaves the discrete space and discretization unchanged.Source: §3.1.1 · p. 7

3.2 Shadow Cells for Non-Matching Faces

  • Shadow cells replace non-matching refinement-edge faces with matching fine–shadow faces, allowing the face loop to process them uniformly.Shadows carry prolonged coarse-cell data and have no independent degrees of freedom.Source: §3.2 · p. 8
  • The implementation assumes 2:1 balance; deeper level jumps would require nested shadows and have not been tested.The authors state that they see no obstruction to extending the construction to deeper jumps.Source: §3.2 · p. 8
  • Shadow refresh prolongs parent data before face sweeps, while shadow restriction returns accumulated contributions to parents when an assembled residual is required.The Krylov iteration can omit both transfers because its split representation preserves shadow consistency; multigrid performs assembly during level restriction.Source: §3.2 · pp. 8–10
  • With the fine-level penalty on each subface, matching fine–shadow forms equal the standard non-matching SIPG contribution, including computed quadrature values.The fine-level face rule maps to the shifted and scaled rule used by a dedicated non-matching kernel.Source: §3.2 · p. 10

4 Duality Pairings on Shadowed Storage

Shadowed storage represents adaptive DG residuals as split dual data and consistent primal vectors, allowing their active-mesh pairings to be computed without shadow assembly. The resulting shadowed conjugate-gradient method is equivalent to preconditioned CG on the assembled active-mesh system when the stated representation conditions hold.

Extended storage

  • Consistent primal vectors prolong parent values to shadows, while split dual residuals store parent-test-function contributions on shadow blocks.The right-hand side has local load integrals on active cells and zeros on shadows.Source: §4 · pp. 11–12

Operator identity

  • Shadow assembly of the split shadowed operator output recovers the exact active-mesh residual for a consistent input vector.The active cells hold volume and matching-face terms, while shadow blocks hold refinement-edge fluxes for their parents.Source: §4 · p. 12

Duality pairing

  • A flat dot product over active and shadow blocks gives the exact active-mesh energy pairing without shadow assembly.This pairs a split dual vector with a consistent primal vector.Source: §4 · pp. 11–12

Krylov iteration

  • Shadowed conjugate gradients pair each split dual residual or operator output with a consistent primal vector, so the iteration omits assembly after operator applications.Vector updates preserve the two representations, and convergence is monitored through the duality pairing rather than the split residual’s Euclidean norm.Source: §4 · pp. 12–13

Equivalence

  • When the preconditioner satisfies the stated consistency and residual-dependence conditions, shadowed CG produces the same iterates as preconditioned CG on the assembled active-mesh system.The local multigrid method satisfies the residual condition by transferring shadow contributions during level restriction.Source: §4 · pp. 12–13

5 Local Multigrid with Adaptively Refined Meshes

The local geometric multigrid V-cycle smooths only cells on each level, with shadows representing refinement-edge coupling in the face loop. Its level transfers and frozen shadow data handle inter-level coupling within this local hierarchy.

5.1 The Local Multigrid Method

  • Local multigrid restricts smoothing to level cells, making V-cycle work proportional to the cells summed over adaptive levels rather than a globally refined hierarchy.Level right-hand sides carry active residuals on active cells and vanish on refined cells; shadow assembly adds edge contributions to coarser parents during restriction.Source: §5.1 · p. 13
  • Level prolongation and its transpose provide Galerkin transfers between nested broken spaces without constraints or interface weighting.On cycle exit, level prolongation also refreshes shadows to restore consistency.Source: §5.1 · p. 13

5.2 Smoother

  • The smoother uses cell-wise block-Jacobi with a tensor-structured block whose inverse is applied by fast diagonalization, requiring no inter-cell communication.The block discards off-cell face couplings; communication within the V-cycle occurs through level-operator face terms and inter-level transfers.Source: §5.2 · p. 14
  • Three varying-weight Richardson updates give better worst-case damping than three block-Jacobi steps with the same optimal constant weight, without changing operator or block-inverse application counts.The factorized Chebyshev form adds no auxiliary search-direction vector or additional full-vector stream.Source: §5.2 · p. 15

5.3 Treatment of the Refinement Edge

  • During a level’s smoothing block, shadows hold fixed coarse-side correction data, beginning at zero before coarse correction and receiving prolonged correction for post-smoothing.These shadows carry the traces into the uniform face loop without masking or operator splitting.Source: §5.3 · p. 15
  • After smoothing, fine-shadow face kernels accumulate refinement-edge fluxes, which shadow assembly transfers to parent residuals alongside level restriction.These are the non-matching face contributions of the active-mesh operator; in multigrid, shadow assembly is fused with restriction.Source: §5.3 · p. 15

6 Numerical Results

The experiments assess solver convergence, GPU kernel performance, and throughput as local refinement depth increases, using degree-three elements on an A100. They test globally and locally refined Cartesian meshes under the stated solver and mesh settings.

6.1 Experimental Settings

  • The experiments use p=3 Poisson problems on 2:1-balanced Cartesian meshes, with a three-update Chebyshev-accelerated block-Jacobi smoother and CG stopping at a 10^-10 preconditioned-residual reduction.Adaptive cells are marked in a radial shell; the coarsest level is solved exactly.Source: §6.1 · p. 16

6.2 Penalty Parameter and Solver Convergence

  • At p=3, the solver uses C=2 (σF=24/hF); CG takes 10–12 iterations, with two added at the first refinement and none as depth increases.The fine–shadow face penalty uses the fine-cell size, doubling the penalty relative to adjacent coarse-side faces; the authors suspect this stiffening explains the two added iterations.Source: §6.2 · pp. 16–17

6.3 The Face Kernel

  • The coalesced face kernel reaches 22.67 GDoF/s and 1650.8 GB/s, or 81% of peak bandwidth, versus 5.60 GDoF/s for face-major ordering.The coalesced evaluation retains 95.9% of its arithmetic-free stream reference; the compared variants and references are summarized in Table 1.Source: §6.3 · p. 18, Table 1 · p. 19
  • At B=32, face throughput peaks at 22.38 GDoF/s, a 49.4% improvement over B=1, while volume throughput retains 96.9% of its maximum.Across measured block sizes, volume throughput varies from 86.48 to 96.12 GDoF/s; Figure 3 compares both kernels.Source: §6.3 · p. 18, Figure 3 · p. 19

6.4 Roofline Analysis

  • The face kernel is memory-bound, reaching 1650.8 GB/s, or 81% of peak memory bandwidth; the volume kernel reaches 8.14 TFLOP/s at 5.00 FLOP/byte.The reported FLOP rates are upper bounds because hardware-counter operation counts include arithmetic the operator does not use; measured bandwidth is unaffected.Source: §6.4 · p. 20, Figure 4 · p. 20

6.5 Operator and Solve Throughput

  • On globally refined meshes, CG takes 10–11 iterations and solve throughput rises from 3.3 to 89.6 MDoF/s.The reported comparison solver reaches 61.86 MDoF/s at p=3, but uses 10^-8 relative tolerance and a vertex-patch smoother, versus 10^-10 and block Jacobi here.Source: §6.5 · p. 21, Table 2 · p. 21
  • At four local refinement levels, the operator retains 71.5% and the smoother 54.9% of their uniform-mesh rates; solve throughput falls from 74.4 to 39.5 MDoF/s.Across 262,144 to 3,936,304 active cells, the operator sustains 11.12–15.55 GDoF/s; the iteration count rises from 10 to 12 after the first refinement and then stays constant.Source: §6.5 · p. 21, Table 3 · p. 21

7 Conclusion

Shadow cells turn refinement interfaces into uniform faces while preserving the standard non-matching formulation. Experiments show efficient GPU execution and constant multigrid iteration counts across the tested refinement depths.

7 Conclusion

  • Projecting the coarse-cell polynomial onto virtual fine-level children removes non-matching interfaces and permits a uniform, dimensionally split face-evaluation kernel.This avoids the irregular addressing and warp divergence associated with standard subface integration.Source: §7 · p. 22
  • Shadow evaluation reproduces the standard non-matching formulation's face integrals and quadrature exactly, preserving its discrete problem and approximation properties.Source: §7 · p. 22
  • A blocked layout and Hermite-type basis limit inter-cell dependencies to immediate face layers, while primal–dual pairing eliminates shadow assembly during Krylov operator application.Source: §7 · p. 22
  • In local geometric multigrid, shadows provide inter-level boundary data and transfer coarse-grid residual corrections.Source: §7 · p. 22
  • 81% of peak device memory bandwidth, within 4.1% of an arithmetic-free reference using the same access pattern; at four refinement levels, throughput is 71.5% of uniform-mesh throughput and multigrid iteration count stays constant.The experiments were performed on an NVIDIA A100.Source: §7 · p. 22
  • The demonstrated scope is a Cartesian background mesh; combining the method with an unfitted or shifted-boundary layer for complex domains is left to future work.Source: §7 · p. 22
Loading 2609.04534v1…