PaperlayerView exact revision on arXiv

Source-linked AI summary

Nodal Discontinuous Galerkin Methods on Graphics Processors

Andreas Klöckner, Tim Warburton, Jeffrey Bridge, Jan S. Hesthaven

arXiv:0901.1024v3math.NA

TL;DR

The paper addresses how to exploit DG’s locality and high-order arithmetic intensity on GPUs while handling GPU-specific partitioning and memory-layout challenges. It develops CUDA-based strategies for nodal DG and evaluates them on a 3D unstructured Maxwell solver, reporting speedups approaching two orders of magnitude and practical guidance for implementation.

  • Problem

    GPU DG implementations must determine effective work partitioning and replace serial processors’ tuned linear-algebra and communication primitives, which may be unavailable or unsuitable on GPUs.

  • Method

    The paper adapts and combines strategies for mapping nodal DG operators and their subtasks onto Nvidia CUDA GPUs, including optimized layouts and memory-access schedules.

  • Results

    DG computations on GPUs achieve speedups just short of two orders of magnitude, while high-order methods reach double-digit percentages of published theoretical peak performance.

  • Takeaways & Limitations

    GPU-oriented DG strategies can substantially increase the problem size and complexity affordable on a given hardware budget.

  • Takeaways & Limitations

    The flux gather remains the DG stage least suited to GPU execution because it is branch-intensive, memory-erratic, and register-demanding.

Abstract

from arXiv · show

Discontinuous Galerkin (DG) methods for the numerical solution of partial differential equations have enjoyed considerable success because they are both flexible and robust: They allow arbitrary unstructured geometries and easy control of accuracy without compromising simulation stability. Lately, another property of DG has been growing in importance: The majority of a DG operator is applied in an element-local way, with weak penalty-based element-to-element coupling. The resulting locality in memory access is one of the factors that enables DG to run on off-the-shelf, massively parallel graphics processors (GPUs). In addition, DG's high-order nature lets it require fewer data points per represented wavelength and hence fewer memory accesses, in exchange for higher arithmetic intensity. Both of these factors work significantly in favor of a GPU implementation of DG. Using a single US$400 Nvidia GTX 280 GPU, we accelerate a solver for Maxwell's equations on a general 3D unstructured grid by a factor of 40 to 60 relative to a serial computation on a current-generation CPU. In many cases, our algorithms exhibit full use of the device's available memory bandwidth. Example computations achieve and surpass 200 gigaflops/s of net application-level floating point work. In this article, we describe and derive the techniques used to reach this level of performance. In addition, we present comprehensive data on the accuracy and runtime behavior of the method.

1. Introduction

DG combines high-order approximation with element decomposition and weak coupling, making it flexible for unstructured geometries and attractive for GPU execution. The paper investigates mapping DG solvers for linear hyperbolic systems onto Nvidia CUDA GPUs, emphasizing high-order arithmetic intensity and practical performance.

  • 1. Introduction: DG combines high-order approximation with Finite-Volume-like surface Riemann solvers between computational elements.This hybrid design combines advantages associated with Finite-Volume and Spectral Element methods.
  • 1. Introduction: GPU execution is motivated by memory bandwidth and latency becoming dominant performance factors in modern computing.Graphics processors offer a different balance between memory access and computational throughput than conventional CPUs.
  • 1. Introduction: The paper addresses partitioning GPU computational work and replacing serial processors’ tuned linear-algebra and communication primitives.These are identified as two central challenges for GPU DG implementations.
  • 1. Introduction: More than an order-of-magnitude speedup is achieved on a single real-world consumer graphics processor versus a CPU implementation of the same method.The paper presents this as one of the first such general finite-element-based solver results.
  • 1. Introduction: High-order methods increase arithmetic intensity, shifting computation toward GPU-favorable compute bandwidth rather than memory bandwidth.The paper attributes a sizable part of the speedup to this property.
  • 1. Introduction: The study focuses on linear hyperbolic conservation laws and tetrahedral elements, leaving nonlinear, elliptic, and parabolic problems for future work.Tetrahedra are selected partly because their DG formulation is comparatively arithmetically intense and memory-efficient.

2. Overview of the Discontinuous Galerkin Method

The nodal DG formulation represents elementwise polynomial solutions at interpolation nodes and expresses volume and surface contributions through local matrices and face fluxes. Its element-local structure permits explicit timestepping and supports a decomposition into GPU-suitable subtasks.

  • 2. Overview of the Discontinuous Galerkin Method: The method approximates a hyperbolic conservation-law solution with local polynomials of maximum degree N on tetrahedral elements.Lagrange basis functions and interpolation nodes define Np local degrees of freedom.
  • 2. Overview of the Discontinuous Galerkin Method: The strong-DG formulation uses volume and surface terms, with numerical fluxes enforcing coupling across element faces.The supplied passages introduce boundary conditions, weak form, and strong-DG reformulation.
  • 2. Overview of the Discontinuous Galerkin Method: The lifting matrix combines face mass application, embedding facial values into a volume vector, and inverse volume-mass application.It converts facial contributions into volume contributions.
  • Overview of the Discontinuous Galerkin Method: Left multiplication by the inverse mass matrix is elementwise, enabling explicit Runge-Kutta timestepping without global communication.This property distinguishes DG from other finite-element methods and simplifies GPU implementation.
  • Implementing DG: DG decomposes naturally into four stages because its operator contains additive volume and surface contributions, with the surface term further split into gather-related work.The decomposition is illustrated in Figure 2.
  • Implementing DG: Nodal DG stores values at interpolation nodes, allowing facial values to be obtained by selecting face nodes from the volume field.This differs from modal DG, where degrees of freedom are expansion coefficients.
  • Implementing DG: Most DG stages are element-local and can often be represented efficiently by dense matrix-vector products on each element.The element-local operations avoid dependence on neighboring elements.

3. The CUDA Parallel Computation Model

CUDA exposes parallelism through threads, warps, blocks, and grids, with shared memory and registers supporting cooperation and data reuse. GPU performance depends on hiding global-memory latency, avoiding bank conflicts, and matching access patterns to alignment and bandwidth requirements.

  • 3. The CUDA Parallel Computation Model: Nvidia GPU hardware provides two levels of parallelism inherited from graphics workloads: pixel-like and primitive-like processing.The architecture targets large numbers of independently processed geometric work units.
  • 3. The CUDA Parallel Computation Model: A multiprocessor executes groups of eight functional units under one instruction decoder, including fused floating-point multiply-add operations.The supplied hardware description emphasizes the execution structure relevant to throughput.
  • 3. The CUDA Parallel Computation Model: Thread blocks share execution hardware, barriers, memory fences, and 16KiB of banked shared memory.Shared-memory banking allows simultaneous access when threads address distinct banks.
  • 3. The CUDA Parallel Computation Model: Grid blocks cannot communicate or rely on ordering during execution; completion of a grid submission provides synchronization.This constrains how multi-stage computations must be scheduled.
  • 3. The CUDA Parallel Computation Model: Global-memory latency reaches several hundred clock cycles, so the GPU schedules other ready warps to hide it.Registers and shared-memory usage limit how many threads remain available for scheduling.
  • 3. The CUDA Parallel Computation Model: For 32-bit global accesses, highest bandwidth requires warp threads to access locations aligned and grouped according to 16-thread boundaries.Access arrangement is therefore a central part of GPU data-layout design.
  • 3. The CUDA Parallel Computation Model: GPU floating-point capacity exceeds its already larger global-memory bandwidth by another order of magnitude, requiring both computational and memory pipelines to remain active.The paper frames algorithm design as keeping both resource paths flowing efficiently.

4. DG on the GPU: Design

The GPU design maps DG’s natural granularities and subtasks onto CUDA through coordinated choices of computation layout, data layout, fetch schedule, and thread-block decomposition. Because subtasks have different output sizes and on-chip-memory needs, separate optimized layouts are preferred despite added fetches.

  • DG on the GPU: Design: The central GPU-mapping decisions concern work partitioning, data layout, and fetch scheduling.These choices determine how DG work units are assigned and how data moves between global and on-chip storage.
  • DG on the GPU: Design: One thread per output aligns computation and storage, while the fetch schedule controls reuse of data in registers or shared memory.Post-computation permutations can require additional shared memory.
  • DG on the GPU: Design: DG exposes granularities from element, face, and system sizes, including Np, Nfp, Nf, and n.The number of elements K also affects work partitioning but is less important in this discussion.
  • DG on the GPU: Design: Moderate-order tetrahedral DOF counts often misalign with GPU batches of 16 and 32, making simple padding waste memory and processing.Figure 3 presents the relevant granularities and alignment issue.
  • DG on the GPU: Design: Different thread-block work factors trade parallelism and shared-memory use against register reuse and sequential processing.The parameters wp, wi, and ws describe work processed in parallel, inline, and sequentially.
  • DG on the GPU: Design: The system size n can influence data and computation layouts, but packed field storage is undesirable for several reasons.The supplied passage introduces this layout choice without detailing all reasons.
  • DG on the GPU: Design: Splitting the operator into subtasks is preferable because each stage can use on-chip memory differently and gather and lift produce different output sizes.Gather outputs NfpNf values per element, whereas lift outputs Np, so one shared layout is suboptimal.

5. DG on the GPU: Implementation

The implementation decomposes DG work into GPU-suitable stages and uses custom layouts, shared memory, texture caching, and partitioning to manage reuse, alignment, conflicts, and scattered access. These choices improve locality while accepting bounded fetch redundancy and hardware-specific constraints.

  • Flux lifting: Flux lifting uses element-local matrix-matrix multiplication followed by elementwise scaling, with field and matrix data arranged for reuse.
  • Flux lifting: Custom GPU algorithms replace vendor BLAS because its performance and alignment requirements are unsuitable for the targeted matrix operations.
  • Flux lifting: Column-major matrix storage improves locality when threads assigned to matrix rows load values from successive columns through texture units.
  • Flux lifting: Microblocking can create double-broadcast bank conflicts at element boundaries, but separating conflict-prone and conflict-free half-warps mitigates their impact.
  • Flux extraction: Partitioning is constrained by shared-memory buffering capacity and the number of block-external faces, which require additional face descriptors.
  • Flux extraction: Flux extraction is the GPU’s least-suited DG stage because its data-driven branches, erratic memory access, and growing register demand hinder execution.
  • Element-local differentiation: Conflict-free shared-memory segmentation introduces fetch redundancy because segments may reload field values from neighboring elements; adjacent blocks and L2 caching may reduce the bandwidth cost.

6. Experimental Results

Experiments evaluate a three-dimensional Maxwell DG solver on an Nvidia GTX 280 for accuracy, performance, component behavior, implementation choices, and problem-size effects. The solver shows expected high-order convergence and substantial GPU acceleration, while performance depends on operator order, optimization strategy, and workload size.

  • Accuracy: Single-precision L2 errors exhibit the expected asymptotic convergence rate h^N+1 before saturation at the limits of precision.Accuracy is measured in a rectangular perfectly conducting vacuum cavity across orders one through nine and deliberately coarse meshes.
  • Overall performance: Speedups range from 24 to 57 over a single-core 3 GHz Intel Core2 Duo E8400 CPU, reaching 48× at order three and 57× at order four.The comparison counts floating-point additions and multiplications and uses single precision for both implementations.
  • Overall performance: Orders three and four are the fastest GPU methods per degree of freedom because their increased floating-point work offsets hardware-granularity effects.Orders one and two achieve lower overall throughput despite lower computational load, while orders three and four also have moderate timestep requirements.
  • Operator components: Element-local differentiation and lifting scale similarly with order and attain the greatest gains, while flux gather peaks at orders three and four.The local operations have the highest arithmetic intensity and most regular access patterns; increasing order eventually makes them dominate runtime.
  • Memory bandwidth: Flux lifting exceeds the published 141.7 GB/s peak at orders five and above, whereas differentiation approaches peak bandwidth and operator assembly does not reach it.The apparent excess for flux lifting is attributed to texture-cache reuse, while assembly remains dominated by global-memory fetches and stores.
  • Implementation choices: Empirical tuning is necessary: poor parameter choices can make computation take about twice as long, and optimization tradeoffs change when implementation tricks are omitted.Field-in-shared performs best for flux lifting and higher-order differentiation; matrix-in-shared is mainly worthwhile at selected low orders, while microblocking helps more broadly.
  • Problem size: Relatively small problems still achieve decent performance, but sufficient floating-point work per timestep is needed to occupy the GPU effectively.Scaling depends on both element count and order N, because order changes the number of flops per degree of freedom.

7. Conclusions

The paper evaluates GPU techniques for discontinuous Galerkin simulations and combines insights from multiple implementations into a final code. The resulting strategies achieve substantial performance and expand the feasible size and complexity of simulations, while future work targets additional precision and problem classes.

  • Implementation and evaluation: The final code combines insights from adapted DG implementations and applies the strategies developed in Sections 4 and 5.It was used to obtain the computational results reported in Section 6.
  • Performance: DG computations on GPUs achieve speed-up factors just short of two orders of magnitude.The strategies also reach double-digit percentages of the hardware's published theoretical peak performance.
  • Practical impact: The speed increase enables substantially larger and more complex simulations within a given hardware budget.The paper contrasts work formerly requiring a roomful of hardware with simulations feasible on a single device or cluster.
  • Practical guidance: The authors provide implementation advice and performance data to help practitioners balance computing performance against implementation effort.The data also support predicting the computational speed of implementations.
  • Future work: Future work extends the approach to double precision and explores nonlinear conservation laws and elliptic problems.These directions broaden hardware support and the range of targeted equations.

A. Index of Notation

The notation index defines symbols for approximation order, element and matrix quantities, derivatives, GPU scheduling, and workload organization. It also specifies interval notation, rounding, and the principal DG and execution parameters used throughout the paper.

  • Notation conventions: ⌈x⌉n denotes x rounded up to the nearest multiple of n, and [a, b⟩ denotes the integers in the half-open interval [a, b).The notation appears in definitions such as nM = ⌈K/KM⌉.
  • DG discretization: N denotes polynomial degree, while Np and Nfp denote local expansion points and facial nodes, respectively.Nf denotes the number of faces in the reference element.
  • Elements and mappings: K is the total number of elements, Dk the kth finite element, and Ψk the local-to-global map for element k.I denotes the unit finite element.
  • Matrices and derivatives: M, MA, L, S∂µ, and D∂µ denote reference mass, face mass, lifting, stiffness, and differentiation matrices.The index ν denotes global derivatives, while µ denotes local derivatives.
  • GPU workload organization: KM is the number of elements in a microblock, and NpM and NfM are its padded volume and face degrees of freedom.nM denotes the total number of microblocks, while MB counts microblocks in a flux-gather block.
  • GPU execution: tx, ty, tz index threads within a block, bx and by index blocks within the execution grid, and T denotes warp scheduling granularity.wp, wi, and ws distinguish parallel, inline, and sequential work units handled by a block.
Loading 0901.1024v3…